HTSeq使用注意事項
簡述:
- HTSeq是轉(zhuǎn)錄組定量分析的軟件烧董,其輸入文件必須有bam(sorted)和GTF文件票顾。
- 一般情況下HTSeq得到的Counts結(jié)果會用于下一步不同樣品間的基因表達量差異分析蛋铆,而不是一個樣品內(nèi)部基因的表達量比較。因此骨田,HTSeq設(shè)置了-a參數(shù)的默認(rèn)值10,來忽略掉比對到多個位置的reads信息长赞,其結(jié)果有利于后續(xù)的差異分析。
- 輸入的GTF文件中不能包含可變剪接信息闽撤,否則HTSeq會認(rèn)為每個可變剪接都是單獨的基因涧卵,導(dǎo)致能比對到多個可變剪接轉(zhuǎn)錄本上的reads的計算結(jié)果是ambiguous,從而不能計算到基因的count中腹尖。即使設(shè)置-i參數(shù)的值為transcript_id柳恐,其結(jié)果一樣是不準(zhǔn)確的,只是得到transcripts的表達量热幔。
- htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m sorted.bam genome.gtf > counts.txt
參數(shù)
-f | --format default: sam 設(shè)置輸入文件的格式乐设,該參數(shù)的值可以是sam或bam。
-r | --order default: name 設(shè)置sam或bam文件的排序方式绎巨,該參數(shù)的值可以是name或pos近尚。前者表示按read名進行排序,后者表示按比對的參考基因組位置進行排序场勤。若測序數(shù)據(jù)是雙末端測序戈锻,當(dāng)輸入sam/bam文件是按pos方式排序的時候,兩端reads的比對結(jié)果在sam/bam文件中一般不是緊鄰的兩行和媳,程序會將reads對的第一個比對結(jié)果放入內(nèi)存格遭,直到讀取到另一端read的比對結(jié)果。因此留瞳,選擇pos可能會導(dǎo)致程序使用較多的內(nèi)存拒迅,它也適合于未排序的sam/bam文件。而pos排序則表示程序認(rèn)為雙末端測序的reads比對結(jié)果在緊鄰的兩行上她倘,也適合于單端測序的比對結(jié)果璧微。很多其它表達量分析軟件要求輸入的sam/bam文件是按pos排序的,但HTSeq推薦使用name排序硬梁,且一般比對軟件的默認(rèn)輸出結(jié)果也是按name進行排序的前硫。
-s | --stranded default: yes 設(shè)置是否是鏈特異性測序。該參數(shù)的值可以是yes,no或reverse荧止。no表示非鏈特異性測序;若是單端測序罩息,yes表示read比對到了基因的正義鏈上;若是雙末端測序瓷炮,yes表示read1比對到了基因正義鏈上,read2比對到基因負義鏈上苍狰;reverse表示雙末端測序情況下與yes值相反的結(jié)果办龄。根據(jù)說明文件的理解淋昭,一般情況下雙末端鏈特異性測序,該參數(shù)的值應(yīng)該選擇reverse英融。
-a | --a default: 10 忽略比對質(zhì)量低于此值的比對結(jié)果歇式。在0.5.4版本以前該參數(shù)默認(rèn)值是0。
-t | --type default: exon 程序會對該指定的feature(gtf/gff文件第三列)進行表達量計算材失,而gtf/gff文件中其它的feature都會被忽略。
-i | --idattr default: gene_id 設(shè)置feature ID是由gtf/gff文件第9列那個標(biāo)簽決定的笼呆;若gtf/gff文件多行具有相同的feature ID旨别,則它們來自同一個feature,程序會計算這些features的表達量之和賦給相應(yīng)的feature ID昼榛。
-m | --mode default: union 設(shè)置表達量計算模式。該參數(shù)的值可以有union, intersection-strict and intersection-nonempty胆屿。這三種模式的選擇請見上面對這3種模式的示意圖偶宫。從圖中可知,對于原核生物纯趋,推薦使用intersection-strict模式吵冒;對于真核生物纯命,推薦使用union模式痹栖。
-o | --samout 輸出一個sam文件,該sam文件的比對結(jié)果中多了一個XF標(biāo)簽疗我,表示該read比對到了某個feature上咆畏。
-q | --quiet 不輸出程序運行的狀態(tài)信息和警告信息吴裤。
分析腳本
cufflinks --no-update-check \
-o Sample_tmp\
-u --max-bundle-frags 2000000 \
-p 6 \
--library-type fr-firststrand \
-b genome/genome.fa -G genome/gene.gtf sorted.bam
usage: htseq-count [options] alignment_file gff_file