BWA-MEM 2013：序列比对算法

首页论文导航 BWA 2009 BWA-SW 2010

摘要

BWA-MEM 是一种新的比对算法，用于将序列读段或组装 contig 比对到大型参考基因组。它能自动在局部比对与端到端比对之间选择，支持双端测序和嵌合比对，对测序错误鲁棒，适用于 70bp 到数 Mb 的序列长度。100bp 读段比对优于当时多种最先进的比对工具。

获取： BWA 软件包组件。github.com/lh3/bwa

1. 引言

大多数短读段比对工具开发于读长约 36bp 的时代。随读段增长，两个新需求变得关键：

长插入缺失——需仿射空位罚分
多个不重叠的局部比对——结构变异或参考组装错误导致

已有长读比对算法各有缺陷：

算法	问题
BWA-SW	100bp 比 Bowtie2 慢，准确率相当；比 Cushaw2 精度低
Bowtie2 / Cushaw2	600bp 以上速度明显下降
GEM	强制端到端，无仿射空位比对

同时，从头组装产生的 contig 长度从几百 bp 到几 Mb，几乎没有算法能在这种跨度下保持高精度并正确处理易位。

2. 单条序列比对

2.1 FMD-index

BWA-MEM 使用 FMD-index（Li, 2012）——同时存储正向基因组 BWT 和反向互补基因组 BWT：

FMD-index(X) = (BWT(X), BWT(X_rc)) (1)

为什么需要双向？ SMEM 搜索需要在 read 上同时向左和向右扩展精确匹配。标准 BWT 只支持从 3'→5' 方向搜索。FMD-index 通过在反向互补链上建 BWT 等效支持了"向左"搜索。代价：内存翻倍，人类基因组约需 5.4 GB。

2.2 SMEM 搜索——种子生成

MEM（Maximal Exact Match）： 读段与参考间精确匹配片段，两端无法再延伸。SMEM（Super-Maximal Exact Match）： 不被任何其他 MEM 完全包含的 MEM。性质：每个位置最多被一个 SMEM 覆盖。

图 1 SMEM 搜索的三趟扫描流程。三趟依次执行：Forward Search 确定右端边界 → Backward Search 找左端扩展点 → LAST 遍历枚举 MEM 并过滤出 SMEM。

三趟扫描详解

Forward Search： 从 FMD-index 根节点开始，沿 read 从 5'→3' 方向逐个碱基做反向搜索。记录每个 position i 的 SA 区间。当区间变为空 → 该位置是某 SMEM 的右端边界。
Backward Search： 从每个右端边界出发，沿 read 从 3'→5' 做正向搜索。标记 SA 区间发生"收缩"的位置（LEP）。
LAST（轻量后缀数组遍历）： 对每个 LEP，在 SA 上快速找到最大扩展长度。输出所有 MEM → 过滤得到 SMEM。

2.3 Re-seeding

真正比对区域可能完全没有 SMEM。解决方案：

设 SMEM 长度 l、出现次数 k。
若 l > minSeedLen × 1.5（默认 28.5bp），在该 SMEM 中点处：
找覆盖中点、出现 ≥ k+1 次的最长精确匹配 → 新种子

参数： -r 1.5（再播种阈值），-c 500（丢弃出现超 500 次的 MEM）。

2.4 种子成链与过滤

共线且距离相近的种子构成链。贪心建链后过滤：短链被长链在 query 上覆盖超 50% 且短链比长链短至少 38bp → 丢弃。参数 -D 0.5。

2.5 种子延伸（Banded Affine-gap DP）

种子按（链长度 → 种子长度）排序。从最优到最差依次处理：若已被已有比对覆盖则跳过，否则 banded DP 延伸。

DP 参数：

参数	含义	默认值
-A	匹配得分	1
-B	错配罚分	4
-O	gap open 罚分	6
-E	gap extension 罚分	1
-L	clip penalty	5

Z-dropoff

bestScore − currentScore(x, y) > Z + |x − y| × p_gapExt (2)

Z 默认 100。加入 |x − y| × p_gapExt 项使坐标差大时阈值自动放宽——避免惩罚单边 long gap。

	BLAST X-dropoff	BWA-MEM Z-dropoff
阈值	固定 X	Z + \|x−y\| × p_gapExt
单边 gap	惩罚	不惩罚
生物学含义	假设一致	允许结构变异

Local vs End-to-End 自动选择

if bestScore_toEnd ≥ bestScore_local − clipPenalty → 端到端
else → 局部（soft-clip） (3)

意义：读尾有真实变异时端到端模式可穿过避免 reference bias；嵌合读段时局部模式只比对匹配部分。

3. 双端读段比对

3.1 批次处理

256K 对/批：加载 FMD-index → 单端比对 → 估计 μ, σ² → 配对 → 清除索引加载 2-bit 参考 → mate rescue。

3.2 Mate Rescue

估算 μ, σ² → 每端 top 100 hit → 窗口 [μ−4σ, μ+4σ] → SSE2 SW 搜索

记录第二好的 SW 分数——用于检测串联重复中的错误比对。参数：-S 跳过 rescue，-m 50 最大轮数。

3.3 配对评分

S_ij = S_i + S_j − min{ −a · log₄[P(d_ij)], U } (4)

符号	含义
S_i, S_j	两端 SW 分数
d_ij	插入片段长度
P(d)	正态分布 P(插入 ≥ d)
−a·log₄[P(d)]	插入偏离惩罚（匹配分空间）
U	未配对惩罚（默认 17）

为什么 log₄？ SW 分数解释为 log odds ratio 时，DNA 是 4 碱基系统，以 4 为底使插入惩罚与 SW 分在同一量纲。参数 -U 17。

条件	含义	结果
−a·log₄[P(d)] < U	插入片段合理	强制配对
−a·log₄[P(d)] ≥ U	插入异常大	扣 U 分，允许不配对

4. 结果

4.1 模拟与真实数据性能

条件	结果
SE 准确率	NovoAlign 最佳；BWA-MEM 接近；高于 GEM、Cushaw2
PE 准确率	BWA-MEM 接近 NovoAlign
速度 100bp SE	与 GEM、Bowtie2 相当
速度 650bp PE	比 Bowtie2 / Cushaw2 快约 6×

长读段优势原因：

SMEM 种子化： 一次遍历找到所有显著种子，复杂度与读段长度亚线性
Banded DP： O(L) 与 query 长度线性相关

4.2 E. coli 菌株比对

BWA-MEM vs nucmer，E. coli K-12 (4.6Mb) → 536 菌株：

指标	nucmer	BWA-MEM
时间	25s	131s
报告差异数	105,505	104,321
重叠	–	102,241

重叠 102,241 个。独有差异多在高分化短区域。只有 BWA-MEM 能扩展到人类基因组（nucmer O(n²) 无法扩展）。

4.3 瓶颈分析

关键结论： 短读段（~100bp）瓶颈 = Seeding（FM-index 搜索）；长读段（>1kb）瓶颈 = Banded DP（SW 延伸）。后者是 BWA-MEM2 的优化方向。

4.4 未来改进

SSE2 向量化的 banded DP（BWA-MEM2 实现）
限制 DP 延伸区域，跳过已被长精确匹配覆盖的区域

5. 算法全流程

FASTQ 读段

↓

FMD-index 载入（~5.4 GB）

↓

① Forward Search（5'→3'，确定右端边界）

↓

② Backward Search（3'→5'，找 LEP）

↓

③ LAST 遍历（找 MEM → 过滤 SMEM）

↓

Re-seeding（长 SMEM 中点再播种）

↓

贪心成链 + 链过滤

↓

Banded DP（Z-dropoff + Local/End 抉择）

↓

双端配对（μ,σ² 估计 + mate rescue）

↓

MAPQ 计算

↓

SAM 输出