轉(zhuǎn)錄組學(xué)習(xí)一（軟件安裝）
轉(zhuǎn)錄組學(xué)習(xí)二（數(shù)據(jù)下載）
轉(zhuǎn)錄組學(xué)習(xí)三（數(shù)據(jù)質(zhì)控）
轉(zhuǎn)錄組學(xué)習(xí)四（參考基因組及gtf注釋探究）
轉(zhuǎn)錄組學(xué)習(xí)五（reads的比對(duì)與samtools排序）
轉(zhuǎn)錄組學(xué)習(xí)六（reads計(jì)數(shù)與標(biāo)準(zhǔn)化）
轉(zhuǎn)錄組學(xué)習(xí)七（差異基因分析）
轉(zhuǎn)錄組學(xué)習(xí)八（功能富集分析）

任務(wù)

• 學(xué)習(xí)了解各個(gè)reads計(jì)數(shù)轮纫，及標(biāo)準(zhǔn)化的原理，如RPKM/FPKM/TPM的統(tǒng)計(jì)學(xué)原理焚鲜；
• 了解reads計(jì)數(shù)的各個(gè)軟件掌唾，用入門(mén)的htseq-count軟件對(duì)每個(gè)樣本內(nèi)生成關(guān)于表達(dá)量的文件；
• 用腳本合并所有的樣本表達(dá)量文件為表達(dá)矩陣忿磅；
• 對(duì)表達(dá)矩陣在R軟件里簡(jiǎn)單的摸索糯彬，例如求平均值，方差等葱她；
• 看一些重要的生物學(xué)意義的特殊基因的表達(dá)情況如何撩扒，比如GAPDH,β-ACTIN等。

<font color=orange>reads計(jì)數(shù)的原理</font>

參考徐洲更在reads計(jì)數(shù)里的解釋?zhuān)?br> reads的計(jì)數(shù)定量主要可分為三個(gè)水平：基因水平吨些、轉(zhuǎn)錄組水平搓谆、外顯子水平。

• <font color=brown>基因水平</font>：常用的軟件包括HTSeq-count, featureCounts, BEDTools, Qualimap, Rsubread, GenomicRange等豪墅。在轉(zhuǎn)錄組學(xué)習(xí)五中的那篇重要的文章對(duì)于有參的expression定量所推薦的是FeatureCounts泉手，無(wú)參的不組裝直接定量的是Salmon-SMEM。這里對(duì)于入門(mén)所用的HtSeq-count但校，解決的就是上一步mapping完得到的比對(duì)結(jié)果Sam/Bam文件里的read和有參基因組的overlap區(qū)域來(lái)判斷這些reads到底是屬于哪個(gè)基因的螃诅，然后再進(jìn)一步對(duì)這些屬于這個(gè)基因的reads進(jìn)行計(jì)數(shù)，而出現(xiàn)了overlap區(qū)域就會(huì)有不同的情況状囱，對(duì)于HTSeq的三種對(duì)reads計(jì)數(shù)處理的模式如圖：（大多數(shù)情況下术裸，作者推薦union模式）
  

  image
  
  使用完之后，就會(huì)對(duì)每個(gè)樣本輸出一個(gè)表達(dá)量的文件亭枷，后續(xù)需要對(duì)表達(dá)量文件寫(xiě)腳本合并袭艺。

• <font color=brown>轉(zhuǎn)錄本水平</font>：使用的工具為Cufflinks和StringTie，eXpress叨粘。不同的轉(zhuǎn)錄本之間通常是重疊區(qū)域十分多且相似的猾编，當(dāng)二代讀長(zhǎng)小于轉(zhuǎn)錄本長(zhǎng)度時(shí)，就會(huì)有區(qū)分上的難題升敲。不過(guò)現(xiàn)在有三代測(cè)序答倡，能夠完整的將每個(gè)轉(zhuǎn)錄本測(cè)得序列。

• <font color=brown>外顯子水平</font>：和基因水平定量類(lèi)似驴党，不過(guò)需要提供無(wú)重疊的外顯子區(qū)域的gtf文件(?),分析差異外顯子使用的DEXSeq提供了一個(gè)Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個(gè)任務(wù)

• <font color=brown>alignment-free軟件</font>：省去比對(duì)瘪撇，直接得到read count。運(yùn)行效率更高，但會(huì)有樣本特異性和讀長(zhǎng)偏差（不比對(duì)如何知道這些reads是屬于基因組上的哪個(gè)gene上倔既，進(jìn)而如何知道定量計(jì)數(shù)呢恕曲？后續(xù)可以關(guān)注下其算法相關(guān)。）

<font color=orange>標(biāo)準(zhǔn)化的原理RPKM/FPKM/TPM</font>

之前學(xué)習(xí)genek-轉(zhuǎn)錄組原理篇課程時(shí)記錄的相關(guān)筆記渤涌，也再次回顧一下

image

1. 需要標(biāo)準(zhǔn)化的原因：
   1. 樣本內(nèi)：相對(duì)定量而不是絕對(duì)定量：僅表示在35次抽樣中佩谣，B基因抽樣到了20次，（而不是其表達(dá)了20次实蓬。也不能說(shuō)其會(huì)比geneA表達(dá)量高茸俭，因?yàn)榛虻拈L(zhǎng)度不同，所以落到基因上的reads數(shù)量也不同）
   2. 樣本件：不同樣本的測(cè)序深度也是不同的瞳秽，測(cè)序深度更深的會(huì)得到更多的reads瓣履。
   3. 故，需要標(biāo)準(zhǔn)化练俐。對(duì)基因長(zhǎng)度和測(cè)序深度的不同進(jìn)行標(biāo)準(zhǔn)化袖迎。
2. PKM 流程

   1. gene長(zhǎng)度標(biāo)準(zhǔn)化(真實(shí)轉(zhuǎn)錄序列的長(zhǎng)度，不考慮可變剪接情況)
   2. 測(cè)序深度標(biāo)準(zhǔn)化(除以每個(gè)樣品總的reads數(shù)目腺晾，能夠比對(duì)到參考序列上的reads數(shù)目):
3. RPKM(Reads Per Kilobase Per Million.)就是上面的對(duì)1 million個(gè)reads進(jìn)行的操作
   1. kilobase 基因長(zhǎng)度單位
   2. million 測(cè)序深度的單位
4. FPKM(Fragments Per Kilobase Per Million).建庫(kù)測(cè)序是以片段為單位燕锥。單端測(cè)序兩者等價(jià)。而雙端測(cè)序應(yīng)該用FPKM是更為通用的表示方式
5. TPM：
   1. 長(zhǎng)度標(biāo)準(zhǔn)化悯蝉、與FPKM同归形。
   2. 測(cè)序深度的標(biāo)準(zhǔn)化，(不是如FPKM的除以總reads樹(shù))而是按照長(zhǎng)度標(biāo)準(zhǔn)化之后的樣本求和鼻由。除以求和的結(jié)果暇榴。TPM是相等的
6. gene的表達(dá)量到底代表是什么意思。
   • 絕對(duì)定量：一個(gè)細(xì)胞中蕉世，一定mol的RNA種有多少轉(zhuǎn)錄本
   • 但是建庫(kù)測(cè)序時(shí)并不知道用于建庫(kù)測(cè)序時(shí)候提取了多少個(gè)細(xì)胞或者是總共有多少條轉(zhuǎn)錄本蔼紧。所以只能進(jìn)行相對(duì)定量。
   • 相對(duì)定量：某一個(gè)基因的所有轉(zhuǎn)錄本所占的比例狠轻。
7. 根據(jù)公式：
   image
   • TPM:總的轉(zhuǎn)錄本的豐度奸例，更符合對(duì)相對(duì)表達(dá)量的定義。表達(dá)豐度求和向楼。(reads數(shù)目/基因的長(zhǎng)度=表達(dá)豐富度)
   • 看了許多的文章都在強(qiáng)調(diào)RPKM/FPKM的不合理性跟TPM的相對(duì)合理性查吊，但為什么現(xiàn)在大多數(shù)的相對(duì)計(jì)數(shù)都在用FPKM/RPKM，而不用TPM呢湖蜕，因?yàn)楹罄m(xù)的軟件不支持逻卖。

<font color=orange>HTSeq-count定量</font>

根據(jù)featureCounts or htseq-count?文章里對(duì)HTSeq-count和feature-counts做了比較，瀏覽下來(lái)得到幾點(diǎn):HTseq-count為Python寫(xiě)的軟件昭抒，為歐洲分子實(shí)驗(yàn)室大神（寫(xiě)了另外的DESeq2）所寫(xiě)箭阶，并且htseq-count也被更多人所引用虚茶。，而對(duì)于featureCounts來(lái)說(shuō)仇参，速度更快，20倍快于htseq婆殿；支持SAF的注釋文件诈乒，會(huì)得到比htseq-count更多的計(jì)數(shù)值

• 基本使用：

htseq-count [options] <alignment_file> <gff_file>

• 文獻(xiàn)中測(cè)序的數(shù)據(jù)是雙末端PE-reads，htseq的計(jì)數(shù)需要進(jìn)行按照reads名稱(chēng)進(jìn)行排序

### samtools重新排序
for i in `seq 56 58`; do nohup samtools sort -@ 5 -n ../5_samtools/SRR35899${i}.bam -o SRR35899${i}_nsort.bam & done

• HTSeq-count的基本參數(shù)：
  • -f bam/sam:指定文件格式婆芦，默認(rèn)為sam怕磨，選擇為bam。
  • -r name/pos:利用samtool sort對(duì)數(shù)據(jù)根據(jù)read name或者位置進(jìn)行排序消约，默認(rèn)是name肠鲫。
  • -s yes/no/reverse：數(shù)據(jù)是否來(lái)自于strand-specific assay。DNA是雙鏈的或粮，所以需要判斷到底來(lái)自于哪條鏈导饲。如果選擇了no， 那么每一條read都會(huì)跟正義鏈和反義鏈進(jìn)行比較氯材。默認(rèn)的yes對(duì)于雙端測(cè)序表示第一個(gè)read都在同一個(gè)鏈上渣锦，第二個(gè)read則在另一條鏈上。
  • -a [num]: 最低質(zhì)量氢哮， 剔除低于閾值的read
  • -m 模式：union,intersection-strict,intersection-nonempty對(duì)應(yīng)著上圖的三種模式袋毙。
  • **-i **：GFF attribute to be used as feature ID (default,suitable for Ensembl GTF files: gene_id)
• HTSEQ-count腳本：

for i in `seq 56 58`
do
htseq-count -s no -r name -f bam SRR35899${i}_nsort.bam ../../Database/human/Homo_sapiens.GRCh38.87.chr.gtf > SRR35899${i}_matrix.count 2> SRR35899${i}.log
done

• 奇怪的錯(cuò)誤：

  image
  在跑的時(shí)候總是存在libbz2.so.1.0不存在，奇奇怪怪的linux文庫(kù)問(wèn)題冗尤，在之前安裝samtools=1.6版本時(shí)候也遇到這樣情況听盖。此問(wèn)題有待解決。特么花了兩個(gè)小時(shí)的時(shí)間就在處理這個(gè)問(wèn)題裂七，最后發(fā)現(xiàn)可能是conda軟件在安裝pysam跟htseq時(shí)候會(huì)有庫(kù)的問(wèn)題皆看，換用pip install HTSeq即可解決問(wèn)題。

• 運(yùn)行結(jié)果：

  
  
  image

<font color=orange>featureCounts定量</font>

• 基本使用Usage：

featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ...

### 參考其他命令行
featureCounts -T 4 -t CDS -s -g Name -C -a gtf -o readscount_cdf_out 1.sam 3.sam 5.sam ...


$featureCounts -T 5 -p -t exon -g gene_id -a $mm10_gtf -o  counts.txt   *.bam

<font color=orange>腳本合并各個(gè)表達(dá)量的文件生成為表達(dá)矩陣</font>

• Perl里利用哈希數(shù)組來(lái)存入數(shù)據(jù)碍讯，其中基本知識(shí)點(diǎn)邏輯：
  • 匹配.count文件的文件名悬蔽，作為后續(xù)合并的header，
  • 將每個(gè)基因名做key捉兴，value為一維數(shù)組蝎困，n->[num];
  • 每個(gè)hash依次輸出,如果不存在則復(fù)制為"0"
  • foreach #ARGV-1)
  • line[0]}[line[1];
  • print hash{$id}}),"\n"

#!/usr/bin/perl -w

$header= "Gene";
foreach $i ( 0..@ARGV-1){
####    $ARGV[$i]=~/(.*?)_(.*?).txt/;
    $ARGV[$i]=~/(.*?)_matrix\.count/;
    $header.="\t$1";

    open IN,"$ARGV[$i]" or die "$!";
    while(<IN>){
        chomp;
        @line=split;
        $hash{$line[0]}[$i]=$line[1];
    }
}
print"$header\n";
foreach $gene_id(sort keys %hash){
    foreach $i(0..@ARGV-1){
        if(!$hash{$gene_id}[$i]){
            $hash{$gene_id}[$i]=0;
        }
    }
    print $gene_id,"\t",join("\t",@{$hash{$gene_id}}),"\n";
}

• 結(jié)果：

  
  
  image

<font color=orange>重新對(duì)小鼠的4個(gè)數(shù)據(jù)進(jìn)行比對(duì)-計(jì)數(shù)</font>

### 下載小鼠的hisat2 參考基因組index數(shù)據(jù)
wget -b ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grcm38.tar.gz

### 下載小鼠的gtf注釋數(shù)據(jù)：
wget -b wget ftp://ftp.ensembl.org/pub/release-90/gtf/mus_musculus/Mus_musculus.GRCm38.90.gtf.gz

### 對(duì)小鼠的數(shù)據(jù)進(jìn)行比對(duì)，直接輸出為bam文件倍啥，跳過(guò)輸出的sam文件禾乘。
hisat2 -p 5 -x /sas/supercloud-kong/wangtianpeng/Database/mouse/grcm38/genome -1 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR3589959_1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR3589959_2.fastq.gz 2>/sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR3589959.log | samtools view -Sb - > SRR3589959.bam

### HISAT2比對(duì)的循環(huán)腳本
for i in `seq 60 62`
do
nohup hisat2 -x /sas/supercloud-kong/wangtianpeng/Database/mouse/grcm38/genome -1 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR35899${i}_1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR35899${i}_2.fastq.gz 2>/sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR35899${i}.log | samtools view -Sb - > /sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR35899${i}.bam & 
done

### HTSeq-count計(jì)數(shù)
for i in `seq 59 62` 
do 
nohup htseq-count -s no -r name -f bam -i gene_name ../5_samtools/SRR35899${i}.bam /sas/supercloud-kong/wangtianpeng/Database/mouse/Mus_musculus.GRCm38.90.gtf >SRR35899${i}_matrix.count 2> SRR35899${i}.log & 
done

11/22晚11點(diǎn)。Raw_count矩陣導(dǎo)入到R里的探索放在之后的差異基因分析里面吧虽缕∈寂海回顧知識(shí)點(diǎn)，主要是將HTSeq定量的原理，標(biāo)準(zhǔn)化FPKM/RPKM/TPM的原理伍派，HTSeq-count軟件的使用江耀，還有重新將59~62的小鼠的數(shù)據(jù)跑了一遍∷咧玻花了不少時(shí)間祥国，好在也學(xué)到了不少。慢慢向前走吧晾腔。

參考文檔

1. 【使用HTSeq進(jìn)行有參轉(zhuǎn)錄組的表達(dá)量計(jì)算】http://www.chenlianfu.com/?p=2438
2. 【HOPTOP-reads計(jì)數(shù)】https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484502&idx=1&sn=ae6e08bc8ac9374799436f3d301a09ec&chksm=e9e02df7de97a4e14077ca04ac6727ac06012e7d3c3b1ce189bafd679e31b23507078bc0e548&scene=21#wechat_redirect
3. 【轉(zhuǎn)錄組入門(mén)（6）：reads計(jì)數(shù)】http://www.reibang.com/p/24cf44b610a7
4. 【genek-轉(zhuǎn)錄組原理篇】http://www.genek.tv/
5. 【featureCounts or htseq-count?】http://bioinformatics.cvr.ac.uk/blog/featurecounts-or-htseq-count/
6. 【轉(zhuǎn)錄組HTseq對(duì)基因表達(dá)量進(jìn)行計(jì)數(shù)】http://www.bio-info-trainee.com/244.html