大神一句話,菜鳥跑半年。我不是大神,但我可以縮短你走彎路的半年~
就像歌兒唱的那樣,如果你不知道該往哪兒走,就留在這學點生信好不好~
這里有豆豆和花花的學習歷程,從新手到進階,生信路上有你有我!
豆豆寫于2020.2.4
經(jīng)常出現(xiàn)這么一種情況,為了完成目標,一定會先快速做一遍,然后得到最直接的結(jié)果。等后來就會發(fā)現(xiàn),原來其中還有很多不知道的事情這次就會以bowtie2的結(jié)果為例,來說說為什么會有這種感受
它的官網(wǎng)在:http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
一句話概括這款軟件:
Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences.
【注意這里是基因組而不是參考轉(zhuǎn)錄組,因為它的BWA一樣,是為基因組測序而生;如果是要比對參考轉(zhuǎn)錄組,那么就考慮hisat2吧】
關(guān)于索引有兩種獲得方法:
官網(wǎng)下載:提供了人、小鼠、大鼠的索引下載
自己構(gòu)建:其實自己下載好參考基因組后,自己構(gòu)建也是很方便的。只需要指定三個參數(shù):一個線程數(shù),一個參考基因組的位置,最后就是輸出的索引前綴【一個小技巧就是輸出的前綴直接用物種的縮寫,比如mm10、hg19、ath等】
# 舉個例子
ref=/home/genome/genome.fa
bowtie2-build --threads 5 $ref hg19
可能你會想,現(xiàn)在都是二代PE測序了,為什么還有單端的存在呢?但如果接觸了Chipseq數(shù)據(jù),就能看到大量的單端測序
# 單端,其中 -p指定線程數(shù)
index=/home/genome/hg19
fq=/YOUR_PATH/
bowtie2 -p 5 -x $index -U $fq
# 雙端
index=/home/genome/hg19
fq1=/YOUR_PATH/
fq2=/YOUR_PATH/
bowtie2 -p 5 -x $index -1 $fq1 -2 $fq2
必須注意的是,bowtie2 -x
這里一定要寫好路徑,并且是寫到前綴,于是可以看到index
這個變量就指定到了hg19這里
之前其實一直沒有考慮過這個問題,因為大家都用bowtie2,而且各種文章也是用新的版本。就是“喜新厭舊”嗎?其實這些在官網(wǎng)都有說過,只不過沒有特殊情況,我們一般也不會去看,因為很多的中文教程足夠使用。
第一代版本發(fā)布于2009年,當時二代測序還沒有興盛,因此它的目標是比對20-50bp的短reads。但是后來隨著測序通量(每天測序得到的堿基數(shù)更多)和讀長的增大(每條read能測到的長度更長),bowtie1不能滿足日益發(fā)展的比對需求,才有了二代。
主要的區(qū)別有:
對于長度大于50bp的reads,bowtie2的速度更快,更靈敏,使用的內(nèi)存也少;但如果小于50bp,依然是bowtie1更快
bowtie1只尋找沒有g(shù)ap的比對,而bowtie2支持了跨gap的比對,gap的數(shù)量和長度都沒有限制
bowtie2支持局部比對(local alignment)和雙端比對(end-to-end alignment),雙端比對這一點和bowtie1一樣
bowtie2對read長度沒有限制,而bowtie1最長1000bp
bowtie2允許比對參考基因組的未知堿基(N),而bowtie1不可以
對于雙端測序的reads比對,bowtie2更靈活,因此一會在結(jié)果解釋中會看到,雙端測序的合理比對和不合理比對
以往都是關(guān)注最后一行,也就是最后計算出來的比對率,看著差不多也就得了
雙端結(jié)果如下:
單端結(jié)果如下:
需要注意的是,這些具體的比對信息,bowtie2將他存儲在了標準錯誤輸出中,也就是說我們需要指定2>align.log
這樣來保存結(jié)果
就以上面??雙端結(jié)果為例
其中以----
為分割線,分成了三大部分
第一行:說明一共有12965647對reads進行了比對【注意這里因為是PE測序,所以是一對一對的,如果論條的話,就應(yīng)該乘以2】
接著,重點關(guān)注一個名詞:aligned concordantly
,直譯就是比對結(jié)果是不是合理。如果read1和read2都同時比對上,并且比對后的結(jié)果也符合邏輯,那么就是合理的;如果read1和read2能同時比對上,但它們各自比對的位置和本身read1、2之間的間隔差距太大,或者它們比對的方向壓根就是一樣的,這樣就是不合理(因為本身測序得到的read1和read2應(yīng)該比對到不同的鏈)
分割的第一部分:在共有12965647對reads中
3805028 (29.35%)對reads沒有合理的比對
1212421 (9.35%)對reads合理比對了一次
7948198 (61.30%)對reads合理比對了多次
分割的第二部分:在沒有合理比對的3805028對reads中
535030 (14.06%)對reads 是雙端比對但不合理
因此這里要注意了:不合理比對不代表這對reads不重要,它包含了三種情況:都比對上但比對不合理(就是這里的535030)、兩個read都沒比對上、只有一個read比對上
分割的第三部分:在不合理比對的3805028對reads中,除去雙端比對但不合理的535030對,還剩余6539996條(這里沒有寫3269998對,是因為這些read可能就不是配對的了),可以看圖中第三部分使用的詞語是mates而不是paired,mates在這里指的就是由不同read組成的群體
4021442 (61.49%) 條一次沒比對上
1314280 (20.10%) 條比對上一次
1204274 (18.41%) 條比對上多次
如何計算?
其實也就是統(tǒng)計了比對上的條數(shù):雙端比對的結(jié)果(比對1次及多次,注意要乘以2)+ 單端比對的結(jié)果(比對1次及多次)
# 來自分割的第一部分
1212421 x 2 + 7948198 x 2
# 來自分割的第二部分
535030 x 2
# 來自分割的第三部分
1314280 + 1204274
# 最后總共比對的條數(shù)是
21909852
# 全部條數(shù)是
12965647 x 2 = 25931294
# 最后的比例就是
21909852/25931294*100%=84.49%
bowtie2默認會同時尋找concordant
和discordant
的比對,除非設(shè)置--no-discordant
和bowtie2類似,舉一反三,hisat2的結(jié)果也是這么理解