序列比对（27）BWT算法

原创已于 2024-08-29 11:15:37 修改 · 4.2k 阅读

46 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#bwt

于 2019-10-24 11:03:52 首次发布

序列算法专栏收录该内容

32 篇文章

订阅专栏

本文深入探讨BWT（Burrows-Wheeler Transform）算法，包括编码和解码过程。BWT通过数据转换使得相似字符相邻，常用于生物信息学中的序列比对。文中提供了两种实现代码，并解释了如何通过L列恢复原始字符串。

本文介绍了BWT算法。

bwa是目前最流行的二代测序比对工具，其中就用到了BWT算法。BWT（Burrows-Wheeler Transform）算法是一种数据转换算法，它将一个字符串中的相似字符放在相邻的位置，以便于后续的压缩。

简要回顾

BWT算法可以分为编码和解码两部分。编码后，原始字符串中的相似字符会处在比较相邻的位置；解码就是将编码后的字符串重新恢复成原始字符串的过程。BWT的一个特点就是经过编码后的字符串可以完全恢复成原始字符串。

BWT编码分为以下几步：

输入一个字符串 $s$ ，假设其中所有字符都介于 $a-z\text{a-z}$ 之间。
在 $s$ 的末尾加上一个标记字符，该字符要比 $s$ 中的所有字符都要小。比如 $$\text{\textdollar}$ 字符。这样将末尾加上标记的新字符串记为 $s^{'}$ 。
重复地将 $s^{'}$ 中的最后一个字符转移到开头，每转移一次就得到一个新的字符串。
将上一步得到的所有新字符串从小到大排序，排序后的字符串数组记为 $M$ 。
$M$ 中每个字符串的第一个字符构成 $F$ 列， $M$ 中每个字符串的最后一个字符构成 $L$ 列。
输出 $L$ 列。

在这里插入图片描述

BWT解码分为以下几步：
（至于为什么这样处理就可以恢复原始字符串，下文会说明。）

输入 $L$ 列。
对 $L$ 列进行从小到大排序得到 $F$ 列。
$L$ 列的第一个字符是原始字符串的最后一个字符。
根据 $L$ 列的字符 $L_i$ 找到 $F$ 列中的相同字符 $F_j$ ，然后得到 $F_j$ 所在行的最后一个字符 $L_j$ 。将 $L_j$ 记录下来。
重复上面一步，直到 $F_j$ 等于标记字符为止。
按照上述步骤找到的各个 $L_j$ 进行反向排列，得到字符串 $r$ 。
输出字符串 $r$ 。

用图示表示就是：
Step1：输入 $L$ 列；对 $L$ 列进行从小到大排序得到 $F$ 列。
在这里插入图片描述

Step2：循环解码。
根据 $L$ 列的字符 $L_i$ 找到 $F$ 列中的相同字符 $F_j$ ，然后得到 $F_j$ 所在行的最后一个字符 $L_j$ 。将 $L_j$ 记录下来；重复上面一步，直到 $F_j$ 等于标记字符为止；按照上述步骤找到的各个 $L_j$ 进行反向排列，得到字符串 $r$ 。

循环解码的具体步骤如下：
Step2-1：
$L$ 列的第一个字符 $c$ 是原始字符串的最后一个字符。
在这里插入图片描述

Step2-2：
根据 $L$ 列中第一行 $L_1$ 的 $c$ 这个字符找到 $F$ 列中同样字符所在的行。由于 $F$ 列中只有一个 $c$ ，所以就是最后一行，对应的 $F_6$ 。
根据 $F_6$ 找到同一行的 $L_6$ ，即字符 $b$ 。
在这里插入图片描述

Step2-3：
根据 $L_6$ 的 $b$ 这个字符找到 $F$ 列中同样字符所在的行。由于 $F$ 列中有两个 $b$ ，那选择哪一个呢？
一般地，如果 $Li(i=1,2,…,n)L_i (i=1,2,\ldots,n)$ 这个字符在 $F$ 列中出现过多次，分别是 $Fj,Fj+1,…F_j, F_{j+1}, \ldots$ ；并且假设 $L_1$ 到 $L_i$ 里这 $i$ 个字符中一共有 $k$ 个( $\le i$ )字符等于 $L_i$ 这个字符，那么我们要为 $L_i$ 在 $F$ 列中找的对应的字符就是 $F_{j+k-1}$ 。
比如对应到这一步， $i = 6$ ， $L_6$ 是字符 $b$ ，而 $F$ 列中共有 $2$ 个 $b$ ： $F_4$ 和 $F_5$ ， $L_6$ 对应哪一个呢？我们看到 $L_1$ 到 $L_6$ 中共有 $2$ 个 $b$ ，按照上面说的，选择 $F_{4+2-1}$ ，即 $F_5$ 。
根据 $F_5$ 找到同一行的 $L_5$ ，即字符 $a$ 。
在这里插入图片描述

Step2-4：
根据 $L_5$ 的 $a$ 这个字符找到 $F$ 列中同样字符所在的行。由于 $F$ 列中有两个 $a$ ，那选择哪一个呢？
同样地， $F$ 列中共有 $2$ 个 $a$ ： $F_2$ 和 $F_3$ ；而 $L_1$ 到 $L_5$ 中共有 $2$ 个 $a$ ，按照上面说的，选择 $F_{2+2-1}$ ，即 $F_3$ 。
根据 $F_3$ 找到同一行的 $L_3$ ，即字符 $b$ 。
在这里插入图片描述

Step2-5：
根据 $L_3$ 的 $b$ 这个字符找到 $F$ 列中同样字符所在的行。 $F$ 列中共有 $2$ 个 $b$ ： $F_4$ 和 $F_5$ ；而 $L_1$ 到 $L_3$ 中共有 $1$ 个 $b$ ，按照上面说的，选择 $F_{4+1-1}$ ，即 $F_4$ 。
根据 $F_4$ 找到同一行的 $L_4$ ，即字符 $a$ 。
在这里插入图片描述

Step2-6：
根据 $L_4$ 的 $a$ 这个字符找到 $F$ 列中同样字符所在的行。 $F$ 列中共有 $2$ 个 $a$ ： $F_2$ 和 $F_3$ ；而 $L_1$ 到 $L_4$ 中共有 $1$ 个 $a$ ，按照上面说的，选择 $F_{2+1-1}$ ，即 $F_2$ 。
根据 $F_2$ 找到同一行的 $L_2$ ，即字符 $$\text{\textdollar}$ 。至此，解码结束。
在这里插入图片描述

Step3：输出原始字符串。
在这里插入图片描述

关键步骤

现在我们来说明为什么按照上述解码过程就可以恢复原始字符串。关键就是要解答两个问题：

为什么要在 $F$ 列中找和 $L_i$ 相同的字符？
如果 $F$ 列中有多个字符和 $L_i$ 相同，怎么办？

问题一：为什么要在 $F$ 列中找和 $L_i$ 相同的字符？

我们重新看BWT编码中的“循环转移”这一步。我们将某一行字符串的Latter String定义为其在“循环转移”这一步中的下一行字符串；而将某一行字符串的Former String定义为其在“循环转移”这一步中的上一行字符串。
在这里插入图片描述

我们可以看出，某一行字符串的最后一个字符是其Latter String的第一个字符；某一行字符串的最后一个字符和其Latter String的最后一个字符的关系是：在原始字符串中，上述两个字符紧挨在一起并且Latter String的最后一个字符排在前面。

比如，根据定义，在“循环转移”这一步中，第一行的 Latter String是第二行；在 $M$ 数组中，第一行的 Latter String是第六行。在“循环转移”这一步中，第一行字符串的最后一个字符是 $c$ ，第一行字符串的 Latter String（也就是第二行字符串）的最后一个字符是 $b$ ，二者的关系是：在原始字符串中， $b$ 排在 $c$ 的前面一位。

现在我们可以回答问题了：我们在 $F$ 列中找和 $L_i$ 相同的字符，就可以找到 $L_i$ 所在行的字符串的Latter String。如上所述，在原始字符串中，Latter String的最后一个字符排在 $L_i$ 前面。

比如， $L_1$ 和 $F_6$ 相同，那么 $L_1$ 所在的行（也就是 $M$ 数组的第一行）的 Latter String 就是 $F_6$ 所在的行（也就是 $M$ 数组的第六行）。

并且，我们让Latter String的最后一个字符成为新的 $L_{i^’}$ ，不断重复这个过程，就可以恢复原始字符串。
在这里插入图片描述

问题二：如果 $F$ 列中有多个字符和 $L_i$ 相同，怎么办？

首先我们要说明一个规律： $M$ 数组中，以某字符结尾的一个字符串在以该字符结尾的所有字符串中的相对顺序和该字符串的Latter String在以该字符串开头的所有字符串中的相对顺序是一样的。

比如， $M$ 数组中，以字符 $b$ 结尾的字符串共有两行，第三行和第六行，第六行的相对顺序是 $2$ 。而以字符 $b$ 开头的字符串也是两行，第四行和第五行，第五行的相对顺序是 $2$ 。而第五行正是第六行的 Latter String，二者的相对顺序是一样的。
在这里插入图片描述

因此，在 $M$ 数组中，假设 $L_1$ 到 $L_i$ 这 $i$ 个字符中一共有 $k$ 个( $\le i$ )字符等于 $L_i$ 这个字符，那么 $L_i$ 所在行的字符串在以 $L_i$ 这个字符结尾的所有字符串中的相对顺序就是 $k$ ，从而 $L_i$ 所在行的字符串的Latter String在以 $L_i$ 这个字符开头的所有字符串中的相对顺序也是 $k$ 。假设 $L_i$ 这个字符在 $F$ 列中出现过多次，分别是 $Fj,Fj+1,…F_j, F_{j+1}, \ldots$ ，那么以 $L_i$ 这个字符开头的所有字符串所在的行就是第 $j,j+1,…j,j+1,\ldots$ 行；其中相对顺序为 $k$ 的行是 $j + k - 1$ 行。也就是说Latter String所在的行是第 $j + k - 1$ 行。
在这里插入图片描述
红色方框为所有以b开头的字符串；绿色方块为所有以b结尾的字符串。黑色箭头指向的是对应的Latter String。

比如，在 $M$ 数组中， $L_6=b$ ， $L_1$ 到 $L_6$ 共有 $2$ 个字符是 $b$ ，所以 $L_6$ 所在的第六行字符串在所有以 $b$ 结尾的字符串中的相对顺序是 $2$ 。从而其 Latter String在所有以 $b$ 开头的字符串中的相对顺序也是 $2$ ，也就是第 $5$ 行。

综上所述， $L_i$ 所在行的字符串的Latter String是第 $j + k - 1$ 行。也就是说，当 $F$ 列中有多个字符和 $L_i$ 相同时，我们可以通过计算从 $L_1$ 到 $L_i$ 共有几个字符和 $L_i$ 相同从而知道Latter String是第几行。

实现代码一

代码的实现有几点要说明：

在编码过程中，“循环转移”这一步中会产生 $n$ 个新字符串，对这些字符串进行排序得到 $M$ 数组。为了存储“循环转移”这一步中产生的各个新字符串，空间利用为 $O(n^2)$ 。更具体地，需要占用 $n^2$ 个字节，但是如果我们只存储每个新字符串的一部分，即开头字符到标记字符的那一部分，总共可以节省约一半的空间而不影响后续的排序。
解码的实现有两种方法，都是为了解决如何根据 $L_i$ 找到同字符的 $F_j$ 。本小节的方法是构建一个“跳转数组”，这个方法的巧妙之处是预先在线性时间内计算出每个 $Li(i=1,2,…,n)L_i(i=1,2,\ldots,n)$ 之前有几个字符和它相同，从而减少根据 $L_i$ 查找 $F_j$ 所需的总的计算量。该方法的缺点是“跳转数组”需要额外的空间开销。

效果如下：
在这里插入图片描述

具体的实现代码如下:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSTR 1000
#define MARKER '$'
#define NUM_ALPHA 26

int comp(const void *s, const void *t) {
    return strcmp(*(char**)s, *(char**) t);    /* 注意这里 *(char**) 的用法 */
}

/* the last char of s is not MARKER */
char* bwtEncode(char *s, const int n) {
    char *L;
    char **M;
    int i, j, l;
    if ((M = (char**) malloc(sizeof(char*) * n)) == NULL || \
        (L = (char*) malloc(sizeof(char) * (n + 2))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < n; i++) {
        if ((M[i] = (char*) malloc(sizeof(char) * (i + 2))) == NULL) {   /* 只需保存开头到标记字符的那一部分字符串，这样总共可以节省大约一半的空间 */
            fputs("Error: out of space!\n", stderr);
            exit(1);       
        }
        for (j = 0; j < i + 1; j++)
            M[i][j] = s[n - 1 - i + j];   /* 这里的字符串没有存储 MARKER */
        M[i][i + 1] = '\0';
    }
    qsort(M, n, sizeof(M[0]), comp);   /* 对旋转后的多个字符串排序 */
    for (i = 0, L[0] = s[n - 1]; i < n; i++) {
        if ((l = strlen(M[i])) < n)
            L[i + 1] = s[n - 1 - l];
        else
            L[i + 1] = MARKER;
    }
    L[n + 1] = '\0';
    for (i = 0; i < n; i++)
        free(M[i]);
    free(M);
    return L;
}

char* bwtDecode(char *L, const int n) {
    int i;
    int *a, *b;
    char *r;   /* original string. */
    int pos;
    if ((a = (int*) calloc(NUM_ALPHA + 1, sizeof(int))) == NULL || \
        (b = (int*) calloc(n, sizeof(int))) == NULL || \
        (r = (char*) malloc(sizeof(char) * (n + 1))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);        
    }
    for (i = 0; i < n; i++) {  /* L列中每种字符的个数 */
        if (L[i] == MARKER) 
            a[0]++;
        else
            a[L[i] - 'a' + 1]++;
    }
    for (i = 1; i < NUM_ALPHA + 1; i++) {   /* F列中排在每种字符前面的其他字符的个数 */
        a[i] += a[i - 1];
    }
    for (i = 0; i < n; i++) {    /* L列中每个字符跳转到F列中的位置 */
        if (L[i] == MARKER)
            b[i] = 0;
        else
            b[i] = a[L[i] - 'a']++;
    }
    for (i = 0, pos = 0; i < n; i++) {
        r[n - 1 - i] = L[pos];
        pos = b[pos];
    }
    r[n] = '\0';
    free(a);
    free(b);
    return r;
}

int main(void) {
    char s[MAXSTR];
    char *L;
    int n;
    char *r;
    printf("input str: ");
    fgets(s, MAXSTR - 1, stdin);
    n = strlen(s);
    if (s[n - 1] == '\n')
        s[--n] = '\0';
    L = bwtEncode(s, n);
    printf("The L column: %s\n", L);
    r = bwtDecode(L, ++n);
    printf("The original str: %s\n", r);
    free(L);
    free(r);
}

实现代码二

如上面一小节所说的，根据 $L_i$ 找到同字符的 $F_j$ 还有另一种实现方法，就是每遇到一个 $L_i$ ，就计算一次 $L_i$ 前面有几个字符和它相同。理论上，这种查找方法需要 $O(n^2)$ 的时间，而上一小节的方法只需要线性时间。但是，这种方法也有好处，就是不需要构建“跳转数组”，从而节省空间。

效果如下：
在这里插入图片描述

具体代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSTR 1000
#define MARKER '$'
#define NUM_ALPHA 26

int comp(const void *s, const void *t) {
    return strcmp(*(char**)s, *(char**) t);    /* 注意这里 *(char**) 的用法 */
}

/* the last char of s is not MARKER */
char* bwtEncode(char *s, const int n) {
    char *L;
    char **M;
    int i, j, l;
    if ((M = (char**) malloc(sizeof(char*) * n)) == NULL || \
        (L = (char*) malloc(sizeof(char) * (n + 2))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < n; i++) {
        if ((M[i] = (char*) malloc(sizeof(char) * (i + 2))) == NULL) {   /* 只需保存开头到标记字符的那一部分字符串，这样总共可以节省大约一半的空间 */
            fputs("Error: out of space!\n", stderr);
            exit(1);       
        }
        for (j = 0; j < i + 1; j++)
            M[i][j] = s[n - 1 - i + j];   /* 这里的字符串没有存储 MARKER */
        M[i][i + 1] = '\0';
    }
    qsort(M, n, sizeof(M[0]), comp);   /* 对旋转后的多个字符串排序 */
    for (i = 0, L[0] = s[n - 1]; i < n; i++) {
        if ((l = strlen(M[i])) < n)
            L[i + 1] = s[n - 1 - l];
        else
            L[i + 1] = MARKER;
    }
    L[n + 1] = '\0';
    for (i = 0; i < n; i++)
        free(M[i]);
    free(M);
    return L;
}

char* bwtDecode(char *L, const int n) {
    int i, j, k;
    int *a;
    char *r;   /* original string. */
    int pos;
    if ((a = (int*) calloc(NUM_ALPHA + 1, sizeof(int))) == NULL || \
        (r = (char*) malloc(sizeof(char) * (n + 1))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);        
    }
    for (i = 0; i < n; i++) {  /* L列中每种字符的个数 */
        if (L[i] == MARKER) 
            a[0]++;
        else
            a[L[i] - 'a' + 1]++;
    }
    for (i = 1; i < NUM_ALPHA + 1; i++) {   /* F列中排在每种字符前面的其他字符的个数 */
        a[i] += a[i - 1];
    }
    for (i = 0, pos = 0; i < n; i++) {
        r[n - 1 - i] = L[pos];
        for (j = 0, k = 0; j < pos; j++)
            if (L[j] == L[pos])
                k++;
        pos = a[L[pos] - 'a'] + k;
    }
    r[n] = '\0';
    free(a);
    return r;
}

int main(void) {
    char s[MAXSTR];
    char *L;
    int n;
    char *r;
    printf("input str: ");
    fgets(s, MAXSTR - 1, stdin);
    n = strlen(s);
    if (s[n - 1] == '\n')
        s[--n] = '\0';
    L = bwtEncode(s, n);
    printf("The L column: %s\n", L);
    r = bwtDecode(L, ++n);
    printf("The original str: %s\n", r);
    free(L);
    free(r);
}