RNA-seq分析htseq-count的使用
RNA-seq分析htseq-count的使用
HTSeq作為一款可以處理高通量數據的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人攜手推出HTSeq — A Python framework to work with high-throughput sequencing data。自發布以來就備受廣大分析人員青睞，其提供了許多功能給那些熟悉python的大佬們去自信修改使用，同時也兼顧著給小白們提供了兩個可以拿來可用的可執行文件 htseq-count（計數） 和 htseq-qa（質量分析）。
這里需要注意的是HTSeq作為read counts的計數軟件，承接的是上游比對軟件對于clean data給出的比對結果即bam文件（由sam文件sort得到），和HTSeq能行使同樣作用的還有類似于GFold，bedtools等軟件，我會在最后做一個基本的結果比對。
附manual
附油管視頻講解
HTSeq的安裝
?HTSeq安裝
?
HTSeq使用注意事項
HTSeq的使用
#這里承接的是上游hisat2比對軟件得到的bam文件，sort by pos, 所以需要重新sort| 1 2 3 | samtools?sort?-n yourfile.bam > yourfile_name.bam ? htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/reference/hisat2_reference/Homo_sapiens.GRCh38.86.chr_patch_hapl_scaff.gtf > counts.txt |
| 1 2 3 4 5 6 7 8 9 10 11 | # 命令參數 -f | --format default: sam 設置輸入文件的格式，該參數的值可以是sam或bam。 -r | --order default: name 設置sam或bam文件的排序方式，該參數的值可以是name或pos。前者表示按read名進行排序，后者表示按比對的參考基因組位置進行排序。若測序數據是雙末端測序，當輸入sam/bam文件是按pos方式排序的時候，兩端reads的比對結果在sam/bam文件中一般不是緊鄰的兩行，程序會將reads對的第一個比對結果放入內存，直到讀取到另一端read的比對結果。因此，選擇pos可能會導致程序使用較多的內存，它也適合于未排序的sam/bam文件。而pos排序則表示程序認為雙末端測序的reads比對結果在緊鄰的兩行上，也適合于單端測序的比對結果。很多其它表達量分析軟件要求輸入的sam/bam文件是按pos排序的，但HTSeq推薦使用name排序，且一般比對軟件的默認輸出結果也是按name進行排序的。 -s | --stranded default: yes 設置是否是鏈特異性測序。該參數的值可以是yes,no或reverse。no表示非鏈特異性測序；若是單端測序，yes表示read比對到了基因的正義鏈上；若是雙末端測序，yes表示read1比對到了基因正義鏈上，read2比對到基因負義鏈上；reverse表示雙末端測序情況下與yes值相反的結果。根據說明文件的理解，一般情況下雙末端鏈特異性測序，該參數的值應該選擇reverse（本人暫時沒有測試該參數）。 -a | --a default: 10 忽略比對質量低于此值的比對結果。在0.5.4版本以前該參數默認值是0。 -t | --type default: exon 程序會對該指定的feature（gtf/gff文件第三列）進行表達量計算，而gtf/gff文件中其它的feature都會被忽略。 -i | --idattr default: gene_id 設置feature ID是由gtf/gff文件第9列那個標簽決定的；若gtf/gff文件多行具有相同的feature ID，則它們來自同一個feature，程序會計算這些features的表達量之和賦給相應的feature ID。 -m | --mode default: union 設置表達量計算模式。該參數的值可以有union, intersection-strict and intersection-nonempty。這三種模式的選擇請見上面對這3種模式的示意圖。從圖中可知，對于原核生物，推薦使用intersection-strict模式；對于真核生物，推薦使用union模式。 -o | --samout 輸出一個sam文件，該sam文件的比對結果中多了一個XF標簽，表示該read比對到了某個feature上。 -q | --quiet 不輸出程序運行的狀態信息和警告信息。 -h | --help 輸出幫助信息。 |
htseq-count 的三種比對模式
union, intersection-strict and intersection-nonempty 對照示意圖可以選擇自己需要的模式
我這里使用intersection_nonempty
?
mode
HTSeq的輸出
HTSeq將Count結果輸出到標準輸出，其結果示例如下：| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | head?counts.txt ENSG00000000003 0 ENSG00000000005 0 ENSG00000000419 1171 ENSG00000000457 563 ENSG00000000460 703 ENSG00000000938 0 ENSG00000000971 1 ENSG00000001036 925 ENSG00000001084 1468 ENSG00000001167 2997 ? tail?count.txt ENSG00000283696 18 ENSG00000283697 0 ENSG00000283698 1 ENSG00000283699 0 ENSG00000283700 0 __no_feature??? 3469791 __ambiguous 630717 __too_low_aQual 1346501 __not_aligned?? 520623 __alignment_not_unique? 2849422 |
GFold：另一個count matrix的提取工具
GFold是一款2012年同濟大學的研究組發表在Bioinformatics 上的軟件，旨在通過對于相對基因變化找出RNA-seq中表達差異的基因，同時也可以用作read count的計數。
安裝
gfold.V1.1.4.tar.gzdownload解壓后即可使用
使用
| 1 2 | gfold count -ann hg19Ref.gtf -tag sample1.sam -o sample1.read_cnt gfold count -ann hg19Ref.gtf -tag sample2.sam -o sample2.read_cnt |
輸出
output文件包含五列：
#說明很詳細，這里不再翻譯| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 | GeneSymbol: For BED file, this is the 4'th column. For GPF file, this is the first column. For GTF format, this corresponds to 'gene_id' if it exists, 'NA' otherwise. ? GeneName: For BED file, this is always 'NA'. For GPF file, this is the 12'th column. For GTF format, this corresponds to 'gene_name' if it exists, 'NA' otherwise. ? Read Count: The number of reads mapped to this gene. ? Gene exon length: The length sum of all the exons of this gene. ? RPKM:（#這里需要注意但是雙端測序技術還未普及，這里未使用FPKM，況且RPKM和FPKM也不是能很好的代表基因表達水平 ） The expression level of this gene in RPKM. |
output文件示例：
| 1 2 3 4 5 6 7 8 9 10 11 | head example.read_cnt ENSG00000000003 TSPAN6? 0?? 4535??? 0 ENSG00000000005 TNMD??? 0?? 1610??? 0 ENSG00000000419 DPM1??? 1588??? 1207??? 27.4411 ENSG00000000457 SCYL3?? 1344??? 6883??? 4.07267 ENSG00000000460 C1orf112??? 1334??? 5967??? 4.66292 ENSG00000000938 FGR 0?? 3474??? 0 ENSG00000000971 CFH 2?? 8145??? 0.0051215 ENSG00000001036 FUCA2?? 1427??? 2793??? 10.6564 ENSG00000001084 GCLC??? 2462??? 8463??? 6.06767 ENSG00000001167 NFYA??? 5123??? 3811??? 28.0378 |
此處使用示例bam文件or sam文件和HTSeq的輸入文件一致，但是結果出入還是較大的，此處僅作說明，不加以推薦。
Bedtools ：再一個count matrix的提取工具
bedtools是一個極其老牌的數據處理軟件了，由猶他大學一個實驗室開發，我也是看了生信菜鳥團Jimmy的一篇文章才知道也可以用來計數的。
安裝
| 1 2 | wget https://github.com/arq5x/bedtools2/releases/download/v2.26.0/bedtools-2.26.0.tar.gz tar?zxvf bedtools-2.26.0.tar.gz |
使用
| 1 | bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bed?file.bed >?read.count.txt |
| 1 2 3 4 5 | # 注意，此處的bed文件需要自己處理得到的，需要四列，第一列為chrN，第二列第三列為基因位置，第四列為基因名。類似于： chr1 0?? 10000?? ivl1 chr1 10000?? 20000?? ivl2 chr1 20000?? 30000?? ivl3 chr1 30000?? 40000?? ivl4 |
輸出
?
?
?
標簽:?linux,?RNA-seq
好文要頂?已關注?收藏該文??
總結
以上是生活随笔為你收集整理的RNA-seq分析htseq-count的使用的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 造车的男人们
- 下一篇： 刘强东在耶鲁北京中心演讲(2015-8)