BWT算法原理解读

  • A+
所属分类:生物信息学

测序数据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的计算公式如下:

      LF(pos,na)=C(na)+pre\_na(pos,na)

上述公式中,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矩阵直接看出来。为了计算更长的字符串,给出的转换公式如下:

      \\L(pW)=LF(L(W),p) \\U(pW) = LF (U(W)+1, p)-1

如何理解这个公式呢?

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

      \\L(pW)=LF(L(W),p) \\U(pW) = LF (U(W), p)

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

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

  1. L(W)\leq pos\leq U(W) 。
  2. 如果第L(W)行的字符串的末位不是p,则必有:pos\_min>L(W) 。
  3. pre\_na(pos\_min, p)=pre\_na(L(W), p)
  4. 当字符串 pW 存在且唯一时,必然有 pos\_min=pos\_max 。

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

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

同理可以得到:LF(pos\_max, p)=LF(U(W), p) 。

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

但是,注意到一个情况,当 pW 不存在于BWT矩阵中时,理论上是不存在L值和U值的,但是根据这个公式,会计算出L值和U值都等于 LF(L(W),p) [ 因为此时 pre\_na(L(W), p)=pre\_na(U(W), p) ],这显然是有问题的。因此对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算法采用了循环构造的方式,是否会产生这种假阳性:匹配到的目标子串的起点是母串的末位,目标子串的第二个字符串是母串的第一位?