序列比对（27）BWT算法
本文介绍了BWT算法。
bwa是目前最流行的二代测序比对工具，其中就用到了BWT算法。BWT（Burrows-Wheeler Transform）算法是一种数据转换算法，它将一个字符串中的相似字符放在相邻的位置，以便于后续的压缩。
简要回顾
BWT算法可以分为编码和解码两部分。编码后，原始字符串中的相似字符会处在比较相邻的位置；解码就是将编码后的字符串重新恢复成原始字符串的过程。BWT的一个特点就是经过编码后的字符串可以完全恢复成原始字符串。
用图示表示就是：
Step2：循环解码。
循环解码的具体步骤如下： Step2-1：
Step2-2：
Step2-3：
Step2-4：
Step2-5：
Step2-6：
Step3：输出原始字符串。
关键步骤
现在我们来说明为什么按照上述解码过程就可以恢复原始字符串。关键就是要解答两个问题：
我们重新看BWT编码中的“循环转移”这一步。我们将某一行字符串的Latter String定义为其在“循环转移”这一步中的下一行字符串；而将某一行字符串的Former String定义为其在“循环转移”这一步中的上一行字符串。
我们可以看出，某一行字符串的最后一个字符是其Latter String的第一个字符；某一行字符串的最后一个字符和其Latter String的最后一个字符的关系是：在原始字符串中，上述两个字符紧挨在一起并且Latter String的最后一个字符排在前面。
红色方框为所有以b开头的字符串；绿色方块为所有以b结尾的字符串。黑色箭头指向的是对应的Latter String。
实现代码一
代码的实现有几点要说明：
效果如下：
具体的实现代码如下:
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#define MAXSTR 1000
5#define MARKER '$'
6#define NUM_ALPHA 26
7
8int comp(const void *s, const void *t) {
9 return strcmp(*(char**)s, *(char**) t); /* 注意这里 *(char**) 的用法 */
10}
11
12/* the last char of s is not MARKER */
13char* bwtEncode(char *s, const int n) {
14 char *L;
15 char **M;
16 int i, j, l;
17 if ((M = (char**) malloc(sizeof(char*) * n)) == NULL || \
18 (L = (char*) malloc(sizeof(char) * (n + 2))) == NULL) {
19 fputs("Error: out of space!\n", stderr);
20 exit(1);
21 }
22 for (i = 0; i < n; i++) {
23 if ((M[i] = (char*) malloc(sizeof(char) * (i + 2))) == NULL) { /* 只需保存开头到标记字符的那一部分字符串，这样总共可以节省大约一半的空间 */
24 fputs("Error: out of space!\n", stderr);
25 exit(1);
26 }
27 for (j = 0; j < i + 1; j++)
28 M[i][j] = s[n - 1 - i + j]; /* 这里的字符串没有存储 MARKER */
29 M[i][i + 1] = '\0';
30 }
31 qsort(M, n, sizeof(M[0]), comp); /* 对旋转后的多个字符串排序 */
32 for (i = 0, L[0] = s[n - 1]; i < n; i++) {
33 if ((l = strlen(M[i])) < n)
34 L[i + 1] = s[n - 1 - l];
35 else
36 L[i + 1] = MARKER;
37 }
38 L[n + 1] = '\0';
39 for (i = 0; i < n; i++)
40 free(M[i]);
41 free(M);
42 return L;
43}
44
45char* bwtDecode(char *L, const int n) {
46 int i;
47 int *a, *b;
48 char *r; /* original string. */
49 int pos;
50 if ((a = (int*) calloc(NUM_ALPHA + 1, sizeof(int))) == NULL || \
51 (b = (int*) calloc(n, sizeof(int))) == NULL || \
52 (r = (char*) malloc(sizeof(char) * (n + 1))) == NULL) {
53 fputs("Error: out of space!\n", stderr);
54 exit(1);
55 }
56 for (i = 0; i < n; i++) { /* L列中每种字符的个数 */
57 if (L[i] == MARKER)
58 a[0]++;
59 else
60 a[L[i] - 'a' + 1]++;
61 }
62 for (i = 1; i < NUM_ALPHA + 1; i++) { /* F列中排在每种字符前面的其他字符的个数 */
63 a[i] += a[i - 1];
64 }
65 for (i = 0; i < n; i++) { /* L列中每个字符跳转到F列中的位置 */
66 if (L[i] == MARKER)
67 b[i] = 0;
68 else
69 b[i] = a[L[i] - 'a']++;
70 }
71 for (i = 0, pos = 0; i < n; i++) {
72 r[n - 1 - i] = L[pos];
73 pos = b[pos];
74 }
75 r[n] = '\0';
76 free(a);
77 free(b);
78 return r;
79}
80
81int main(void) {
82 char s[MAXSTR];
83 char *L;
84 int n;
85 char *r;
86 printf("input str: ");
87 fgets(s, MAXSTR - 1, stdin);
88 n = strlen(s);
89 if (s[n - 1] == '\n')
90 s[--n] = '\0';
91 L = bwtEncode(s, n);
92 printf("The L column: %s\n", L);
93 r = bwtDecode(L, ++n);
94 printf("The original str: %s\n", r);
95 free(L);
96 free(r);
97}实现代码二
效果如下：
具体代码如下：
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#define MAXSTR 1000
5#define MARKER '$'
6#define NUM_ALPHA 26
7
8int comp(const void *s, const void *t) {
9 return strcmp(*(char**)s, *(char**) t); /* 注意这里 *(char**) 的用法 */
10}
11
12/* the last char of s is not MARKER */
13char* bwtEncode(char *s, const int n) {
14 char *L;
15 char **M;
16 int i, j, l;
17 if ((M = (char**) malloc(sizeof(char*) * n)) == NULL || \
18 (L = (char*) malloc(sizeof(char) * (n + 2))) == NULL) {
19 fputs("Error: out of space!\n", stderr);
20 exit(1);
21 }
22 for (i = 0; i < n; i++) {
23 if ((M[i] = (char*) malloc(sizeof(char) * (i + 2))) == NULL) { /* 只需保存开头到标记字符的那一部分字符串，这样总共可以节省大约一半的空间 */
24 fputs("Error: out of space!\n", stderr);
25 exit(1);
26 }
27 for (j = 0; j < i + 1; j++)
28 M[i][j] = s[n - 1 - i + j]; /* 这里的字符串没有存储 MARKER */
29 M[i][i + 1] = '\0';
30 }
31 qsort(M, n, sizeof(M[0]), comp); /* 对旋转后的多个字符串排序 */
32 for (i = 0, L[0] = s[n - 1]; i < n; i++) {
33 if ((l = strlen(M[i])) < n)
34 L[i + 1] = s[n - 1 - l];
35 else
36 L[i + 1] = MARKER;
37 }
38 L[n + 1] = '\0';
39 for (i = 0; i < n; i++)
40 free(M[i]);
41 free(M);
42 return L;
43}
44
45char* bwtDecode(char *L, const int n) {
46 int i, j, k;
47 int *a;
48 char *r; /* original string. */
49 int pos;
50 if ((a = (int*) calloc(NUM_ALPHA + 1, sizeof(int))) == NULL || \
51 (r = (char*) malloc(sizeof(char) * (n + 1))) == NULL) {
52 fputs("Error: out of space!\n", stderr);
53 exit(1);
54 }
55 for (i = 0; i < n; i++) { /* L列中每种字符的个数 */
56 if (L[i] == MARKER)
57 a[0]++;
58 else
59 a[L[i] - 'a' + 1]++;
60 }
61 for (i = 1; i < NUM_ALPHA + 1; i++) { /* F列中排在每种字符前面的其他字符的个数 */
62 a[i] += a[i - 1];
63 }
64 for (i = 0, pos = 0; i < n; i++) {
65 r[n - 1 - i] = L[pos];
66 for (j = 0, k = 0; j < pos; j++)
67 if (L[j] == L[pos])
68 k++;
69 pos = a[L[pos] - 'a'] + k;
70 }
71 r[n] = '\0';
72 free(a);
73 return r;
74}
75
76int main(void) {
77 char s[MAXSTR];
78 char *L;
79 int n;
80 char *r;
81 printf("input str: ");
82 fgets(s, MAXSTR - 1, stdin);
83 n = strlen(s);
84 if (s[n - 1] == '\n')
85 s[--n] = '\0';
86 L = bwtEncode(s, n);
87 printf("The L column: %s\n", L);
88 r = bwtDecode(L, ++n);
89 printf("The original str: %s\n", r);
90 free(L);
91 free(r);
92}小结
本文比较详细地介绍了BWT算法，包括编码和解码过程，并且给出了解释。随后，给出了实现代码。遗留的一个问题是为什么经过BWT编码后，原始字符串中的相似字符会到相邻的位置？这个有待后续学习了。
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2019-10-24，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录