高通量測(cè)序數(shù)據(jù)處理學(xué)習(xí)記錄(二):Read Counts的提取

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)
  1. HTSeq是對(duì)有參考基因組的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)進(jìn)行表達(dá)量分析的压怠,其輸入文件必須有SAM和GTF文件。
  2. 一般情況下HTSeq得到的Counts結(jié)果會(huì)用于下一步不同樣品間的基因表達(dá)量差異分析飞苇,而不是一個(gè)樣品內(nèi)部基因的表達(dá)量比較。因此蜗顽,HTSeq設(shè)置了-a參數(shù)的默認(rèn)值10布卡,來(lái)忽略掉比對(duì)到多個(gè)位置的reads信息,其結(jié)果有利于后續(xù)的差異分析雇盖。
  3. 輸入的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
mode

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


日常Bob鎮(zhèn)樓
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個(gè)濱河市菱皆,隨后出現(xiàn)的幾起案子须误,更是在濱河造成了極大的恐慌,老刑警劉巖仇轻,帶你破解...
    沈念sama閱讀 206,013評(píng)論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件京痢,死亡現(xiàn)場(chǎng)離奇詭異,居然都是意外死亡篷店,警方通過(guò)查閱死者的電腦和手機(jī)祭椰,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評(píng)論 2 382
  • 文/潘曉璐 我一進(jìn)店門(mén),熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)疲陕,“玉大人方淤,你說(shuō)我怎么就攤上這事√阊辏” “怎么了携茂?”我有些...
    開(kāi)封第一講書(shū)人閱讀 152,370評(píng)論 0 342
  • 文/不壞的土叔 我叫張陵,是天一觀的道長(zhǎng)诅岩。 經(jīng)常有香客問(wèn)我讳苦,道長(zhǎng),這世上最難降的妖魔是什么吩谦? 我笑而不...
    開(kāi)封第一講書(shū)人閱讀 55,168評(píng)論 1 278
  • 正文 為了忘掉前任鸳谜,我火速辦了婚禮,結(jié)果婚禮上式廷,老公的妹妹穿的比我還像新娘咐扭。我一直安慰自己,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 64,153評(píng)論 5 371
  • 文/花漫 我一把揭開(kāi)白布草描。 她就那樣靜靜地躺著览绿,像睡著了一般。 火紅的嫁衣襯著肌膚如雪穗慕。 梳的紋絲不亂的頭發(fā)上饿敲,一...
    開(kāi)封第一講書(shū)人閱讀 48,954評(píng)論 1 283
  • 那天,我揣著相機(jī)與錄音逛绵,去河邊找鬼怀各。 笑死,一個(gè)胖子當(dāng)著我的面吹牛术浪,可吹牛的內(nèi)容都是我干的瓢对。 我是一名探鬼主播,決...
    沈念sama閱讀 38,271評(píng)論 3 399
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼胰苏,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼硕蛹!你這毒婦竟也來(lái)了?” 一聲冷哼從身側(cè)響起硕并,我...
    開(kāi)封第一講書(shū)人閱讀 36,916評(píng)論 0 259
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤法焰,失蹤者是張志新(化名)和其女友劉穎,沒(méi)想到半個(gè)月后倔毙,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體埃仪,經(jīng)...
    沈念sama閱讀 43,382評(píng)論 1 300
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,877評(píng)論 2 323
  • 正文 我和宋清朗相戀三年陕赃,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了卵蛉。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 37,989評(píng)論 1 333
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡么库,死狀恐怖傻丝,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情廊散,我是刑警寧澤桑滩,帶...
    沈念sama閱讀 33,624評(píng)論 4 322
  • 正文 年R本政府宣布,位于F島的核電站允睹,受9級(jí)特大地震影響运准,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜缭受,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,209評(píng)論 3 307
  • 文/蒙蒙 一胁澳、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧米者,春花似錦韭畸、人聲如沸宇智。這莊子的主人今日做“春日...
    開(kāi)封第一講書(shū)人閱讀 30,199評(píng)論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)随橘。三九已至,卻和暖如春锦庸,著一層夾襖步出監(jiān)牢的瞬間机蔗,已是汗流浹背。 一陣腳步聲響...
    開(kāi)封第一講書(shū)人閱讀 31,418評(píng)論 1 260
  • 我被黑心中介騙來(lái)泰國(guó)打工甘萧, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留萝嘁,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 45,401評(píng)論 2 352
  • 正文 我出身青樓扬卷,卻偏偏與公主長(zhǎng)得像牙言,于是被迫代替她去往敵國(guó)和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子怪得,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,700評(píng)論 2 345

推薦閱讀更多精彩內(nèi)容