BWT算法原理解读

测序数据alignment有一些不错的算法，其中Burrows–Wheeler transform算法（简称BWT）是非常高效的一种。本文简单总结下BWT算法思路和原理。

BWT的计算与还原

BWT计算及还原步骤此处不赘述，大致如下图（图1和图2）所示，详细讲解可参考其他资料。

      

            图1：BWT计算（图片来自宾州大学王凯老师）

      

      

            图2：BWT还原（图片来自宾州大学王凯老师）

BWT还原的步骤很巧妙，背后的原理究竟如何？笔者的观点如下。

[ 我们将BWT计算的第四步的矩阵称为BWT矩阵，将BWT计算的输出结果序列称为BWT序列 ]

先看还原的第一步，对BWT序列排序后，本质上得到的是BWT矩阵的第一列，而我们的BWT序列是BWT矩阵的最后一列，又因为矩阵的每一行都是循环生成的，所以第一列和最后一列是毗邻关系，或者说，将第一列结合到最后一列上，矩阵仍然等价（每一行仍然是原始序列的某个循环变换转态），我们可以认为序列关系没被破坏。因此，BWT还原的第一步，将BWT序列加上排序后的序列，等价于将BWT矩阵的第一列移到最后一列后面。此时，若再将得到的列数为2的新序列重新排序，所得到的结构其实就等价于BWT矩阵的前2列，那么再将BWT序列加在前面来构造列数为3的新序列，其实就等价于将BWT矩阵的前2列移到最后一列后面。依次类推，直至构造与原序列相同长度的新序列，也就是将BWT矩阵的前（n-1）列移到最后一列后面。那么，也就是说，这些还原步骤，本质上是利用BWT序列（最后一列信息），构造出了整个BWT矩阵。最后在构造的BWT矩阵中，选择最后一个字符为`$`的就是原序列。

BWT的LF算法

LF算法的作用：对于序列中的某个碱基，如果已知其出现BWT矩阵最后一列中的位置为第pos行，寻找其出现在BWT矩阵第一列中的位置。LF的全称 last to first 也就是指这个意思。

LF的计算公式如下：

      

上述公式中，C(na)可以理解为在BWT矩阵第一列中碱基na第一次出现的位置，pre_na(pos,na)可以理解为在BWT矩阵最后一列中的第pos行之前总共出现过多少个碱基na。

这个公式如何理解呢？

本质上就是先寻找碱基na在第一列中的粗略位置（第一次出现的位置）为C(na)，然后用一个偏移量校正（比如对于碱基A，需要确定是第几个A）。如何求这个偏移量呢？可以证明：如果某个特定位置的碱基在相同碱基中，在最后一列中是第k个出现，则在第一列中也是第k个出现。如：最后一列中第2个出现的C对应于第一列中第2个出现的C。

证明很简单：对于最后一列而言，相同碱基的排序顺序，其实取决于第一列的排序顺序（第一列争气的话，就能往上排）；而对于第一列而言，相同碱基的排序顺序，其实取决于第二列的排序顺序。很巧合，由于是循环序列，这些依赖关系本质上是等价的（对着BWT矩阵多看看哈）。

至此，我们解决了LF的计算，那么如何通过LF算法将BWT序列还原呢？

回顾LF的定义，对于任意碱基如果知道最后一个的位置，就可以求出其在第一列的位置。而我们知道，对于同一行而言，第一个碱基是位于最后一个碱基的后面的。因此可以通过第一个碱基可以反推出其前面的碱基是哪个。如下图所示（图3）：

      

            图3：LF算法还原序列（图片来自宾州大学王凯老师）

BWT搜索字符串

定义 L(W)表示在BWT矩阵中，前缀W出现的最先的位置pos；U(W)表示在BWT矩阵中，前缀W出现的最后的位置pos。

很显然，如果L(W)=U(W)，则搜索到的字符串W的位置唯一；如果L(W)>U(W)，则搜索的字符串W的位置有多个，且连续地从L(W)行至U(W)行都有前缀W（注意：在原序列的位置中并不连续）；如果L(W)<U(W)，显然是矛盾的，说明字符串W不存在（说明预先的假定为错误的）。

当W长度为1时，L(W)和U(W)的值可以从BWT矩阵直接看出来。为了计算更长的字符串，给出的转换公式如下：

      

如何理解这个公式呢？

先看最朴素的想法，转换公式应该写成：

      

设想一下，对于第pos行，如果其前缀为W，而最后一个字符为p，则这一行进行循环变换（将p移到第一列前面），就可以得到pW为前缀的字符串，那么这个循环变换结果在哪一行呢？回想一下LF的作用就知道了，如果能找到这个p在第一列时处于哪一行不就行了嘛！再一次见证了BWT和LF的设计的巧妙所在。

在以W为前缀的字符串中，如果存在第pos行的字符串的末位为p，且令  和  分别表示pos中最小和最大的位置，则有如下几个性质：

1.  。
2. 如果第L(W)行的字符串的末位不是p，则必有： 。
3. 
4. 当字符串 pW 存在且唯一时，必然有  。

性质2成立是因为：在L(W)与pos_min之间，必然不存在字符p。

如果第L(W)行的字符串的末位不是p，则根据性质2及LF的计算公式可知： 。

同理可以得到： 。

因此，上述朴素的公式是合理的。

但是，注意到一个情况，当 pW 不存在于BWT矩阵中时，理论上是不存在L值和U值的，但是根据这个公式，会计算出L值和U值都等于  [ 因为此时  ]，这显然是有问题的。因此对U值的计算引入偏移量-1来抗错。当pW存在时，变形后的公式仍能保证U值计算结果与朴素式的结果相等；当pW 不存在时，变形后的公式计算的U值比朴素的结果小1，这也就使得 U=L-1 ，从而我们能轻易判断出pW是不存在的，因为与我们常识中U>=L的逻辑矛盾！相当于判断了一次 LF(U(W),p) 与 LF(U(W)+1,p)-1是否相等，如果不相等就说明搜索的目标子串不存在。

得到L值和U值后，可以根据还原计算方法，反推出在原序列中的位置。

 

以上便是BWT算法的核心思想。原理的理解，纯属笔者个人观点，如有纰漏，烦请指出。

既然屡清楚了BWT，那么思考简单的问题：由于BWT算法采用了循环构造的方式，是否会产生这种假阳性：匹配到的目标子串的起点是母串的末位，目标子串的第二个字符串是母串的第一位？