序列比对（四） 启发式比对之BWT算法

由于精确比对在处理数据量较大的序列比对问题时，耗时比较长，所以需要一些更快速的比对算法。比如之前为了解决数据库搜索问题，而产生的BLAST算法。我们今天介绍的BWT算法也是为了解决一个精确比对无法解决的问题，或者说在短时间内无法解决的问题。那就是短序列比对到特长序列上，如二代测序数据比对到人类基因组上。由于人类基因组序列长度很长，如果使用精确比对，一条短序列的比对就会耗时很久。这样就必须需要一种快速的比对方法来解决这个棘手的问题。
目前二代测序常用的软件也是基于BWT算法的，比如BWA，Bowtie等等。BWT算法原本是用于数据压缩的，在数据压缩上效率很高。下面我会根据我的理解，基于BWA软件的文章，简单介绍下BWT算法，包括，BWT转换，以及具体的比对方法。

（一）BWT转换

1.BWT编码

原始序列：

ACCGCGCGT

加入$作为序列尾部字符：

0 ACCGCGCGT$

1 CCGCGCGT$A

2 CGCGCGT$AC

3 GCGCGT$ACC

4 CGCGT$ACCG

5 GCGT$ACCGC

6 CGT$ACCGCG

7 GT$ACCGCGC

8 T$ACCGCGCG

9 $ACCGCGCGT

按字典表排序：

0 9 $ACCGCGCGT

1 0 ACCGCGCGT$

2 1 CCGCGCGT$A

3 2 CGCGCGT$AC

4 4 CGCGT$ACCG

5 6 CGT$ACCGCG

6 3 GCGCGT$ACC

7 5 GCGT$ACCGC

8 7 GT$ACCGCGC

9 8 T$ACCGCGCG

BWT序列，就是排序后的矩阵的最后一列字符。

T$ACGGCCCG

如果对BWT序列进行一定的压缩，如T$AC2G3CG，就可以起到很好的压缩效果。

2. BWT解码

可以通过BWT序列（以下成为L列）和排序后矩阵的第一列序列(以下成为F列)，就可以还原原序列，具体的解码如下：
我们需要从L列的第一个字符开始，由于$字符是在所有字符之后，所有BWT序列的第一个字符一定是原序列的最后一个字符，所以我们由此开始

T –> 找到F列中第一个T的位置即（9），则从（9）中发现在L中的相应序号的位置的字符是G，即为T在原序列的前一个位置的字符

以此类推，就可以还原原序列。

解码伪代码：

def LF(r):
    c = BWT[r]
    return C(c) + Occ(c,r) + 1

def BWT_reverse():
    r = 1
    T = ""
    while (BWT[r] != "$"):
        T = T + BWT[r]
        r = LF(r)
    return T

其中，C[c]的定义是，字典序小于字符c的所有字符个数。Occ[c,r]表示在BWT（T）中第r行之前出现字符c的个数。

（二）LP查询，精确匹配

也就是根据查询序列从参考序列中找到相应的比对位置。大致的思路与上面的解码类似，只是开始搜索的位置需要从F列开始，根据查询序列的最后一个字符在F列中的位置（一般会有多个位置）开始不断往前匹配，从而找到相应的位置。

精确匹配伪代码：

def LF(c,r):
    return C(c) + Occ(c,r) + 1

def exactmatch(Q[1,q]):
    c = Q[q]
    sp = C[c] + 1
    ep = C[c+1] + 1
    i = q - 1
    while (ep > sp) and (i > 0):
        c = Q[i]
        sp = LFC(c,sp)
        ep = LFC(c,ep)
        i = i - 1
    return sp, ep

其中，c+1表示字典序比字符c大的下一个字符。

（三）小结

在具体算法实现过程中，比如软件BWA和Bowtie等，还会考虑内存和效率上的取舍，选择合适的方案，主要优化点在提前计算好一些数据，如C[c],Qcc(c,r)以及相应的位置信息等等，当然考虑到内存消耗也不会完全储存，会储存简化版。也会考虑mismatch和indel等等匹配。