RNA-seq分析htseq-count的使用
RNA-seq分析htseq-count的使用
HTSeq作為一款可以處理高通量數(shù)據(jù)的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人攜手推出HTSeq — A Python framework to work with high-throughput sequencing data。自發(fā)布以來就備受廣大分析人員青睞，其提供了許多功能給那些熟悉python的大佬們?nèi)プ孕判薷氖褂?#xff0c;同時(shí)也兼顧著給小白們提供了兩個(gè)可以拿來可用的可執(zhí)行文件 htseq-count（計(jì)數(shù)） 和 htseq-qa（質(zhì)量分析）。
這里需要注意的是HTSeq作為read counts的計(jì)數(shù)軟件，承接的是上游比對(duì)軟件對(duì)于clean data給出的比對(duì)結(jié)果即bam文件（由sam文件sort得到），和HTSeq能行使同樣作用的還有類似于GFold，bedtools等軟件，我會(huì)在最后做一個(gè)基本的結(jié)果比對(duì)。
附manual
附油管視頻講解
HTSeq的安裝
?HTSeq安裝
?
HTSeq使用注意事項(xiàng)
HTSeq的使用
#這里承接的是上游hisat2比對(duì)軟件得到的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 | # 命令參數(shù) -f | --format default: sam 設(shè)置輸入文件的格式，該參數(shù)的值可以是sam或bam。 -r | --order default: name 設(shè)置sam或bam文件的排序方式，該參數(shù)的值可以是name或pos。前者表示按read名進(jìn)行排序，后者表示按比對(duì)的參考基因組位置進(jìn)行排序。若測(cè)序數(shù)據(jù)是雙末端測(cè)序，當(dāng)輸入sam/bam文件是按pos方式排序的時(shí)候，兩端reads的比對(duì)結(jié)果在sam/bam文件中一般不是緊鄰的兩行，程序會(huì)將reads對(duì)的第一個(gè)比對(duì)結(jié)果放入內(nèi)存，直到讀取到另一端read的比對(duì)結(jié)果。因此，選擇pos可能會(huì)導(dǎo)致程序使用較多的內(nèi)存，它也適合于未排序的sam/bam文件。而pos排序則表示程序認(rèn)為雙末端測(cè)序的reads比對(duì)結(jié)果在緊鄰的兩行上，也適合于單端測(cè)序的比對(duì)結(jié)果。很多其它表達(dá)量分析軟件要求輸入的sam/bam文件是按pos排序的，但HTSeq推薦使用name排序，且一般比對(duì)軟件的默認(rèn)輸出結(jié)果也是按name進(jìn)行排序的。 -s | --stranded default: yes 設(shè)置是否是鏈特異性測(cè)序。該參數(shù)的值可以是yes,no或reverse。no表示非鏈特異性測(cè)序；若是單端測(cè)序，yes表示read比對(duì)到了基因的正義鏈上；若是雙末端測(cè)序，yes表示read1比對(duì)到了基因正義鏈上，read2比對(duì)到基因負(fù)義鏈上；reverse表示雙末端測(cè)序情況下與yes值相反的結(jié)果。根據(jù)說明文件的理解，一般情況下雙末端鏈特異性測(cè)序，該參數(shù)的值應(yīng)該選擇reverse（本人暫時(shí)沒有測(cè)試該參數(shù)）。 -a | --a default: 10 忽略比對(duì)質(zhì)量低于此值的比對(duì)結(jié)果。在0.5.4版本以前該參數(shù)默認(rèn)值是0。 -t | --type default: exon 程序會(huì)對(duì)該指定的feature（gtf/gff文件第三列）進(jìn)行表達(dá)量計(jì)算，而gtf/gff文件中其它的feature都會(huì)被忽略。 -i | --idattr default: gene_id 設(shè)置feature ID是由gtf/gff文件第9列那個(gè)標(biāo)簽決定的；若gtf/gff文件多行具有相同的feature ID，則它們來自同一個(gè)feature，程序會(huì)計(jì)算這些features的表達(dá)量之和賦給相應(yīng)的feature ID。 -m | --mode default: union 設(shè)置表達(dá)量計(jì)算模式。該參數(shù)的值可以有union, intersection-strict and intersection-nonempty。這三種模式的選擇請(qǐng)見上面對(duì)這3種模式的示意圖。從圖中可知，對(duì)于原核生物，推薦使用intersection-strict模式；對(duì)于真核生物，推薦使用union模式。 -o | --samout 輸出一個(gè)sam文件，該sam文件的比對(duì)結(jié)果中多了一個(gè)XF標(biāo)簽，表示該read比對(duì)到了某個(gè)feature上。 -q | --quiet 不輸出程序運(yùn)行的狀態(tài)信息和警告信息。 -h | --help 輸出幫助信息。 |
htseq-count 的三種比對(duì)模式
union, intersection-strict and intersection-nonempty 對(duì)照示意圖可以選擇自己需要的模式
我這里使用intersection_nonempty
?
mode
HTSeq的輸出
HTSeq將Count結(jié)果輸出到標(biāo)準(zhǔn)輸出，其結(jié)果示例如下：| 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：另一個(gè)count matrix的提取工具
GFold是一款2012年同濟(jì)大學(xué)的研究組發(fā)表在Bioinformatics 上的軟件，旨在通過對(duì)于相對(duì)基因變化找出RNA-seq中表達(dá)差異的基因，同時(shí)也可以用作read count的計(jì)數(shù)。
安裝
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文件包含五列：
#說明很詳細(xì)，這里不再翻譯| 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:（#這里需要注意但是雙端測(cè)序技術(shù)還未普及，這里未使用FPKM，況且RPKM和FPKM也不是能很好的代表基因表達(dá)水平 ） 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的輸入文件一致，但是結(jié)果出入還是較大的，此處僅作說明，不加以推薦。
Bedtools ：再一個(gè)count matrix的提取工具
bedtools是一個(gè)極其老牌的數(shù)據(jù)處理軟件了，由猶他大學(xué)一個(gè)實(shí)驗(yàn)室開發(fā)，我也是看了生信菜鳥團(tuán)Jimmy的一篇文章才知道也可以用來計(jì)數(shù)的。
安裝
| 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 |
輸出
?
?
?
標(biāo)簽:?linux,?RNA-seq
好文要頂?已關(guān)注?收藏該文??
總結(jié)
以上是生活随笔為你收集整理的RNA-seq分析htseq-count的使用的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 造车的男人们
- 下一篇： 刘强东在耶鲁北京中心演讲(2015-8)