high-throughput seq analysis急先鋒——HTSeq的使用介紹
在所有物是人非的景色里,我最中意你储藐。
正體
HTSeq作為一款可以處理高通量數(shù)據(jù)的python包乙漓,由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人攜手推出HTSeq — A Python framework to work with high-throughput sequencing data吕粗。自發(fā)布以來(lái)就備受廣大分析人員青睞熔掺,其提供了許多功能給那些熟悉python的大佬們?nèi)プ孕判薷氖褂帽ゲ瑫r(shí)也兼顧著給小白們提供了兩個(gè)可以拿來(lái)可用的可執(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的安裝
# 創(chuàng)建存放文件夾
mkdir ~/biosoft/HTseq && cd ~/biosoft/HTseq
# download并解壓
wget https://pypi.python.org/packages/fd/94/b7c8c1dcb7a3c3dcbde66b8d29583df4fa0059d88cc3592f62d15ef539a2/HTSeq-0.9.1.tar.gz#md5=fc71e021bf284a68f5ac7533a57641ac
tar zxvf HTSeq-0.9.1.tar.gz
cd HTSeq-0.9.1/
#使用python命令安裝肺素,此處注意,install --user參數(shù)最好用上深浮,除非你可以獲取root權(quán)限
python setup.py build
python setup.py install --user
# add bin/ to your PATH
vim .bashrc
PATH=/home/path_to/.local/bin:$PATH
source .bashrc
HTSeq使用注意事項(xiàng)
- HTSeq是對(duì)有參考基因組的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)進(jìn)行表達(dá)量分析的压怠,其輸入文件必須有SAM和GTF文件。
- 一般情況下HTSeq得到的Counts結(jié)果會(huì)用于下一步不同樣品間的基因表達(dá)量差異分析飞苇,而不是一個(gè)樣品內(nèi)部基因的表達(dá)量比較。因此蜗顽,HTSeq設(shè)置了-a參數(shù)的默認(rèn)值10布卡,來(lái)忽略掉比對(duì)到多個(gè)位置的reads信息,其結(jié)果有利于后續(xù)的差異分析雇盖。
- 輸入的GTF文件中不能包含可變剪接信息忿等,否則HTSeq會(huì)認(rèn)為每個(gè)可變剪接都是單獨(dú)的基因,導(dǎo)致能比對(duì)到多個(gè)可變剪接轉(zhuǎn)錄本上的reads的計(jì)算結(jié)果是ambiguous崔挖,從而不能計(jì)算到基因的count中贸街。即使設(shè)置-i參數(shù)的值為transcript_id,其結(jié)果一樣是不準(zhǔn)確的狸相,只是得到transcripts的表達(dá)量薛匪。
HTSeq的使用
#這里承接的是上游hisat2比對(duì)軟件得到的bam文件,sort by pos, 所以需要重新sort
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
# 命令參數(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ù)說(shuō)明文件的理解山叮,一般情況下雙末端鏈特異性測(cè)序著榴,該參數(shù)的值應(yīng)該選擇reverse(本人暫時(shí)沒(méi)有測(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,則它們來(lái)自同一個(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)見(jiàn)上面對(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
HTSeq的輸出
HTSeq將Count結(jié)果輸出到標(biāo)準(zhǔn)輸出色查,其結(jié)果示例如下:
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 上的軟件薯演,旨在通過(guò)對(duì)于相對(duì)基因變化找出RNA-seq中表達(dá)差異的基因,同時(shí)也可以用作read count的計(jì)數(shù)秧了。
安裝
gfold.V1.1.4.tar.gzdownload解壓后即可使用
使用
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文件包含五列:
#說(shuō)明很詳細(xì)跨扮,這里不再翻譯
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文件示例:
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é)果出入還是較大的帝嗡,此處僅作說(shuō)明,不加以推薦璃氢。
Bedtools :再一個(gè)count matrix的提取工具
bedtools是一個(gè)極其老牌的數(shù)據(jù)處理軟件了哟玷,由猶他大學(xué)一個(gè)實(shí)驗(yàn)室開(kāi)發(fā),我也是看了生信菜鳥(niǎo)團(tuán)Jimmy的一篇文章才知道也可以用來(lái)計(jì)數(shù)的一也。
安裝
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
使用
bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bed file.bed > read.count.txt
# 2017-09-12 update
# 注意巢寡,此處的bed文件需要自己處理得到的,需要四列塘秦,第一列為chrN讼渊,第二列第三列為基因位置,第四列為基因名尊剔。類似于:
chr1 0 10000 ivl1
chr1 10000 20000 ivl2
chr1 20000 30000 ivl3
chr1 30000 40000 ivl4
輸出
參考文獻(xiàn):
http://www.chenlianfu.com/?p=2438
http://www.#edu.cn/~zhanglab/GFOLD/index.html
http://www.bio-info-trainee.com/745.html