htseq-count計(jì)數(shù)的相關(guān)內(nèi)容前面在不同的學(xué)習(xí)階段寫(xiě)過(guò)兩次粗蔚，分別是合并htseq-count的結(jié)果得到count matrix和htseq-count的一個(gè)坑乃沙，其中第二篇中關(guān)于“坑”的總結(jié)我覺(jué)得還是挺用的士败。
1. 基因表達(dá)定量的三個(gè)水平
- 基因
- 轉(zhuǎn)錄本
- 外顯子
基因和外顯子的定義明晰村怪，統(tǒng)計(jì)起來(lái)相較于轉(zhuǎn)錄本簡(jiǎn)單析藕；而由于不同的轉(zhuǎn)錄本往往存在外顯子的重疊慷垮，統(tǒng)計(jì)起來(lái)就比較困難了阿纤【涔啵基因水平的定量常見(jiàn)。
2. 四種不同的reads計(jì)數(shù)思路
- 當(dāng)比對(duì)到有注釋的基因組時(shí)欠拾，基于注釋文件統(tǒng)計(jì)reads
- 基于參考基因組的轉(zhuǎn)錄本組裝時(shí)胰锌，如Cufflinks會(huì)提供注釋文件, 且能夠發(fā)現(xiàn)新的基因和轉(zhuǎn)錄本。這種情況下清蚀，也要結(jié)合軟件給的注釋文件計(jì)數(shù)
- 比對(duì)到轉(zhuǎn)錄本序列可以直接計(jì)數(shù)匕荸，不借助注釋文件
- 重頭組裝出轉(zhuǎn)錄本序列，接下來(lái)同3
3. 哪些因素影響了feature內(nèi)的reads數(shù)
- 測(cè)序深度
- feature長(zhǎng)度
- feature復(fù)雜度
- GC偏好
- 序列特異偏好
常將前兩者考慮在標(biāo)準(zhǔn)化之內(nèi)
4. 關(guān)于HTSeq
4.1 如何處理多比對(duì)reads
HTSeq忽略掉這些多比對(duì)reads
4.2 HTSeq的計(jì)數(shù)模式
default: union
4.3 HTSeq的使用
usage: htseq-count [options] alignment_file gff_file
-f {sam,bam} (default: sam)
-r {pos,name} (default: name)
-s {yes,no,reverse} (default: yes) #此處關(guān)于選項(xiàng)-s為我自己的認(rèn)識(shí)枷邪，不一定對(duì)
#數(shù)據(jù)是否來(lái)源于鏈特異性測(cè)序榛搔，鏈特異性是指在建庫(kù)測(cè)序時(shí)，只測(cè)mRNA反轉(zhuǎn)錄出的cDNA序列东揣，而不測(cè)該cDNA序列反向互補(bǔ)的另一條DNA序列践惑；換句話說(shuō)就是，鏈特異性能更準(zhǔn)確反映出mRNA的序列信息
#我們知道在gff/gtf中第7列是+-信息嘶卧，+表示來(lái)源于參考基因組序列正鏈尔觉，-表示參考基因組序列的反向互補(bǔ)鏈
#sam/bam文件的第2列是flag信息，也可以看出比對(duì)到正鏈還是負(fù)鏈
#stranded=no芥吟，無(wú)鏈特異性侦铜，一條reads通過(guò)flag列知道比對(duì)到+還是-鏈后专甩，不管是不是和gff/gtf相匹配，都算是這個(gè)feature中的
#stranded=yes, 且se測(cè)序钉稍，要和gff/gtf相匹配涤躲，才算是這個(gè)feature中的
#stranded=yes, 且pe測(cè)序，要和gff/gtf相匹配贡未，才算是這個(gè)feature中的
#stranded=reverse种樱，是yes的相反，這時(shí)不是和gff/gtf相匹配了俊卤，而是恰好相反嫩挤，可能源于另一種鏈特異性，只測(cè)cDNA序列反向互補(bǔ)的另一條DNA序列
-a MINAQUAL (default: 10)
#忽略比對(duì)質(zhì)量低于此值的比對(duì)結(jié)果
-t FEATURETYPE
#feature type (3rd column in GFF file) to be used, all features of other type are ignored (default, suitable for Ensembl GTF files: exon)
#沒(méi)想到這個(gè)還能自己設(shè)置
-i IDATTR
#GFF attribute to be used as feature ID (default, suitable for Ensembl GTF files: gene_id)
-m {union,intersection-strict,intersection-nonempty} (default: union)
4.4 HTSeq輸出結(jié)果
$ ls *count
SRR3286802.count SRR3286803.count SRR3286804.count SRR3286805.count SRR3286806.count SRR3286807.count
#基于相同gff/gtf得到的計(jì)數(shù)文件消恍，行數(shù)相同岂昭，第一列（基因名）相同
$ wc -l *count
37889 SRR3286802.count
37889 SRR3286803.count
37889 SRR3286804.count
37889 SRR3286805.count
37889 SRR3286806.count
37889 SRR3286807.count
#且最后5列統(tǒng)計(jì)了整個(gè)計(jì)數(shù)過(guò)程沒(méi)有使用到的reads
$ tail -n 5 SRR3286802.count
__no_feature 237560
__ambiguous 1846779
__too_low_aQual 0
__not_aligned 1323985
__alignment_not_unique 2015872
- based on the NH tag in the BAM file, they aligned to more than one place in the reference genome (alignment_not_unique);
- they did not align at all (not_aligned);
- their alignment quality was lower than the user-specified threshold (too_low_aQual);
- their alignment overlapped with more than one gene (ambiguous);
- their alignment did not overlap any gene (no_feature).
4.5 什么情況下使用
因?yàn)槭腔趃ff/gtf的feature來(lái)計(jì)數(shù)，所以比對(duì)策略應(yīng)該是往參考基因組上比對(duì)哺哼。
