轉(zhuǎn)錄組分析要用比對(duì)后的bam或sam文件進(jìn)行后續(xù)計(jì)數(shù)。也就是對(duì)基因/轉(zhuǎn)錄本的區(qū)域上比對(duì)到的reads進(jìn)行計(jì)數(shù)。最后算出來某個(gè)轉(zhuǎn)錄本或者外顯子上比對(duì)上多少個(gè)reads订晌。顯然，如果某個(gè)樣本中比對(duì)到某個(gè)基因轉(zhuǎn)錄本上的reads多蚌吸，那這個(gè)基因的表達(dá)量理論上就高锈拨。在不同樣本之間只用計(jì)數(shù)來比對(duì)似乎比較簡單。但考慮到不同基因長度不一樣羹唠，還有管家基因等奕枢，在計(jì)算表達(dá)矩陣的時(shí)候要更復(fù)雜些。
用于計(jì)數(shù)的軟件有：htseq-count佩微，featurecounts(subread)等缝彬。
一、featurecounts
這個(gè)用于計(jì)數(shù)的軟件時(shí)集成在subread里面的哺眯」惹常可以用作RNAseq數(shù)據(jù)計(jì)數(shù)，也可以用于ChIPseq的數(shù)據(jù)計(jì)數(shù)分析哦奶卓。運(yùn)行它需要兩種文件：
1一疯，BAM/SAM作為樣本輸入，里面有每個(gè)reads比對(duì)到染色體上的位置信息寝杖。
2违施，GFF/GTF/SAF文件作為注釋文件，里面有feature 的屬性瑟幕，染色體位置起始點(diǎn)這類信息磕蒲。
注意的是：第二類文件要和前面BAM/SAM獲得用的注釋fa文件的版本來源一致留潦。比如用了UCSC的hg38的fa文件去比對(duì)，那現(xiàn)在就要用UCSC的gtr來注釋feature辣往。（查看day41復(fù)習(xí)）
此外對(duì)feature還要理解一下兔院，feature就是基因組上的一段序列，在一條染色體上的某一片段站削，有明確起點(diǎn)和終點(diǎn)坊萝，有明確的特定功能——也就是它的大小性質(zhì)已被定義。那么mapping到某個(gè)feature上的reads多少许起，就能反應(yīng)樣本的某個(gè)特定功能的強(qiáng)或弱十偶。metafeature就是一些feature組成的集合。比如exon可以說是feature园细，一個(gè)基因所有exon組成了一個(gè)metafeature惦积。
1.單端數(shù)據(jù)
參數(shù)
-g 默認(rèn)為gene_id，就是需要提供一個(gè)id identifier 來將feature水平的統(tǒng)計(jì)匯總為meta-feature水平的統(tǒng)計(jì)猛频。此參數(shù)必須為gtf上有的列名 狮崩。
-t 默認(rèn)為exon，就是read只有落到這些feature上才會(huì)被統(tǒng)計(jì)到鹿寻。此參數(shù)必須為gtf上有的列名 睦柴。可以設(shè)定為-t gene毡熏，就不按照exon了吧坦敌。
-f 設(shè)定后統(tǒng)計(jì)feature層面的數(shù)據(jù)，如exon-level招刹，否則會(huì)統(tǒng)計(jì)meta-feature層面的數(shù)據(jù)
conda activate py2.7 #激活2.7環(huán)境
project=~/rnaseq #項(xiàng)目文件夾
#step5：count
mkdir $project/count
output_dir=$project/count
GTF=/home/cloudam/reference/subread_index/hg38.knownGene.gtf
featureCounts -a $GTF -o $output_dir/counts.txt $project/align/*.bam
用UCSC里面的knownGene.gtf做出來的注釋恬试，基因名完全看不懂。運(yùn)行速度倒是很快疯暑，每個(gè)樣本都不到1min，最后生成的count.txt文件幾十M哑舒，而且是包含了四個(gè)樣本的呢妇拯。
conda activate py2.7 #激活2.7環(huán)境
project=~/rnaseq #項(xiàng)目文件夾
#step5：count
mkdir $project/count
output_dir=$project/count
GTF=/home/cloudam/reference/subread_index/hg38.ncbiRefSeq.gtf
featureCounts -a $GTF -o $output_dir/counts2.txt $project/align/*.bam
用ncbiRefSeq.gtf做注釋文件，就能得到熟悉的gene名啦洗鸵。
2.雙端數(shù)據(jù)
待續(xù)越锈。。膘滨。甘凭。
