豆豆寫(xiě)于2020.2.4
經(jīng)常出現(xiàn)這么一種情況，為了完成目標(biāo)，一定會(huì)先快速做一遍，然后得到最直接的結(jié)果。等后來(lái)就會(huì)發(fā)現(xiàn)，原來(lái)其中還有很多不知道的事情

這次就會(huì)以bowtie2的結(jié)果為例，來(lái)說(shuō)說(shuō)為什么會(huì)有這種感受

首先還是先回顧下bowtie2的基本知識(shí)吧

它的官網(wǎng)在：http://bowtie-bio./bowtie2/manual.shtml

一句話概括這款軟件：

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences.

關(guān)于使用

使用它首先需要對(duì)參考基因組構(gòu)建索引

【注意這里是基因組而不是參考轉(zhuǎn)錄組，因?yàn)樗腂WA一樣，是為基因組測(cè)序而生；如果是要比對(duì)參考轉(zhuǎn)錄組，那么就考慮hisat2吧】

關(guān)于索引有兩種獲得方法：

• 官網(wǎng)下載：提供了人、小鼠、大鼠的索引下載

  

• 自己構(gòu)建：其實(shí)自己下載好參考基因組后，自己構(gòu)建也是很方便的。只需要指定三個(gè)參數(shù)：一個(gè)線程數(shù)，一個(gè)參考基因組的位置，最后就是輸出的索引前綴【一個(gè)小技巧就是輸出的前綴直接用物種的縮寫(xiě)，比如mm10、hg19、ath等】

  # 舉個(gè)例子
ref=/home/genome/genome.fa
bowtie2-build --threads 5 $ref hg19


有了參考索引，就能進(jìn)行單端或者雙端的比對(duì)

可能你會(huì)想，現(xiàn)在都是二代PE測(cè)序了，為什么還有單端的存在呢？但如果接觸了Chipseq數(shù)據(jù)，就能看到大量的單端測(cè)序

# 單端，其中 -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這里一定要寫(xiě)好路徑，并且是寫(xiě)到前綴，于是可以看到index這個(gè)變量就指定到了hg19這里

既然是第二代，那么它的第一代有啥區(qū)別？

之前其實(shí)一直沒(méi)有考慮過(guò)這個(gè)問(wèn)題，因?yàn)榇蠹叶加胋owtie2，而且各種文章也是用新的版本。就是“喜新厭舊”嗎？其實(shí)這些在官網(wǎng)都有說(shuō)過(guò)，只不過(guò)沒(méi)有特殊情況，我們一般也不會(huì)去看，因?yàn)楹芏嗟闹形慕坛套銐蚴褂谩?/p>

第一代版本發(fā)布于2009年，當(dāng)時(shí)二代測(cè)序還沒(méi)有興盛，因此它的目標(biāo)是比對(duì)20-50bp的短reads。但是后來(lái)隨著測(cè)序通量（每天測(cè)序得到的堿基數(shù)更多）和讀長(zhǎng)的增大（每條read能測(cè)到的長(zhǎng)度更長(zhǎng)），bowtie1不能滿足日益發(fā)展的比對(duì)需求，才有了二代。

主要的區(qū)別有：

• 對(duì)于長(zhǎng)度大于50bp的reads，bowtie2的速度更快，更靈敏，使用的內(nèi)存也少；但如果小于50bp，依然是bowtie1更快

• bowtie1只尋找沒(méi)有g(shù)ap的比對(duì)，而bowtie2支持了跨gap的比對(duì)，gap的數(shù)量和長(zhǎng)度都沒(méi)有限制

• bowtie2支持局部比對(duì)（local alignment）和雙端比對(duì)（end-to-end alignment），雙端比對(duì)這一點(diǎn)和bowtie1一樣

• bowtie2對(duì)read長(zhǎng)度沒(méi)有限制，而bowtie1最長(zhǎng)1000bp

• bowtie2允許比對(duì)參考基因組的未知堿基（N），而bowtie1不可以

• 對(duì)于雙端測(cè)序的reads比對(duì)，bowtie2更靈活，因此一會(huì)在結(jié)果解釋中會(huì)看到，雙端測(cè)序的合理比對(duì)和不合理比對(duì)

重點(diǎn)是結(jié)果如何解讀

以往都是關(guān)注最后一行，也就是最后計(jì)算出來(lái)的比對(duì)率，看著差不多也就得了

看看樣子

雙端結(jié)果如下：

單端結(jié)果如下：

需要注意的是，這些具體的比對(duì)信息，bowtie2將他存儲(chǔ)在了標(biāo)準(zhǔn)錯(cuò)誤輸出中，也就是說(shuō)我們需要指定2>align.log這樣來(lái)保存結(jié)果

逐步解讀

就以上面??雙端結(jié)果為例
其中以----為分割線，分成了三大部分

第一行：說(shuō)明一共有12965647對(duì)reads進(jìn)行了比對(duì)【注意這里因?yàn)槭荘E測(cè)序，所以是一對(duì)一對(duì)的，如果論條的話，就應(yīng)該乘以2】

接著，重點(diǎn)關(guān)注一個(gè)名詞：aligned concordantly，直譯就是比對(duì)結(jié)果是不是合理。如果read1和read2都同時(shí)比對(duì)上，并且比對(duì)后的結(jié)果也符合邏輯，那么就是合理的；如果read1和read2能同時(shí)比對(duì)上，但它們各自比對(duì)的位置和本身read1、2之間的間隔差距太大，或者它們比對(duì)的方向壓根就是一樣的，這樣就是不合理（因?yàn)楸旧頊y(cè)序得到的read1和read2應(yīng)該比對(duì)到不同的鏈）

分割的第一部分：在共有12965647對(duì)reads中

• 3805028 (29.35%)對(duì)reads沒(méi)有合理的比對(duì)

• 1212421 (9.35%)對(duì)reads合理比對(duì)了一次

• 7948198 (61.30%)對(duì)reads合理比對(duì)了多次

分割的第二部分：在沒(méi)有合理比對(duì)的3805028對(duì)reads中

• 535030 (14.06%)對(duì)reads 是雙端比對(duì)但不合理

因此這里要注意了：不合理比對(duì)不代表這對(duì)reads不重要，它包含了三種情況：都比對(duì)上但比對(duì)不合理（就是這里的535030）、兩個(gè)read都沒(méi)比對(duì)上、只有一個(gè)read比對(duì)上

分割的第三部分：在不合理比對(duì)的3805028對(duì)reads中，除去雙端比對(duì)但不合理的535030對(duì)，還剩余6539996條（這里沒(méi)有寫(xiě)3269998對(duì)，是因?yàn)檫@些read可能就不是配對(duì)的了），可以看圖中第三部分使用的詞語(yǔ)是mates而不是paired，mates在這里指的就是由不同read組成的群體

• 4021442 (61.49%) 條一次沒(méi)比對(duì)上

• 1314280 (20.10%) 條比對(duì)上一次

• 1204274 (18.41%) 條比對(duì)上多次

最后計(jì)算：84.49% overall alignment rate

如何計(jì)算?
其實(shí)也就是統(tǒng)計(jì)了比對(duì)上的條數(shù)：

雙端比對(duì)的結(jié)果（比對(duì)1次及多次，注意要乘以2）+ 單端比對(duì)的結(jié)果（比對(duì)1次及多次）

補(bǔ)充

• bowtie2默認(rèn)會(huì)同時(shí)尋找concordant和discordant的比對(duì)，除非設(shè)置--no-discordant

• 和bowtie2類(lèi)似，舉一反三，hisat2的結(jié)果也是這么理解