BWT比对算法

作者: 小潤澤 | 来源:发表于2020-07-07 18:00 被阅读0次

    简介

    BWT算法在多款序列比对软件(BWA,bowtie)中都有涉及,那么对于RNA-seq的2代数据,一般建库长度是单端300bp,双端各150bp左右。

    序列比对

    对于两个序列进行比对,即pairwise alignment,我们可以按比对方式分为全局比对(NW算法)和局部比对(SW算法):


    当然对于两条短序列,可采用上述算法进行比较,但是如果其中一条序列换成了较长的参考基因组序列,而另一条为fq文件的reads序列,采取这样的方式比对就会显得很慢

    因此,为了提高效率,我们首先需要根据参考基因组建立index,然后在Seq 2里面取出一个短序列作为seed,通过seed于ref 的index建立联系,再通过ref index找到附近的序列,然后进行pairwise alignment
    那么,建立ref index是通过说明方法建立的呢?

    1. 第一代软件

    比方说SOAP:



    它把参考基因组分成一个个小片段,并且将这些小片段建立成hash关系,也就是说我在序列比对的时候,只要我获取到reads的内容,我就可以match(翻译为完全比对?)到参考基因组上,并获得其位置信息。

    2.第二代软件

    这里介绍的是BWA,bowtie。而这两款软件除了要对基因组建立index以外,还采用的是BWT算法来实现比对。

    (1).什么是BWT算法


    step1:假设我有一段序列: ACAACG
    那么我们在序列后面加一个“$”符号:

    step2
    之后,“$”符号向后平移一个单位,得到 :

    “$”符号再向后平移一个单位:得到:



    以此类推,最终得到:


    这样的序列的排序表,接下来按照字母顺序进行排序($,A,C,G,T)有:



    这个矩阵我们称之为转换矩阵
    我们根据F列和L列元素的相对位置就可以推导出该序列的全部信息:
    首先,我们单独取出F列和L列


    在F列和L列中找到“$”,连接它们,那么在水平位置对应于的G(蓝色箭头)

    这个G就是“$”前面的元素:



    接下来在F列中找到元素G:



    此时F列的G对应L列的C(蓝色箭头),此时序列为:

    此时按照上述方法,才F列中找到元素C,由于L列中的元素C为第二个C(上面G后面的C视为第一个C),所以应该对于F列的第二个C



    那么此时在L列里面对应元素A:

    接下来L列第三个A对应F列第三个A,因此以此类推就可以最终还原(最终对应于L列的$):



    (2).如何实现BWT比对

    比方说我现在有一条短序列CAA要比对到ACAACG上,首先说明的是在比对过程中,CAA是反向进行比对的,即A->A->C
    由于CAA第一个比对的元素是A,在F列中有三个A,因此我们分类讨论
    选择第一个A:



    此时比对结果为ACA,不符合
    选择第二个A:


    结果为$,不符合

    选择第三个A:



    此时比对结果为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

    相关文章

      网友评论

        本文标题:BWT比对算法

        本文链接:https://www.haomeiwen.com/subject/uxuaqktx.html