阅读文章 · 拾文

关于 Burrows-Wheeler 变换（Burrows-Wheeler Transform，BWT），常见的说法有三种：“BWT 能压缩数据”；“FM-index 的搜索时间与文本长度无关”；“BWA 就是一个 FM-index”。第一句不对：BWT 只是把字符重新排列，长度一个字节都不少，Burrows 和 Wheeler 在原始报告的摘要里就写明 “The transformation does not itself compress the data”。第二句只对计数成立，要报告每次出现的位置还得额外付出采样间隔那么多步。第三句漏掉了大半：BWA-MEM 只用 FM-index 找种子，比对质量靠的是后面的 Smith-Waterman 延伸。

本文回答四个问题：一个排列为什么能只凭一列字符还原（LF 映射）；bzip2 1.0.8 实际是怎样围绕 BWT 搭流水线的；FM-index 怎样用 rank 查询在 $O(m)$ 步内计数、又怎样用采样后缀数组定位；BWA 0.7.18 的源码在哪些地方和论文里的”标准做法”不同。后缀数组本身的构造（倍增、SA-IS）和基于后缀数组的二分搜索由系列的后缀数组负责，这里把后缀数组当作已知输入。

文中所有实验数字都来自同目录的 reproduce/fm.c 和 reproduce/run.sh，指标全部与时钟无关（游程数、熵、LF 步数、压缩后字节数），环境见第三节。

一、定义：排序所有旋转，取最后一列

给定长度为 $n-1$ 的文本，末尾追加一个比所有字符都小、且只出现一次的哨兵 $，得到长度为 $n$ 的 $T$。把 $T$ 的 $n$ 个循环旋转按字典序排序，排成一个 $n \times n$ 的矩阵 $M$；第一列记为 $F$，最后一列记为 $L$。$L$ 就是 BWT 的输出。

以 banana$ 为例的 BWT 矩阵：7 个循环旋转排序后，第一列 F 为 aaabnn，最后一列 L 为 annbaa；每行左侧标出行号和后缀数组值 SA，哨兵之后绕回的字符用灰色表示，第 4 行就是原文本本身

图里有三处值得看：

旋转顺序就是后缀顺序。 因为 $ 唯一且最小，两个旋转在比较到 $ 之前一定已经分出大小，$ 之后绕回来的灰色部分永远不参与比较。所以第 $i$ 行以后缀 $T[\mathrm{SA}[i]..]$ 开头，排序旋转等价于构造后缀数组 $\mathrm{SA}$。Burrows 和 Wheeler 在报告第 4.1 节正是用”追加 EOF 后排序后缀”来实现旋转排序的。
$L$ 可以直接从后缀数组读出：

\[ L[i] = T\big[(\mathrm{SA}[i] - 1) \bmod n\big]. \]

第 $i$ 行的最后一个字符，就是该后缀在原文本中的前一个字符。banana$ 的 $\mathrm{SA} = (6,5,3,1,0,4,2)$，于是 $L = \texttt{annb\$aa}$。 - $F$ 不需要存。 它只是 $T$ 的字符排序，由每个字符的出现次数就能还原。记 $C[c]$ 为 $T$ 中严格小于 $c$ 的字符个数，则 $F$ 中 $c$ 占据的行是 $[C[c],\ C[c+1])$。本例中 $C[\$]=0,\ C[a]=1,\ C[b]=4,\ C[n]=5$。

这个变换由 David Wheeler 在 1983 年于 AT&T 贝尔实验室发现，但当时没有发表；公开文本是 Burrows 与 Wheeler 在 1994 年 5 月 10 日发布的 DEC SRC Research Report 124（报告第 1 节自述）。报告里的原始定义没有哨兵：对原串的 $N$ 个旋转排序，输出 $L$ 以及原串所在的行号 $I$。bzip2 沿用的就是这种不带哨兵的形式，行号叫 origPtr（第三节）。

二、LF 映射与逆变换

只拿到 $L$ 为什么能还原 $T$？关键是 $L$ 和 $F$ 之间的一条对应关系。

LF 映射

引理（LF 映射）。 对任意字符 $c$，$L$ 中第 $k$ 个 $c$（按行号从上往下数，从 0 开始）与 $F$ 中第 $k$ 个 $c$ 是 $T$ 中同一个字符。于是第 $i$ 行的最后一个字符 $c = L[i]$ 在 $F$ 中位于第

\[ \mathrm{LF}(i) = C[c] + \mathrm{occ}(c, i) \]

行，其中 $\mathrm{occ}(c, i)$ 是 $L[0..i-1]$ 中 $c$ 的个数，即 $L$ 上的 rank 查询。

证明。 设 $i < j$ 且 $L[i] = L[j] = c$，第 $i$、$j$ 行的旋转分别是 $R_i < R_j$。把它们各自向右转一位，得到 $cR_i[0..n-2]$ 和 $cR_j[0..n-2]$，这两个旋转都以 $c$ 开头，而且恰好是 $L$ 中这两个 $c$ 在 $F$ 中所在的行。它们首字符相同，大小由 $R_i[0..n-2]$ 与 $R_j[0..n-2]$ 决定。两个不同的旋转不可能在前 $n-1$ 个字符上全等（否则剩下的一个字符也相等，两者就是同一个旋转），因此 $R_i < R_j$ 已经在前 $n-1$ 个字符上决出，右转后仍保持 $cR_i[0..n-2] < cR_j[0..n-2]$。所以以 $c$ 结尾的行，右转后按原来的相对顺序落在 $F$ 中 $c$ 的那一段里，第 $k$ 个对应第 $k$ 个。$\blacksquare$

LF 映射示意：左列是 L，右列是 F，每个字符带有它在本列同字母中的序号，连线把 L 中第 k 个 c 连到 F 中第 k 个 c；同一字母的连线互不交叉，底部给出 LF 公式和 C 数组的取值

图中同色的线互不交叉，这就是”相对顺序不变”。例如 $L[5] = a$ 是 $L$ 中第 1 个 a（第 0 个在第 0 行），$\mathrm{LF}(5) = C[a] + \mathrm{occ}(a, 5) = 1 + 1 = 2$，它是 $F$ 中第 1 个 a。

这个公式并不是 FM-index 的发明。Burrows 和 Wheeler 的报告第 4.2 节已经给出逆变换的两遍扫描：第一遍算出 $P[i]$（$L[i]$ 在 $L[0..i-1]$ 中出现的次数）和 $C$，第二遍用 $T[i] = P[i] + C[L[i]]$ 跳转。$P[i]$ 就是 $\mathrm{occ}(L[i], i)$。FM-index 的贡献在于：不再为每个位置存一个 $P[i]$，而是在压缩后的 $L$ 上支持任意 $(c, i)$ 的 rank 查询（第四、五节）。

逆变换

第 0 行是 $ 开头的那一行，它的 $L[0]$ 是 $T$ 中 $ 前面的字符，也就是文本的最后一个字符。每走一步 LF，就回到原文本中的前一个位置：

对应的代码（摘自 reproduce/fm.c，符号 0 表示 $，其余字节映射为 1 到 256）：

static int lf(const FM *f, int i)
{
    int c = f->L[i];
    return f->C[c] + occ(f, c, i);
}

/* Recover t[0..n-1] (t[n-1] = '$') by walking LF from row 0, the '$' row. */
static void fm_inverse(const FM *f, int *t)
{
    int i = 0;
    t[f->n - 1] = 0;
    for (int j = f->n - 2; j >= 0; j--) {
        t[j] = f->L[i];
        i = lf(f, i);
    }
}

逆变换共 $n-1$ 步，每步一次 $C$ 查表和一次 $\mathrm{occ}$。若像原始报告那样预先算好整张 $P$ 数组，每步 $O(1)$，总时间 $O(n)$，代价是每个位置多存一个整数。bzip2 在这里提供了两档：默认的解压内存约为 $100\text{k} + 4 \times$ 块大小；-s 改用 16 位加 4 位拼成的数组（decompress.c 中的 ll16 与 ll4），降到 $100\text{k} + 2.5 \times$ 块大小，手册说明代价是速度大约减半。

没有哨兵时（bzip2 的形式），若原串是某个子串的整数次重复，会出现完全相同的旋转。报告第 2 节指出逆变换依然正确，只是 LF 序列会重复访问同一批行。reproduce/fm.c 的测试专门生成了周期串，用原始报告的 $P$、$C$ 两遍扫描做无哨兵逆变换，2000 轮全部还原成功。

三、BWT 之后为什么好压缩：bzip2 1.0.8 的实际流水线

聚簇效应

报告第 3 节用英文里的 the 解释：所有以 he 开头的旋转在排序后挨在一起，它们的最后一个字符大多是 t。$L$ 中某一段因此集中出现少数几种字符，而 move-to-front（MTF，Bentley、Sleator、Tarjan 与 Wei，CACM 1986）编码正好擅长这种局部性：每个字符输出它在最近使用列表中的位置，再把它移到表头，连续重复的字符就变成一串 0。报告里的例子是：初始列表为 $(a, b, c, r)$，对 $L = \texttt{caraab}$ 编码得到 $(2, 1, 3, 1, 0, 3)$。

这种聚簇能实测。fm stats 对输入做与 bzip2 相同的无哨兵循环 BWT，统计 $L$ 的游程数 $r$，再对 $L$ 做 MTF（列表只含实际出现的字节，与 bzip2 一致），统计 0 的比例和 MTF 输出的零阶经验熵 $H_0$：

实验口径：两个文本文件取自 bzip2 1.0.8 源码包（run.sh 下载并校验 SHA-256）；随机 DNA 为均匀的 ACGT，种子 1；高度重复 DNA 是同一段 10,000 碱基随机序列的 100 份拷贝，每个碱基以 1‰ 概率替换为随机碱基，种子 1。最后一列是系统自带 bzip2 1.0.8 的 -9 输出字节数乘 8 再除以输入字节数。环境为 Intel Core i9-12900K、WSL2 内核 6.6.87.2、GCC 16.1.1，编译参数 -O2 -Wall -Wextra；程序不计时，连续运行 3 次输出逐字节一致。

读这张表时注意三点：

文本上，BWT 把游程数压到原来的约三分之一到四成，MTF 输出七成以上是 0，零阶熵从约 5 bit/字节降到约 2 bit/字节。bzip2 真正的熵编码（下文的零游程编码加多表 Huffman）比零阶熵估计还要再好一截，最终是 1.69 和 1.92 bit/字节。
随机数据上 BWT 毫无帮助。 随机 DNA 的游程数几乎不变，bzip2 输出 2.19 bit/碱基，比直接用 2 bit 打包还大。BWT 利用的是上下文可预测性，没有可预测性就没有收益。
重复数据上游程数随重复度骤降。 100 份近似拷贝让 $n/r$ 达到 79，这正是第八节 r-index 一类结构的出发点。

压缩率的理论解释来自 Manzini（JACM 2001）：他在不假设信源模型的前提下，证明了 BWT 加 MTF 的原始算法以及加游程编码的变体，其压缩率在最坏情况下都能用输入的 $k$ 阶经验熵 $H_k$ 界住，对任意 $k \ge 0$ 成立。此前的分析都假设输入来自有限阶马尔可夫信源。

bzip2 的流水线

bzip2 的手册自己说得很清楚：“bzip2 is not research work”，它是把现有想法工程化；手册列出的四份核心文献是 Burrows-Wheeler 报告、Hirschberg 与 LeLewer 的前缀码解码、Wheeler 的多表 Huffman 程序 bred3，以及 Bentley 与 Sedgewick 的字符串排序论文。按 1.0.8 源码，一个块的处理顺序如下：

flowchart TD
    A["input bytes"] --> B["RLE1 in bzlib.c: a run of 4..255 equal bytes becomes 4 bytes + one count byte"]
    B --> C["block: at most 100000 x k - 19 bytes after RLE1, k = 1..9"]
    C --> D["BZ2_blockSort: sort cyclic rotations, record origPtr"]
    D --> E["generateMTFValues: MTF over symbols in use, zero runs coded as RUNA / RUNB, then EOB"]
    E --> F["sendMTFValues: 2 to 6 Huffman tables, a selector for every 50 symbols"]
    F --> G["block header: magic, block CRC, randomised bit, origPtr in 24 bits"]

与常见的”BWT、MTF、RLE、Huffman 四步”相比，有几处出入：

游程编码有两次，而且第一次在 BWT 之前。 bzlib.c 的 ADD_CHAR_TO_BLOCK 和 add_pair_to_block() 把 4 到 255 个相同字节写成 4 个字节加一个计数字节（值为长度减 4）。BWT 之后的”游程编码”只针对 MTF 输出中的 0，并且合并在 generateMTFValues() 里：一段 0 的长度用 BZ_RUNA、BZ_RUNB 两个符号按双射二进制（RUNA 为 1、RUNB 为 2，低位在前）写出。这正是原始报告第 5 节建议的”用一个表示游程长度的码代替一串 0”。
MTF 列表只包含块中出现过的字节（unseqToSeq 映射），字母表大小是实际用到的符号数加 RUNA、RUNB 和块结束符 EOB。
Huffman 不是一张表。 sendMTFValues() 按 MTF 符号数选 2 到 6 张码表（少于 200 个符号用 2 张，至少 2400 个用 6 张），每 50 个符号（BZ_G_SIZE）用一个 selector 指定码表，码长上限 23（BZ_MAX_CODE_LEN）。
块大小以十万字节为单位。 -1 到 -9 对应 100,000 到 900,000 字节（默认 -9），bzlib.c 中的上限是 nblockMAX = 100000 * blockSize100k - 19，而且这是 RLE1 之后的字节数。手册给出的内存估算是压缩 $400\text{k} + 8 \times$ 块大小、解压 $100\text{k} + 4 \times$ 块大小（-s 时为 $2.5 \times$）。

块越大，BWT 能看到的上下文越多。原始报告的 Table 2 给出了 Calgary 语料 book1 在不同块大小下的压缩率：1 KB 时 4.34 bit/字符，64 KB 时 3.00，256 KB 时 2.68，750 KB（整个文件）时 2.49；在 1 亿字节级的 Hector 语料上，块增大到约 1 亿字节时降到 2.01。bzip2 手册在介绍 -1 到 -9 时也提醒，更大的块收益递减很快，“大部分压缩来自前两三百 KB 的块大小”。

旋转排序：bzip2 自己的实现

bzip2 没有使用通用的后缀数组构造算法。blocksort.c 的 BZ2_blockSort() 分两条路：

块小于 10,000 字节时直接用 fallbackSort()，源码注释称之为 “exponential radix sort”，受 Manber-Myers 后缀数组构造启发，也就是倍增思想。
否则用 mainSort()：先按前两个字节基数排序，再对每个桶用 mainQSort3()（注释写明实现的是 Bentley 与 Sedgewick 的三路字符串快排）加 mainSimpleSort()（希尔排序）。排序过程中消耗一个工作量预算 budget = nblock * ((wfact-1) / 3)，预算耗尽说明输入”too repetitive”，改走 fallbackSort()。

这套”快路径加倍增兜底”的设计是 0.9.5 引入的，此前 bzip2 靠随机化扰动块内容来避开坏情况；compress.c 注释写明，从 0.9.5 起随机化位总是写 0，解码端仍保留处理能力以兼容旧文件。排序完成后，origPtr 就是 ptr[i] == 0 的那个 $i$，以 24 位写入块头。

四、FM-index：用 rank 查询做 backward search

Ferragina 和 Manzini 在 FOCS 2000 的论文 “Opportunistic Data Structures with Applications” 中提出了后来被称为 FM-index 的结构。摘要给出的界是：文本 $T[1,u]$ 以每字符 $O(H_k(T)) + o(1)$ 比特存储（对任意固定 $k$），查找模式 $P[1,p]$ 的全部 $occ$ 次出现耗时 $O(p + occ\log^{\epsilon} u)$。期刊版（JACM 2005，“Indexing Compressed Text”）把第一种结构的空间写成 $5nH_k(T) + o(n)$ 比特、时间 $O(p + occ \log^{1+\epsilon} n)$。“opportunistic” 指文本越可压缩，索引越小，而查询不显著变慢。

它由三部分组成：$C$ 数组（$\sigma$ 个整数）、支持 $\mathrm{occ}(c, i)$ 的 $L$（第五节讨论怎么存）、以及用于定位的采样后缀数组（第六节）。原文本不必另存，第二节的逆变换就能把它还原出来。

区间不变量

模式 $P[0..m-1]$ 在 $T$ 中的出现，一一对应于 BWT 矩阵中以 $P$ 开头的行，这些行在排序矩阵里连续，构成一个区间 $[lo, hi)$。backward search 从 $P$ 的最后一个字符开始，每次在左边加一个字符，维护不变量：

处理完 $P[k..m-1]$ 后，$[lo, hi)$ 恰好是以 $P[k..m-1]$ 开头的行。

加上字符 $c = P[k-1]$ 时：以 $cP[k..m-1]$ 开头的行，左转一位后以 $P[k..m-1]$ 开头、以 $c$ 结尾；反过来，区间内 $L[i] = c$ 的行经 LF 右转一位后正好以 $cP[k..m-1]$ 开头。所以新区间就是 $\{\mathrm{LF}(i) : lo \le i < hi,\ L[i] = c\}$。由 LF 引理，同一字符的 LF 保序，这些像是 $F$ 中 $c$ 段里连续的一截，起点跳过了区间之前的 $\mathrm{occ}(c, lo)$ 个 $c$：

\[ lo' = C[c] + \mathrm{occ}(c, lo), \qquad hi' = C[c] + \mathrm{occ}(c, hi). \]

$lo' \ge hi'$ 时模式不出现。$P$ 不含 $，所以”旋转以 $P$ 开头”与”后缀以 $P$ 开头”等价，匹配不会绕过文本末尾。

图中每一栏橙色框出的 $L$ 字符，就是下一步会被 LF 带走的那些行。例如第二栏区间 $[1,4)$ 内有两个 n（第 1、2 行），加入 n 后它们被映射到 $F$ 中 n 段的第 5、6 行。最终区间 $[2,4)$ 含 2 行，ana 出现 2 次；整个过程只查了 $C$ 和 $\mathrm{occ}$，没有读原文本。

对应代码（摘自 reproduce/fm.c）：

/* Backward search: rows [*lo, *hi) are the suffixes prefixed by p[0..m-1]. */
static int fm_count(const FM *f, const int *p, int m, int *lo, int *hi)
{
    int s = 0, e = f->n;
    for (int k = m - 1; k >= 0 && s < e; k--) {
        int c = p[k];
        s = f->C[c] + occ(f, c, s);
        e = f->C[c] + occ(f, c, e);
    }
    *lo = s; *hi = e;
    return e > s ? e - s : 0;
}

对拍

./fm test 2000 生成 2000 段随机文本（字母表大小取 1、2、4、26、255，长度 1 到 3000，三分之一是周期串），每段查 60 个模式（一半取自文本的子串，一半随机生成、可能包含文本中没有的字符），逐一与朴素的逐位置比较计数对照，并用两种采样方式定位全部出现、排序后与朴素结果比较；同时检查带哨兵的逆变换、无哨兵的循环逆变换，以及长度不超过 400 时用朴素旋转排序核对后缀数组。结果：

rounds=2000 queries=120000 located=30160824 inverse=2000 cyclic_inverse=2000 naive_sa_checks=516 failures=0

同一程序用 -fsanitize=address,undefined 编译后跑 300 轮，没有报错。

复杂度要看清楚是什么的复杂度

计数要 $m$ 轮，每轮两次 $\mathrm{occ}$，总共 $2m$ 次 rank 查询。若 rank 是 $O(1)$，计数就是 $O(m)$，与 $n$ 无关。作为对照，后缀数组上的二分查找需要 $O(m\log n)$ 次字符比较，借助 LCP 信息可降到 $O(m + \log n)$（Manber 与 Myers，SIAM J. Comput. 1993；细节见后缀数组）。

“与 $n$ 无关”只说明操作次数，没有说明每次操作的代价。每一轮的 $lo$、$hi$ 可能落在 $L$ 的任意位置，文本大到远超缓存时，一次 rank 基本就是一次缓存未命中。BWA-MEM2（Vasimuddin 等，IPDPS 2019）在保持输出与 BWA-MEM 完全一致的前提下重写了几个核心内核，摘要列出的手段包括改善缓存复用、软件预取和 SIMD，端到端计算时间在单线程上最多快 3.5 倍（论文数据，未在本站复现）。FM-index 在实际系统中的主要成本是访存，而不是指令数。

五、occ 怎么存：checkpoint、位向量 rank 与小波树

最直接的办法是为每个位置、每个字符存一个计数，空间 $\sigma n$ 个整数，远大于文本本身。实际实现都在”存多少计数”和”查询时扫多少字符”之间取舍。

checkpoint 加扫描

每隔 $b$ 个位置存一份全部字符的累计计数，查询 $\mathrm{occ}(c, i)$ 时取 $\lfloor i/b \rfloor$ 处的计数，再扫描不超过 $b$ 个字符。reproduce/fm.c 取 $b = 64$，这是教学用的写法。

BWA 0.7.18 针对 DNA 把这一思路做到了缓存行粒度。bwt.h 定义 OCC_INTERVAL 为 128（OCC_INTV_SHIFT 7），宏 bwt_occ_intv() 表明存储布局是：每 128 个碱基，先放 4 个 bwtint_t（uint64_t）计数，紧跟这 128 个碱基的 2 bit 编码，合计 $32 + 32 = 64$ 字节，也就是每个碱基 4 bit。bwt_occ() 取出 checkpoint 后，用 __occ_aux() 以位运算（SWAR）数 64 位字中某个 2 bit 符号的个数，并非调用 __builtin_popcountll。$ 不存入 BWT，而是记下它所在的行号 primary，查询时对行号做一次修正。

位向量 rank 与小波树

字母表是二元时，$\mathrm{occ}$ 就是位向量上的 $\mathrm{rank}_1$。Jacobson（FOCS 1989）给出了用 $o(n)$ 比特额外目录实现常数时间 rank 的构造：大块存绝对计数，小块存相对计数，块内用查表或 popcount 收尾。

一般字母表可以用小波树（wavelet tree，Grossi、Gupta 与 Vitter，SODA 2003）化归为位向量：根节点把字母表一分为二，用一个位向量记录每个字符属于哪一半，再把两半的子序列分别递归下去。

L = annbaa 上的小波树：根节点按 {,a} 与 {b,n} 划分，位向量为 0111000；左孩子序列 aaa 按 {} 与 {a} 划分，位向量为 1011；右孩子序列 nnb 按 {b} 与 {n} 划分，位向量为 110。查询 occ(a, 5) 时，先在根节点前 5 位上求 rank0 得 2，再在左孩子前 2 位上求 rank1 得 1

查询 $\mathrm{occ}(a, 5)$：a 属于根的左半，根节点前 5 位 01110 中有 2 个 0，转到左孩子的前 2 位；a 属于左孩子的右半，前 2 位 10 中有 1 个 1，答案是 1。核对：$L[0..4] = \texttt{annb\$}$ 中确实只有一个 a。每层一次位向量 rank，树高 $\lceil\log_2\sigma\rceil$，所以 $\mathrm{occ}$ 是 $O(\log\sigma)$，backward search 变为 $O(m\log\sigma)$。位向量总长 $n\lceil\log_2\sigma\rceil$ 比特，Navarro 的综述 “Wavelet Trees for All”（JDA 2014）讨论了用压缩位向量把空间进一步降到接近 $nH_0$ 的做法，以及小波树在 FM-index 之外的用途。

DNA 的 $\sigma = 4$，BWA 这种”2 bit 打包加 checkpoint”的专用布局比小波树少一层间接访问；字母表大的场景（自然语言文本、蛋白质序列）才更需要小波树。

六、定位：采样后缀数组

backward search 给出的是行区间，要报告出现位置，还需要这些行的 $\mathrm{SA}$ 值。完整的后缀数组要 $n\lceil\log_2 n\rceil$ 比特，比压缩后的 BWT 大得多。解决办法是只存一部分 $\mathrm{SA}$ 值，其余的靠 LF 走过去：$\mathrm{LF}(i)$ 所在行的后缀恰好从 $\mathrm{SA}[i] - 1$ 开始，所以从第 $i$ 行走 $k$ 步 LF 到达一个存了值的行 $j$，就有

\[ \mathrm{SA}[i] = \mathrm{SA}[j] + k. \]

按文本位置采样

Ferragina 和 Manzini 的做法相当于只保留文本位置是 $d$ 的倍数的后缀，即 $\mathrm{SA}[i] \bmod d = 0$ 的行。每走一步 LF，文本位置减一，因此至多 $d - 1$ 步必然碰到一个采样点。代价是要用一个 $n$ 比特的位向量标记哪些行被采样，再用 rank 算出样本在数组中的下标。

这张图接着第四节的例子：ana 的区间是第 2、3 行。第 2 行走了 3 步，恰好达到 $d - 1$ 的上界。代码（摘自 reproduce/fm.c）：

static int ts_locate(const TextSample *s, const FM *f, int i, int *steps)
{
    int k = 0;
    while (!(s->bits[i / 64] >> (i % 64) & 1)) { i = lf(f, i); k++; }
    uint64_t below = s->bits[i / 64] & ((1ULL << (i % 64)) - 1);
    *steps = k;
    return s->val[s->blk[i / 64] + __builtin_popcountll(below)] + k;
}

BWA 的做法：按行号采样

BWA 0.7.18 没有这样做。bwtindex.c 的 bwa_idx_build() 调用 bwt_cal_sa(bwt, 32)，只为行号是 32 的倍数的行保存 $\mathrm{SA}$ 值；查询函数如下（源码原样摘录）：

/* BWA v0.7.18, bwt.c, bwt_sa() */
bwtint_t bwt_sa(const bwt_t *bwt, bwtint_t k)
{
    bwtint_t sa = 0, mask = bwt->sa_intv - 1;
    while (k & mask) {
        ++sa;
        k = bwt_invPsi(bwt, k);
    }
    /* without setting bwt->sa[0] = -1, the following line should be
       changed to (sa + bwt->sa[k/bwt->sa_intv]) % (bwt->seq_len + 1) */
    return sa + bwt->sa[k/bwt->sa_intv];
}

好处是不需要位向量：判断一行是否被采样只要看行号的低位，样本下标就是 $k / d$。坏处是失去了 $d - 1$ 步的保证，LF 走到哪些行事先无法控制。注释里的取模也不是多余的：行走可能穿过 $ 所在的行，从文本位置 0 绕回末尾，BWA 用把 sa[0] 置为 $-1$ 的技巧抵消这一绕回，reproduce/fm.c 的 rs_locate() 则直接对 $n$ 取模。

两种方案在 $d = 32$、样本数相同时的实测（对全部 $n$ 行逐一定位，统计 LF 步数）：

按文本位置采样的平均值正好是 $(d-1)/2 = 15.5$，最大值是 $d - 1 = 31$，这是由构造保证的。按行号采样的平均步数约为 30，接近两倍；最大值达到数百步。一个粗略的解释是：若把 LF 访问到的行看成随机行，每一步命中采样行的概率是 $1/d$，步数近似服从几何分布，均值为 $d - 1 = 31$，尾部概率约为 $(1 - 1/d)^k$，在 $10^6$ 行中出现几百步的行并不意外。这个近似与实测相符，但它不是证明。

所以这是一笔空间换步数的账。以 64 位样本、$d = 32$ 计，样本本身占每行 2 bit；按文本位置采样还要再加每行 1 bit 的位向量和它的 rank 目录，BWA 省掉了这部分，代价是平均多走一倍的 LF 步、且没有最坏情况上界。BWA-MEM 对出现次数过多的种子另有限制：bwa mem -c 默认 500，“skip seeds with more than INT occurrences”（fastmap.c），不会为高重复种子逐个定位。

七、基因组比对：BWA 与 Bowtie 实际用了什么

谱系

FM-index 进入基因组学并非始于某一个工具。Lam 等人的 BWT-SW（Bioinformatics 2008）用压缩索引做局部比对；2009 年，面向二代测序短读段的三个工具几乎同时发表：

Bowtie（Langmead、Trapnell、Pop 与 Salzberg，Genome Biology，2009 年 3 月）：摘要称对人类基因组每 CPU 小时比对超过 2500 万条读段，内存约 1.3 GB，并提出了允许错配的质量感知回溯（quality-aware backtracking）。这是 2009 年硬件上的论文数据。
BWA（Li 与 Durbin，Bioinformatics，2009 年 5 月在线发表）：在 backward search 上做允许错配和空位的回溯搜索，摘要称比 MAQ 快约 10 到 20 倍、准确度相近。
SOAP2（R. Li 等，Bioinformatics，2009 年 6 月在线发表），同样基于 BWT 索引。

之后的代表性演进是 Li 在 2012 年提出的 FMD-index（Bioinformatics 2012，fermi 组装器论文），它支持序列的双向扩展，并给出了查找全部超级最大精确匹配（super-maximal exact match，SMEM）的快速算法；BWA-MEM（Li，arXiv:1303.3997，2013，预印本，未经同行评审）用 SMEM 作种子。摘要称它适用于 70 bp 到几 Mb 的序列，并能自动在局部与端到端比对之间选择。

BWA 0.7.18 源码里的几个事实

索引的是正反两条链。 bwa_idx_build() 第一次调用 bns_fasta2bntseq(fp, prefix, 0) 时 for_only = 0，把反向互补序列接在正向序列后面，再对拼接结果建一个 BWT。这就是 FMD-index 能双向扩展的基础：在同一个索引里，一条链上的向右扩展等价于另一条链上的向左扩展。
构造算法按长度切换。 未指定 -a 时，拼接后长度 l_pac 超过 50,000,000 用 bwtsw（bwt_gen.c，文件头标注 2004 年 Wong Chi Kwong 的 BWT 构造代码），否则用 is（is.c，来自 Yuta Mori 2008 年的 sais-lite，即 SA-IS）。
SA 采样间隔固定为 32，采样方式是上一节的按行号；occ 的 checkpoint 间隔为 128，布局见第五节。按这两个常数，.bwt 约占每个索引碱基 4 bit，.sa 约占 2 bit（64 位样本每 32 行一个）；这是从源码常量推算的，未实测。
BWA-MEM 的流程：mem_collect_intv() 调用 bwt_smem1a() 收集 SMEM 区间，默认最小种子长度 19（bwamem.c 中 min_seed_len = 19），用 bwt_sa() 把区间里的行换成参考序列坐标，成链后调用 ksw_extend2() 做带仿射空位罚分的延伸。FM-index 负责的是种子，不是整条比对。

八、争论与开放问题

种子该用 FM-index 还是哈希

FM-index 的长处是能在很小的内存里枚举任意长度的精确匹配，SMEM 这种变长种子离不开它。但同一位作者后来的 minimap2（Li，Bioinformatics 2018）改用 minimizer 哈希索引做稀疏的定长种子，摘要称它在准确度相当时比主流短读段比对工具快 3 到 4 倍，比长读段比对工具快 30 倍以上。另一边，BWA-MEM2（2019）选择保留 FM-index 和 BWA-MEM 的输出，把力气花在缓存和 SIMD 上。两条路线各有论文数据支持，它们的对比依赖各自选的数据集和准确度口径，目前没有一组公认的基准能裁决”种子索引应该用哪种”。从机制上能说的是：两类种子都要求精确匹配，测序错误会把精确匹配切短，错误率越高，能达到 19 bp 这类最小种子长度的 SMEM 越少；两者的差别在于变长种子要靠 FM-index 逐字符扩展，而定长 minimizer 只需一次哈希查表。

bzip2 作者对自己格式的反思

bzip2 手册里有一节 Seward 对文件格式的反思，他认为几处设计在事后看来并不必要：第一步的 RLE1 “entirely irrelevant”，原本是为了防止排序遇到长串重复字符时退化，而原始报告的 Q6a、Q6b 步骤已经说明块排序本身可以处理这种情况；随机化机制也没有必要；块的压缩后长度没有记录在流里，导致解压实现很绕；CRC32 可以换成更快的 Adler-32。他的结论是 bzip2 的格式”在我充分理解其性能后果之前就冻结了”。今天的 bzip2 流水线因此不能当作”BWT 压缩器的最优设计”来读，它是一个格式冻结后的工程产物。原始报告第 6 节也承认，逐字节的 MTF 本身会带来损失，块再大也不会逼近理论最优。

BWT 游程数 $r$ 与 LZ77 的关系

第三节的表里，高度重复 DNA 的 $n/r$ 达到 79。以 $r$ 为空间单位的索引由此而来：Gagie、Navarro 与 Prezza 的 r-index（SODA 2018，期刊版 JACM 2020）在 $O(r)$ 个字的空间内支持定位，适合同一物种大量基因组组成的高重复集合。长期悬而未决的是 $r$ 与其他重复度量之间的关系：Kempa 与 Kociumaka 在 FOCS 2020 的 “Resolution of the Burrows-Wheeler Transform Conjecture” 摘要中写道，此前对 $r$ 没有已知的非平凡上界，而几乎所有其他压缩方法都已证明与 LZ77 解析的短语数 $z$ 相差至多多对数因子；他们证明了对任意文本 $r = O(z\log^2 n)$。这把 BWT 类索引的空间与 LZ77 可压缩性联系了起来，但对具体数据集而言，$r$ 和 $z$ 的实际比值仍要靠测量。

九、工程上容易出错的地方

十、参考资料

规范与文档

bzip2 1.0.8 手册（源码包中的 manual.xml、bzip2.txt）：块大小与内存估算、-s 选项、文件格式反思一节、延伸阅读列表。

源码

bzip2 1.0.8（https://sourceware.org/pub/bzip2/bzip2-1.0.8.tar.gz，SHA-256 ab5a03176ee106d3f0fa90e381da478ddae405918153cca248e682cd0c4a2269）：bzlib.c 的 ADD_CHAR_TO_BLOCK、add_pair_to_block()、nblockMAX；compress.c 的 BZ2_compressBlock()、generateMTFValues()、sendMTFValues()；blocksort.c 的 BZ2_blockSort()、mainSort()、mainQSort3()、fallbackSort()；bzlib_private.h 的 BZ_RUNA、BZ_RUNB、BZ_N_GROUPS、BZ_G_SIZE；decompress.c 的 smallDecompress 路径。
BWA v0.7.18（https://github.com/lh3/bwa，tag v0.7.18）：bwt.h 的 OCC_INTERVAL、bwt_occ_intv()、bwtint_t；bwt.c 的 bwt_occ()、__occ_aux()、bwt_sa()、bwt_cal_sa()；bwtindex.c 的 bwa_idx_build()；bntseq.c 的 bns_fasta2bntseq()；bwamem.c 的 mem_collect_intv()、min_seed_len；fastmap.c 的 -c 选项。

核心论文

M. Burrows, D. J. Wheeler, “A Block-sorting Lossless Data Compression Algorithm”, DEC Systems Research Center, Research Report 124, 1994-05-10.
P. Ferragina, G. Manzini, “Opportunistic Data Structures with Applications”, FOCS 2000, 390–398.
P. Ferragina, G. Manzini, “Indexing Compressed Text”, Journal of the ACM 52(4), 2005, 552–581.
G. Manzini, “An Analysis of the Burrows-Wheeler Transform”, Journal of the ACM 48(3), 2001, 407–430.
H. Li, R. Durbin, “Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform”, Bioinformatics 25(14), 2009, 1754–1760.
B. Langmead, C. Trapnell, M. Pop, S. L. Salzberg, “Ultrafast and Memory-efficient Alignment of Short DNA Sequences to the Human Genome”, Genome Biology 10(3), 2009, R25.

其他论文

J. L. Bentley, D. D. Sleator, R. E. Tarjan, V. K. Wei, “A Locally Adaptive Data Compression Scheme”, Communications of the ACM 29(4), 1986, 320–330.
J. L. Bentley, R. Sedgewick, “Fast Algorithms for Sorting and Searching Strings”, SODA 1997.
U. Manber, G. Myers, “Suffix Arrays: A New Method for On-line String Searches”, SIAM Journal on Computing 22(5), 1993, 935–948.
G. Jacobson, “Space-efficient Static Trees and Graphs”, FOCS 1989, 549–554.
R. Grossi, A. Gupta, J. S. Vitter, “High-order Entropy-compressed Text Indexes”, SODA 2003, 841–850.
G. Navarro, “Wavelet Trees for All”, Journal of Discrete Algorithms 25, 2014, 2–20.
T. W. Lam, W. K. Sung, S. L. Tam, C. K. Wong, S. M. Yiu, “Compressed Indexing and Local Alignment of DNA”, Bioinformatics 24(6), 2008, 791–797.
R. Li, C. Yu, Y. Li, T.-W. Lam, S.-M. Yiu, K. Kristiansen, et al., “SOAP2: an Improved Ultrafast Tool for Short Read Alignment”, Bioinformatics 25(15), 2009, 1966–1967.
H. Li, “Exploring Single-sample SNP and INDEL Calling with Whole-genome de novo Assembly”, Bioinformatics 28(14), 2012, 1838–1844.
H. Li, “Aligning Sequence Reads, Clone Sequences and Assembly Contigs with BWA-MEM”, arXiv:1303.3997, 2013（预印本）。
H. Li, “Minimap2: Pairwise Alignment for Nucleotide Sequences”, Bioinformatics 34(18), 2018, 3094–3100.
Md. Vasimuddin, S. Misra, H. Li, S. Aluru, “Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems”, IPDPS 2019, 314–324.
T. Gagie, G. Navarro, N. Prezza, “Optimal-Time Text Indexing in BWT-runs Bounded Space”, SODA 2018, 1459–1477；期刊版 “Fully Functional Suffix Trees and Optimal Text Searching in BWT-Runs Bounded Space”, Journal of the ACM 67(1), 2020.
D. Kempa, T. Kociumaka, “Resolution of the Burrows-Wheeler Transform Conjecture”, FOCS 2020, 1002–1013.

实验

reproduce/fm.c：BWT、逆变换（带哨兵与无哨兵）、checkpoint occ、backward search、两种 SA 采样，以及对拍、统计与数据生成；编译 gcc -O2 -Wall -Wextra -o fm fm.c -lm。
reproduce/run.sh：下载并校验 bzip2 源码包、生成两份 DNA 数据、运行全部实验；参考输出在 reproduce/results/stats.txt。
reproduce/draw_figures.py：从 banana$ 计算并生成本文全部分步 SVG 图。

系列导航： - 上一篇：AC 自动机：失败链接、输出链接与转移表布局 - 下一篇：编辑距离与模糊匹配：Wagner-Fischer、位并行与 Levenshtein 自动机

读完这篇，下一步读什么

优先读同系列或同问题的下一篇，把单篇消费变成主题集群。

2026-05-21 · algorithms

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