BWT比对算法
简介
BWT算法在多款序列比对软件(BWA,bowtie)中都有涉及,那么对于RNA-seq的2代数据,一般建库长度是单端300bp,双端各150bp左右。
序列比对
对于两个序列进行比对,即pairwise alignment,我们可以按比对方式分为全局比对(NW算法)和局部比对(SW算法):
![](https://img.haomeiwen.com/i19396348/8f5c7924a6b4c03d.png)
当然对于两条短序列,可采用上述算法进行比较,但是如果其中一条序列换成了较长的参考基因组序列,而另一条为fq文件的reads序列,采取这样的方式比对就会显得很慢
![](https://img.haomeiwen.com/i19396348/cad7b3557916f4a6.png)
因此,为了提高效率,我们首先需要根据参考基因组建立index,然后在Seq 2里面取出一个短序列作为seed,通过seed于ref 的index建立联系,再通过ref index找到附近的序列,然后进行pairwise alignment
那么,建立ref index是通过说明方法建立的呢?
1. 第一代软件
比方说SOAP:
![](https://img.haomeiwen.com/i19396348/e5bcb583dad91721.png)
它把参考基因组分成一个个小片段,并且将这些小片段建立成hash关系,也就是说我在序列比对的时候,只要我获取到reads的内容,我就可以match(翻译为完全比对?)到参考基因组上,并获得其位置信息。
2.第二代软件
这里介绍的是BWA,bowtie。而这两款软件除了要对基因组建立index以外,还采用的是BWT算法来实现比对。
(1).什么是BWT算法
![](https://img.haomeiwen.com/i19396348/e5c6865145cd8652.png)
step1:假设我有一段序列: ACAACG
那么我们在序列后面加一个“$”符号:
![](https://img.haomeiwen.com/i19396348/9867f4f60850c02f.png)
step2:
之后,“$”符号向后平移一个单位,得到 :
![](https://img.haomeiwen.com/i19396348/feeb69d7370ac2d5.png)
“$”符号再向后平移一个单位:得到:
![](https://img.haomeiwen.com/i19396348/10ee357160cce5e9.png)
以此类推,最终得到:
![](https://img.haomeiwen.com/i19396348/6f38e8b4dbd1079e.png)
这样的序列的排序表,接下来按照字母顺序进行排序($,A,C,G,T)有:
![](https://img.haomeiwen.com/i19396348/a2f975a205fa7065.png)
这个矩阵我们称之为转换矩阵
我们根据F列和L列元素的相对位置就可以推导出该序列的全部信息:
首先,我们单独取出F列和L列
![](https://img.haomeiwen.com/i19396348/e4031096b50b4372.png)
在F列和L列中找到“$”,连接它们,那么在水平位置对应于的G(蓝色箭头)
这个G就是“$”前面的元素:
![](https://img.haomeiwen.com/i19396348/7e590d381283e6d5.png)
接下来在F列中找到元素G:
![](https://img.haomeiwen.com/i19396348/9f6149c6b9d2b340.png)
此时F列的G对应L列的C(蓝色箭头),此时序列为:
![](https://img.haomeiwen.com/i19396348/8d3a3e1e8b28c53a.png)
此时按照上述方法,才F列中找到元素C,由于L列中的元素C为第二个C(上面G后面的C视为第一个C),所以应该对于F列的第二个C
![](https://img.haomeiwen.com/i19396348/78b3960f36d837b6.png)
那么此时在L列里面对应元素A:
![](https://img.haomeiwen.com/i19396348/e7afdb6be45083a9.png)
接下来L列第三个A对应F列第三个A,因此以此类推就可以最终还原(最终对应于L列的$):
![](https://img.haomeiwen.com/i19396348/a3167444e4d60f96.png)
![](https://img.haomeiwen.com/i19396348/79771a4fa4616ca4.png)
(2).如何实现BWT比对
比方说我现在有一条短序列CAA要比对到ACAACG上,首先说明的是在比对过程中,CAA是反向进行比对的,即A->A->C
由于CAA第一个比对的元素是A,在F列中有三个A,因此我们分类讨论
选择第一个A:
![](https://img.haomeiwen.com/i19396348/f295d7b07b0a1af6.png)
此时比对结果为ACA,不符合
选择第二个A:
![](https://img.haomeiwen.com/i19396348/ea7705802b00d399.png)
结果为$,不符合
选择第三个A:
![](https://img.haomeiwen.com/i19396348/86eb7ecbacddad17.png)
此时比对结果为CAA,符合原序列
在现实中,CAA可以想象为fq文件reads,ACAACG可以想象为参考基因组,并且在比对中会存在mismatch和gap这种情况的发生,因此软件会对每一种比对的情况进行打分,择优处理
目前 bowtie 是没有gap这一项设定的
参考:https://www.bilibili.com/video/av15743137/
刘小乐哈佛大学课程:https://www.bilibili.com/video/BV1De411s71p
https://www.bilibili.com/video/BV1d4411E7uS?from=search&seid=7482313674070699968