BWT 与 FM-index：从 bzip2 到基因组比对

1994 年，Michael Burrows 和 David Wheeler 在 DEC 的一份技术报告里提出了一种令人 拍案叫绝的字符串变换方法：把原始字符串的所有循环旋转按字典序排列，取最后一列， 就得到了一个对压缩极其友好的新序列。这个变换后来以他们的名字命名为 Burrows-Wheeler Transform（BWT）。

十年后，Paolo Ferragina 和 Giovanni Manzini 在 BWT 的基础上构建出了 FM-index， 一种支持在压缩空间内完成精确模式匹配的全文索引结构。这项工作不仅拿到了 ESA 2000 的最佳论文奖，更在 2009 年被 Heng Li 引入基因组学领域——BWA 工具的 诞生彻底改变了二代测序数据的比对范式。

本文将从 BWT 的定义出发，层层递进地讲解其可逆性证明、在 bzip2 中的压缩流水线、 FM-index 的查询机制、小波树对 rank 查询的加速、SA 采样的定位策略，以及在基因组 比对中的实际应用。文中还包含一份完整的 C 语言实现，约 200 行，可直接编译运行。

一、BWT 的定义

给定一个长度为 n 的字符串 T，我们在末尾追加一个字典序最小的终止符 $ （不出现在原始字母表中），得到 T$ = T[0..n]。BWT 的构造过程如下：

1. 生成 T$ 的所有 n+1 个循环旋转（cyclic rotation）。
2. 将这些旋转按字典序排序，形成一个 (n+1) x (n+1) 的矩阵 M，称为 BWT 矩阵（也叫 Burrows-Wheeler 矩阵）。
3. 取矩阵 M 的最后一列 L，即为 BWT(T)。

以字符串 banana 为例：

T$ = banana$

所有循环旋转：           排序后（BWT 矩阵 M）：
  banana$                  $banana          -> L[0] = a
  anana$b                  a$banan          -> L[1] = n
  nana$ba                  ana$ban          -> L[2] = n
  ana$ban                  anana$b          -> L[3] = b
  na$bana                  banana$          -> L[4] = $
  a$banan                  na$bana          -> L[5] = a
  $banana                  nana$ba          -> L[6] = a

BWT("banana") = "annb$aa"

下图展示了完整的变换流程：

我们把 BWT 矩阵的第一列记为 F 列（First column），最后一列记为 L 列（Last column）。F 列就是 T$ 中所有字符排序后的结果，可以直接从 字符频次得到；L 列就是 BWT 的输出。

形式化定义。 令 SA[0..n] 为 T$ 的后缀数组，则：

BWT[i] = T$[SA[i] - 1]    若 SA[i] > 0
BWT[i] = T$[n]  = $        若 SA[i] = 0

换言之，BWT 的每个位置记录的是后缀数组中对应后缀的前一个字符。这个等价定义 揭示了 BWT 与后缀数组之间的深层联系：给定后缀数组，BWT 可以在 O(n) 时间内 计算出来；反过来，给定 BWT，也可以在 O(n) 时间内恢复后缀数组（通过 LF-mapping 反复迭代）。

二、BWT 的可逆性证明

BWT 最令人惊叹的性质是：仅凭最后一列 L，就可以完整地恢复原始字符串 T。 这在直觉上很不可思议——我们丢弃了矩阵中绝大部分信息，只保留了一列， 却没有任何信息损失。

证明依赖于一个关键引理：LF-mapping 性质。

引理（LF-mapping）。 对于 BWT 矩阵中的任意字符 c，c 在 L 列中第 k 次出现 （从 0 计数）所在的行，与 c 在 F 列中第 k 次出现所在的行，对应的是同一个字符 c 在 T$ 中的同一次出现。

也就是说，相同字符在 F 列和 L 列中保持相同的相对顺序。

证明。 考虑 L 列中出现的某个字符 c。L[i] = c 意味着 BWT 矩阵第 i 行的 循环旋转以 c 结尾。将 c 移到该旋转的开头，得到的就是一个以 c 开头的旋转。 由于 BWT 矩阵是按字典序排列的，所有以 c 开头的旋转在矩阵中是连续的 （它们构成了 F 列中 c 的那一段）。

关键观察：对于 L 列中的两个相同字符 L[i] = L[j] = c（i < j），将它们各自移到 开头后，所得到的两个旋转的第二个字符开始的后续部分分别等于原来第 i 行和第 j 行 的完整旋转。由于第 i 行在第 j 行之前（字典序更小），移位后以 c 开头的两个旋转 也保持了相同的相对顺序。因此 c 在 L 列中的第 k 次出现映射到 c 在 F 列中的第 k 次出现。证毕。

有了 LF-mapping，我们可以定义映射函数：

LF(i) = C[c] + Occ(c, i)

其中：

• c = L[i]，即 L 列第 i 个字符。
• C[c] 是 F 列中严格小于 c 的字符总数（即字符 c 在 F 列中首次出现的位置）。
• Occ(c, i) 是 L[0..i-1] 中字符 c 的出现次数（即 c 在 L 列前 i 个位置中的 rank）。

逆变换算法。 从 L 列恢复 T 的过程如下：

算法 InverseBWT(L):
    计算 C[] 和 Occ 表
    i = 0                  // F 列中 $ 所在行（即第 0 行）
    for j = n-1 downto 0:
        T[j] = L[i]
        i = LF(i)
    return T

从 F 列中 $ 所在的行（第 0 行）开始，每次读取 L[i] 得到原始字符串的一个字符（倒序）， 然后通过 LF-mapping 跳到下一行。这个过程恰好遍历了整个字符串，因为每行 被恰好访问一次。

时间复杂度。 如果预先构建 Occ 表，逆变换只需 O(n) 时间。

三、BWT 与数据压缩：bzip2 流水线

BWT 本身并不压缩数据，但它把相似的上下文聚集在一起，使得输出中相同字符 倾向于连续出现。这种”聚簇效应”是 BWT 在压缩领域的核心价值。

bzip2 的压缩流水线由四个阶段组成：

原始数据
  |
  v
[1] BWT（Burrows-Wheeler 变换）
  |     将相似上下文聚集，产生大量连续重复字符
  v
[2] MTF（Move-To-Front 编码）
  |     将连续重复转化为大量小整数（0 和 1 居多）
  v
[3] RLE（游程编码）
  |     压缩连续的零
  v
[4] 熵编码（Huffman 编码）
        最终压缩

BWT 的聚簇效应

为什么 BWT 的输出对压缩友好？考虑英文文本中的 th。所有后缀以 h 开头 且前一个字符为 t 的位置，在 BWT 矩阵中对应的行是连续的（因为它们的后缀 排在一起），而这些行的 L 列值都是 t。类似地，he 前面经常出现 t， in 前面经常出现某些固定字符。这种局部聚簇使得同一字符在 L 列中大量 连续出现。

Move-To-Front 编码

MTF 维护一个初始为 [0, 1, 2, …, 255] 的字符列表。对于输入的每个字符， 输出该字符在列表中的位置，然后把它移到列表头部。如果输入中相同字符连续出现， MTF 的输出就是大量的 0。

BWT 输出：  a a a a b a a a c c a a
MTF 输出：  97 0 0 0 98 1 0 0 99 0 2 0

游程编码与 Huffman

RLE 对连续的零进行特殊编码，大幅减小数据量。最后的 Huffman 编码根据符号 频率分配变长码字，完成最终压缩。

bzip2 的块大小默认为 900KB，这意味着 BWT 矩阵的构建在每个块内独立进行。 块大小越大，BWT 的聚簇效果越好，但内存消耗也更大（BWT 矩阵的排序需要 约 8 倍块大小的内存）。

解压流程

解压是压缩的精确逆过程：

压缩数据 -> Huffman 解码 -> RLE 解码 -> MTF 逆变换 -> BWT 逆变换 -> 原始数据

BWT 逆变换利用上节介绍的 LF-mapping 在 O(n) 时间内完成，不需要重建 整个 BWT 矩阵。

四、FM-index：全文索引的革命

2000 年，Ferragina 和 Manzini 发表了 FM-index（Full-text index in Minute space），在 BWT 的基础上构建了一种支持精确模式匹配的压缩索引。

FM-index 的核心思想：BWT 矩阵的 F 列和 L 列之间的 LF-mapping 关系， 不仅可以用于逆变换，还可以用于字符串搜索。

数据结构组成

FM-index 由以下组件构成：

1. BWT 数组 L[0..n]：文本 T$ 的 BWT 变换结果。
2. C 数组：C[c] = F 列中字典序严格小于 c 的字符总数。对于大小为 σ 的字母表，C 数组只需 O(σ) 空间。
3. Occ 表（occurrence table）：Occ(c, i) = L[0..i-1] 中字符 c 的 出现次数。这是支持 rank 查询的核心。

C 数组的计算

// 统计字符频率
int freq[SIGMA] = {0};
for (int i = 0; i <= n; i++)
    freq[L[i]]++;

// 前缀和得到 C 数组
C[0] = 0;
for (int c = 1; c < SIGMA; c++)
    C[c] = C[c-1] + freq[c-1];

Occ 表的存储

朴素方案是存储一个 (σ x (n+1)) 的二维数组，空间 O(σn)。对于 DNA 序列 （σ = 5，包括 $ACGT），这是可以接受的。对于一般文本（σ = 256）， 需要更紧凑的方案。

常见的做法是分块存储：每隔 k 个位置存一个完整的 Occ 快照（checkpoint）， 查询时从最近的 checkpoint 开始扫描。k 的选择是空间与时间的权衡：

• k = 1：O(σn) 空间，O(1) 查询。
• k = 64：O(σn/64) 空间，O(64) 查询（但可以用 popcount 加速到 O(1)）。
• k = log²n：O(σn/log²n) 空间，O(log²n) 查询。

五、Backward Search：O(m) 精确匹配

FM-index 的精确匹配算法称为 backward search（逆向搜索）。给定模式串 P[0..m-1]，算法从右向左逐字符扫描 P，在每一步维护一个 SA 区间 [lo, hi)， 表示 BWT 矩阵中所有以当前已处理后缀为前缀的行的范围。

算法 BackwardSearch(P[0..m-1]):
    lo = 0
    hi = n + 1
    for i = m-1 downto 0:
        c = P[i]
        lo = C[c] + Occ(c, lo)
        hi = C[c] + Occ(c, hi)
        if lo >= hi:
            return "未找到"
    return "出现 hi - lo 次，SA 区间 [lo, hi)"

正确性说明。 初始时 [lo, hi) = [0, n+1) 表示整个 SA 范围。处理 P[m-1] 后，区间缩小为所有以 P[m-1] 开头的后缀。每处理一个新字符 P[i]，区间 进一步缩小为所有以 P[i..m-1] 为前缀的后缀。

每一步的更新公式基于 LF-mapping 的性质：

• C[c] 跳过所有首字符小于 c 的后缀。
• Occ(c, lo) 跳过在当前区间之前、L 列中为 c 的位置。
• Occ(c, hi) 确定区间的上界。

时间复杂度。 m 次迭代，每次一次 C 查询（O(1)）和两次 Occ 查询。 如果 Occ 查询为 O(1)（通过预计算或 popcount），则总时间为 O(m)。 这个复杂度与文本长度 n 无关，是 FM-index 最令人惊艳的特性。

对比后缀数组。 后缀数组 + LCP 数组的二分查找需要 O(m log n) 时间， FM-index 将其降到了 O(m)。

搜索示例

在 BWT(“banana”) = “annb$aa” 上搜索模式 “ana”：

初始：lo=0, hi=7

处理 P[2]='a': lo = C['a'] + Occ('a',0) = 1 + 0 = 1
               hi = C['a'] + Occ('a',7) = 1 + 3 = 4
               区间 [1,4)  -> 3 个以 'a' 开头的后缀

处理 P[1]='n': lo = C['n'] + Occ('n',1) = 5 + 0 = 5
               hi = C['n'] + Occ('n',4) = 5 + 2 = 7
               区间 [5,7)  -> 2 个以 "na" 开头的后缀

处理 P[0]='a': lo = C['a'] + Occ('a',5) = 1 + 1 = 2
               hi = C['a'] + Occ('a',7) = 1 + 3 = 4
               区间 [2,4)  -> 2 个以 "ana" 开头的后缀

结果："ana" 在 "banana" 中出现 2 次。

六、小波树与 Rank 查询

当字母表较大时（如 ASCII 的 σ=256，或 Unicode），朴素 Occ 表的空间消耗 变得不可接受。小波树（wavelet tree） 是解决这个问题的优雅方案。

基本结构

小波树是一棵二叉树，每个节点对应字母表的一个子集。根节点对应整个字母表， 左子树对应字典序较小的一半字符，右子树对应较大的一半。每个节点存储一个 位向量（bit vector），其中第 i 位表示该节点序列中第 i 个字符属于左子集（0） 还是右子集（1）。

字母表 {a, b, c, d}
分割：{a,b} | {c,d}

        root: [a,b,c,d,a,c,b,a,d]
        bits: [0,0,1,1,0,1,0,0,1]
       /                          \
  {a,b}: [a,b,a,b,a]         {c,d}: [c,d,c,d]
  bits:  [0,1,0,1,0]         bits:  [0,1,0,1]
   /       \                   /       \
 {a}:[a,a,a] {b}:[b,b]    {c}:[c,c] {d}:[d,d]

Rank 查询

要计算 Occ(c, i)（即 rank_c(i)），从根节点开始向下遍历：

1. 在当前节点的位向量上做 rank 查询，更新位置 i。
2. 如果 c 属于左子集，向左走，i = rank_0(i)。
3. 如果 c 属于右子集，向右走，i = rank_1(i)。
4. 到达叶节点时，i 就是答案。

树高为 O(log σ)，每层做一次位向量 rank（O(1) 时间），所以总时间为 O(log σ)。

位向量的 Rank 实现

位向量上的 rank 查询是整个结构的基石。经典实现使用两层索引：

大块（superblock）：每 256 位一个，存储前缀 popcount（绝对值）
小块（block）：    每 64 位一个，存储相对于大块的 popcount（相对值）
查询时：rank_1(i) = superblock[i/256] + block[i/64] + popcount(word & mask)

利用硬件 popcount 指令（x86 的 POPCNT，GCC 的 __builtin_popcountll）， 每次 rank 查询只需常数时间。

空间分析

小波树总共存储 n * ceil(log σ) 位的位向量，加上 rank 的辅助结构 （每个位向量额外 o(n) 位）。总空间为 n * log σ + o(n log σ) 位， 接近信息论下界 n * H_0(T)。

对于 FM-index，使用小波树后，backward search 的每步时间从 O(1) 变为 O(log σ)，总时间为 O(m log σ)。在实践中，对于 DNA 序列 （σ = 4 或 5），log σ 只有 2-3，几乎可以忽略。

七、C 语言实现：BWT 构建与 FM-index 搜索

下面给出一份完整的 C 语言实现，包含 BWT 的朴素构建、FM-index 的查询， 以及逆变换。代码可以直接用 gcc -O2 -o bwt bwt.c 编译运行。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAX_N    100000
#define SIGMA    256

/* --------------- 后缀数组（朴素构建） --------------- */

static int sa[MAX_N + 1];
static int sa_n;
static const char *sa_text;

static int sa_cmp(const void *a, const void *b)
{
    int i = *(const int *)a;
    int j = *(const int *)b;
    return strcmp(sa_text + i, sa_text + j);
}

static void build_suffix_array(const char *t, int n, int *out_sa)
{
    sa_text = t;
    sa_n = n;
    for (int i = 0; i <= n; i++)
        out_sa[i] = i;
    qsort(out_sa, n + 1, sizeof(int), sa_cmp);
}

/* --------------- BWT 构建 --------------- */

static char bwt[MAX_N + 2];

static void build_bwt(const char *t, int n, int *suf_arr, char *out_bwt)
{
    for (int i = 0; i <= n; i++) {
        if (suf_arr[i] == 0)
            out_bwt[i] = t[n]; /* 即 '$' */
        else
            out_bwt[i] = t[suf_arr[i] - 1];
    }
    out_bwt[n + 1] = '\0';
}

/* --------------- FM-index 结构 --------------- */

typedef struct {
    int  n;                     /* 文本长度（含 $） */
    char *L;                    /* BWT 数组 */
    int  C[SIGMA];              /* C 数组 */
    int  (*occ)[SIGMA];         /* Occ[i][c] = L[0..i-1] 中 c 的出现次数 */
} FMIndex;

static void fm_build(FMIndex *fm, const char *bwt_str, int n)
{
    fm->n = n;

    /* 复制 BWT */
    fm->L = (char *)malloc(n + 1);
    memcpy(fm->L, bwt_str, n + 1);

    /* 分配并计算 Occ 表 */
    fm->occ = (int (*)[SIGMA])calloc((size_t)(n + 2), sizeof(int[SIGMA]));

    /* occ[0][c] = 0 对所有 c */
    for (int i = 0; i < n + 1; i++) {
        unsigned char c = (unsigned char)fm->L[i];
        for (int a = 0; a < SIGMA; a++)
            fm->occ[i + 1][a] = fm->occ[i][a];
        fm->occ[i + 1][c]++;
    }

    /* 计算 C 数组 */
    int freq[SIGMA] = {0};
    for (int i = 0; i <= n; i++)
        freq[(unsigned char)fm->L[i]]++;

    fm->C[0] = 0;
    for (int c = 1; c < SIGMA; c++)
        fm->C[c] = fm->C[c - 1] + freq[c - 1];
}

static void fm_free(FMIndex *fm)
{
    free(fm->L);
    free(fm->occ);
}

/* Backward search: 返回匹配次数，lo/hi 为 SA 区间 */
static int fm_search(const FMIndex *fm, const char *pattern,
                     int *out_lo, int *out_hi)
{
    int m = (int)strlen(pattern);
    int lo = 0;
    int hi = fm->n + 1;

    for (int i = m - 1; i >= 0; i--) {
        unsigned char c = (unsigned char)pattern[i];
        lo = fm->C[c] + fm->occ[lo][c];
        hi = fm->C[c] + fm->occ[hi][c];
        if (lo >= hi)
            return 0;
    }

    if (out_lo) *out_lo = lo;
    if (out_hi) *out_hi = hi;
    return hi - lo;
}

/* --------------- BWT 逆变换 --------------- */

static void inverse_bwt(const FMIndex *fm, char *out_text)
{
    int n = fm->n;

    /* 从 F 列中 '$' 所在的行（即第 0 行）开始 */
    int i = 0;
    for (int j = n - 1; j >= 0; j--) {
        out_text[j] = fm->L[i];
        unsigned char c = (unsigned char)fm->L[i];
        i = fm->C[c] + fm->occ[i][c];
    }
    out_text[n] = '\0';
}

/* --------------- 定位：通过后缀数组回查 --------------- */

static void fm_locate(const FMIndex *fm, const int *suf_arr,
                      const char *pattern)
{
    int lo, hi;
    int count = fm_search(fm, pattern, &lo, &hi);

    printf("Pattern \"%s\": %d occurrence(s)", pattern, count);
    if (count > 0) {
        printf(" at position(s):");
        for (int i = lo; i < hi; i++)
            printf(" %d", suf_arr[i]);
    }
    printf("\n");
}

/* --------------- 主程序 --------------- */

int main(int argc, char *argv[])
{
    const char *text = (argc > 1) ? argv[1] : "banana";
    int n = (int)strlen(text);

    /* 在末尾添加 $ */
    char *t = (char *)malloc(n + 2);
    memcpy(t, text, n);
    t[n] = '$';
    t[n + 1] = '\0';

    printf("Text: %s\n", text);
    printf("T$:   %s\n", t);

    /* 构建后缀数组 */
    build_suffix_array(t, n, sa);

    printf("SA:  ");
    for (int i = 0; i <= n; i++)
        printf(" %d", sa[i]);
    printf("\n");

    /* 构建 BWT */
    build_bwt(t, n, sa, bwt);
    printf("BWT:  %s\n", bwt);

    /* 构建 FM-index */
    FMIndex fm;
    fm_build(&fm, bwt, n);

    /* 测试搜索 */
    const char *patterns[] = {"ana", "an", "na", "ban", "xyz", NULL};
    for (int i = 0; patterns[i]; i++)
        fm_locate(&fm, sa, patterns[i]);

    /* 逆变换验证 */
    char *recovered = (char *)malloc(n + 1);
    inverse_bwt(&fm, recovered);
    printf("Recovered: %s\n", recovered);
    printf("Match: %s\n", strcmp(text, recovered) == 0 ? "YES" : "NO");

    /* 清理 */
    fm_free(&fm);
    free(t);
    free(recovered);
    return 0;
}

编译与运行：

gcc -O2 -o bwt bwt.c
./bwt banana

预期输出：

Text: banana
T$:   banana$
SA:   6 5 3 1 0 4 2
BWT:  annb$aa
Pattern "ana": 2 occurrence(s) at position(s): 3 1
Pattern "an": 2 occurrence(s) at position(s): 3 1
Pattern "na": 2 occurrence(s) at position(s): 4 2
Pattern "ban": 1 occurrence(s) at position(s): 0
Pattern "xyz": 0 occurrence(s)
Recovered: banana
Match: YES

代码说明

后缀数组构建。 这里使用了 qsort + strcmp 的朴素 O(n^2 log n) 方法， 仅用于演示。生产环境应使用 SA-IS 算法（O(n) 时间，O(n) 空间）或 DC3/Skew 算法。

Occ 表。 预计算了完整的二维 Occ 表，每个位置存储所有字符的累计频次。 这在 σ 较小时（如 DNA 的 5 个字符）非常高效，但对于大字母表 需要换用小波树或分块存储。

BWT 构建。 在实际应用中，BWT 通常直接从后缀数组推导，而后缀数组 由 SA-IS 在线性时间内构建。libdivsufsort 是广泛使用的高效实现。

八、SA 采样与位置定位

FM-index 的 backward search 返回的是 SA 区间 [lo, hi)，告诉我们模式串 出现了多少次。但要知道每次出现的具体位置，还需要查询后缀数组的值 SA[i]。

如果完整存储后缀数组，空间为 O(n log n) 位，这与 FM-index 的压缩目标矛盾。 解决方案是 SA 采样（SA sampling）。

采样策略

每隔 d 个位置存储一个后缀数组值。对于未被采样的位置 i，通过反复应用 LF-mapping 向前走，直到到达一个被采样的位置 j：

SA[i] = SA[j] + steps

其中 steps 是从 i 到 j 经过的 LF-mapping 步数。

算法 LocateSA(i):
    steps = 0
    while i 不是采样位置:
        i = LF(i)
        steps++
    return sampled_SA[i] + steps

空间与时间权衡

在实际基因组比对中，d = 32 或 d = 64 是常见选择。对于人类基因组 （n ≈ 3 x 10^9），完整 SA 需要 12 GB（每个值 4 字节）， d = 32 的采样只需约 375 MB。

ISA 采样

另一种策略是采样逆后缀数组（ISA）：每隔 d 个文本位置存储其在 SA 中的 排名。ISA 采样的好处是它与 BWT 的某些操作（如从文本位置出发的查询）更兼容。 BWA 同时使用了 SA 采样和 ISA 采样。

九、基因组比对：BWA-MEM 算法概述

BWA（Burrows-Wheeler Aligner）是目前使用最广泛的短读段比对工具之一。 BWA-MEM 是其最新的算法变体，适用于 70bp 到数 Mbp 的读段。

基因组比对的挑战

二代测序（NGS）每次运行产生数亿条短读段（reads），每条长 100-300bp。 这些读段需要被比对到参考基因组（如人类基因组，约 30 亿碱基对）上。 核心挑战在于：

• 数据量巨大：一次测序产生 TB 级数据。
• 允许错配和间隔（mismatch/indel）：测序有错误率。
• 重复序列：基因组中存在大量重复区域。
• 速度要求：比对需要在合理时间内完成。

BWA-MEM 的工作流程

BWA-MEM 的比对分为三个主要阶段：

阶段一：种子生成（Seeding）。 使用 FM-index 的 backward search 在参考基因组的 BWT 上寻找 超级最大精确匹配（SMEM，Super-Maximal Exact Match）。 SMEM 是读段与参考基因组之间的最长精确匹配片段，且在两端都无法继续延伸。

算法从读段的每个位置出发，利用 backward search 逐步延伸，记录每个位置能 达到的最长匹配。当匹配长度低于阈值（默认 19bp）时跳过。

阶段二：种子延伸（Extension）。 对每个种子（SMEM），使用带 affine gap penalty 的 Smith-Waterman 算法 进行局部比对延伸。BWA-MEM 使用 SSE2 指令集加速的条带化（banded） Smith-Waterman 实现（ksw 库）。

阶段三：输出与后处理。 生成 SAM/BAM 格式的比对结果，包括： - 比对位置和方向 - CIGAR 字符串（描述比对关系） - 比对质量（MAPQ） - 配对信息（paired-end reads）

FM-index 在 BWA 中的具体使用

BWA 对参考基因组构建 FM-index 的过程：

参考基因组 -> SA-IS 构建后缀数组 -> 推导 BWT -> 构建 Occ 表 -> SA 采样

人类基因组的 FM-index 占用约 3.2 GB 内存（BWT + Occ 表 + SA 采样）， 远小于完整后缀数组的 12 GB。

在搜索过程中，BWA 使用的 Occ 表采用 128 位为一个块 的分块存储方案， 每个块存储 4 个 DNA 碱基的累计频次。利用 __builtin_popcountll 指令， 每次 Occ 查询只需约 10 个时钟周期。

为什么不用哈希或后缀数组？

FM-index 的优势在于空间紧凑且支持变长模式匹配（不限于固定 k-mer 长度）。 但在三代长读段时代，minimap2 基于 minimizer 的哈希索引因其更好的种子 生成策略而逐渐成为主流。

十、性能基准：FM-index 对比后缀数组

理论复杂度的差异在实际运行中表现如何？下面给出一组基准测试的典型结果。

索引构建

FM-index 的构建需要先构建后缀数组，因此构建时间略长。但磁盘占用显著更小。

查询性能

FM-index 在存在性和计数查询上有明显优势，但定位操作因 SA 采样的开销 而略逊色。在大多数应用中（如基因组比对），计数和存在性查询远多于定位， 因此 FM-index 的综合性能更优。

Cache 行为

FM-index 的 backward search 在访问 Occ 表时存在较差的缓存局部性—— 每步可能跳到 Occ 表的任意位置。相比之下，后缀数组的二分查找虽然也有 随机访问，但最终几步通常在缓存行内。

实际工程中，通过将 BWT 和 Occ 表交错存储（interleaved layout）， 可以改善缓存命中率。BWA 和 Bowtie2 都采用了这种优化。

十一、工程实践中的陷阱

在实际部署 BWT 和 FM-index 时，有许多教科书不会提到的工程细节。 下表汇总了常见的陷阱及其应对方案。

十二、真实世界的应用

bzip2

bzip2 是 BWT 在压缩领域最成功的应用。由 Julian Seward 于 1996 年发布， 至今仍是 Linux 发行版中的标准工具。

bzip2 的核心参数是块大小（-1 到 -9，对应 100KB 到 900KB）。 更大的块意味着更好的压缩比，但也需要更多内存：

压缩内存   ≈ 块大小 x 8
解压内存   ≈ 块大小 x 4

bzip2 的压缩比通常优于 gzip（DEFLATE 算法），但劣于 xz（LZMA 算法）。 在压缩/解压速度上，bzip2 也慢于 gzip。bzip2 的主要优势在于解压速度 相对合理（因为 BWT 逆变换是线性的），以及良好的压缩比/内存权衡。

值得一提的是，bzip2 的后缀数组排序曾是其最大的性能瓶颈。早期版本使用 qsort，后来切换到了更高效的算法（类似于 Bentley-Sedgewick 的多键快排）。

Bowtie / Bowtie2

Bowtie 是 2009 年发表在 Genome Biology 上的短读段比对工具， 是 FM-index 在基因组学中的第一个大规模应用。Bowtie 使用 FM-index 支持精确匹配和允许少量错配的近似匹配（通过回溯搜索）。

Bowtie2 是其后续版本，增加了对间隔（indel）的支持，使用类似 BWA 的 种子-延伸策略。Bowtie2 的 FM-index 实现针对 DNA 的 4 字符字母表 做了深度优化：Occ 表使用 2-bit 编码存储 DNA 碱基，每 128 个碱基 一个 checkpoint。

BWA / BWA-MEM

BWA 由 Heng Li 开发，是目前引用量最高的生物信息学工具之一 （截至 2025 年超过 60,000 次引用）。BWA-MEM 算法自 2013 年以来 成为短读段比对的事实标准。

BWA 的 FM-index 实现有几个值得注意的工程决策：

• 双向 BWT：同时存储正向和反向文本的 BWT，支持正向和反向的 backward search，便于实现 SMEM 的高效枚举。
• 2-bit 编码：DNA 碱基使用 2 位编码（A=00, C=01, G=10, T=11）， $ 和 N 特殊处理。
• 交错存储：BWT 和 Occ checkpoint 交错存储在同一数组中， 一次缓存行加载即可获取所需的所有信息。

其他应用

• SDSL（Succinct Data Structure Library）：C++ 库，提供了 BWT、FM-index、小波树、压缩后缀数组等数据结构的高质量实现。
• r-index：基于 BWT 的 run-length 编码构建的索引，空间与 BWT 的 游程数成正比，适合高度重复的文本（如泛基因组数据）。
• pbzip2：bzip2 的并行版本，利用块独立性在多核上并行压缩/解压。

十三、我的一些思考

BWT 是我在学习字符串算法过程中遇到的最优美的构造之一。它的美在于 一种深刻的”无中生有”——通过一次看似毫无目的的矩阵变换，同时获得了 压缩友好性和可搜索性，这两个看起来毫不相关的性质。

我第一次实现 BWT 逆变换时，花了整整一个下午才真正理解 LF-mapping 为什么成立。那个”相同字符在 F 列和 L 列中保持相对顺序”的证明， 文字读起来很自然，但要在脑海中真正”看到”它，需要反复地画图、 举例、用手指在矩阵上一行一行地追踪。我建议每个学习 BWT 的人 都亲手在纸上做一遍这个过程。

FM-index 的 backward search 同样令人赞叹。传统的子串搜索 （如 KMP、Boyer-Moore）沿着文本从左到右滑动窗口，而 backward search 完全颠覆了这个思路：它不看文本本身，只操作一个经过变换的数组和 一些辅助表，就完成了搜索。更惊人的是，搜索时间只与模式串长度有关， 与文本长度完全无关。

从工程角度看，FM-index 的设计哲学值得玩味。它本质上是一种 “用计算换空间”的策略：通过复杂的编码和查表操作，在压缩表示上 直接完成原本需要在完整数据上进行的操作。这个思路在今天的内存墙时代 尤其有价值——当数据大到放不进内存时，压缩索引不只是节省了存储， 更是通过减少 I/O 而提升了实际性能。

在基因组学领域，BWT 和 FM-index 的影响怎么强调都不为过。 BWA 和 Bowtie 的出现使得在一台普通工作站上处理整个人类基因组 的测序数据成为可能。没有这些工具，二代测序革命不可能以如此快的 速度推进。

但技术总在演进。三代长读段测序（PacBio、Nanopore）的兴起使得 比对工具的设计发生了转变。minimap2 基于 minimizer 哈希的种子策略 在长读段场景下表现更好，因为长读段的高错误率使得 FM-index 的 精确匹配种子变得不够实用。不过，对于短读段和高精度的变异检测， BWA-MEM 仍然是首选。

BWT 和 FM-index 给我最大的启发是：好的算法设计往往来自对问题的 重新表述（reformulation）。BWT 把”压缩”重新表述为”排列使得相似 上下文聚集”；FM-index 把”搜索”重新表述为”在变换域上做区间收缩”。 这种能力——找到一个新的视角，使得原来困难的问题变得自然—— 是算法研究中最珍贵的素质。

最后说一个有趣的事实：BWT 和后缀数组之间存在一种”对偶”关系。 后缀数组记录的是”每个位置的后缀排在第几”，BWT 记录的是 “排在第几的后缀的前一个字符是什么”。看似不同的信息，却可以 互相推导，这是字符串算法中最精巧的对称性之一。

参考文献

1. Burrows, M. and Wheeler, D. J. (1994). A Block-sorting Lossless Data Compression Algorithm. DEC Systems Research Center Technical Report 124.
2. Ferragina, P. and Manzini, G. (2000). Opportunistic Data Structures with Applications. FOCS 2000, pp. 390-398.
3. Ferragina, P. and Manzini, G. (2005). Indexing Compressed Text. Journal of the ACM, 52(4):552-581.
4. Li, H. and Durbin, R. (2009). Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics, 25(14):1754-1760.
5. Li, H. (2013). Aligning Sequence Reads, Clone Sequences and Assembly Contigs with BWA-MEM. arXiv:1303.3997.
6. Langmead, B. et al. (2009). Ultrafast and Memory-efficient Alignment of Short DNA Sequences to the Human Genome. Genome Biology, 10:R25.
7. Grossi, R., Gupta, A. and Vitter, J. S. (2003). High-order Entropy-compressed Text Indexes. SODA 2003, pp. 841-850.
8. Navarro, G. (2014). Wavelet Trees for All. Journal of Discrete Algorithms, 25:2-20.
9. Nong, G., Zhang, S. and Chan, W. H. (2009). Linear Suffix Array Construction by Almost Pure Induced-Sorting. DCC 2009, pp. 193-202.
10. Gagie, T., Navarro, G. and Prezza, N. (2020). Fully Functional Suffix Trees and Optimal Text Searching in BWT-Runs Bounded Space. Journal of the ACM, 67(1):1-54.

上一篇: AC 自动机 下一篇: 编辑距离与模糊匹配 相关阅读: - 后缀数组 - Huffman 编码与 DEFLATE