基于 BWT 的快速准确长读段比对

Fast and accurate long-read alignment with Burrows-Wheeler transform
Heng Li, Richard Durbin · Wellcome Trust Sanger InstituteBioinformatics 2010; 26(5): 589–595DOI: 10.1093/bioinformatics/btp698

摘要

动机:已有短读段比对器对 >200bp 的读段效率低下。长读段场景下仅 BLAT 和 SSAHA2 可用,但按碱基/时间计算远慢于短读段比对器。

结果:提出 BWA-SW(Burrows-Wheeler Aligner's Smith-Waterman Alignment),将长达 1Mb 的序列比对到大型参考数据库。准确率与 SSAHA2 相当,比 BLAT 更准,速度数倍到数十倍更快。

获取: BWA 组件,GPL 开源。bio-bwa.sourceforge.net
原文

Motivation: Many programs for aligning short sequencing reads to a reference genome have been developed in the last 2 years. Most of them are very efficient for short reads but inefficient or not applicable for reads >200 bp because the algorithms are heavily and specifically tuned for short queries with low sequencing error rate. For longer reads, hashing-based software such as BLAT and SSAHA2 remain the only choices. Nonetheless, these methods are substantially slower than short-read aligners in terms of aligned bases per unit time.

Results: We designed and implemented a new algorithm, Burrows-Wheeler Aligner's Smith-Waterman Alignment (BWA-SW), to align long sequences up to 1 Mb against a large sequence database (e.g. the human genome) with a few gigabytes of memory. The algorithm is as accurate as SSAHA2, more accurate than BLAT, and is several to tens of times faster than both.

1. 引言

自 2000 年起,从 MegaBLAST、SSAHA2、BLAT 到 PatternHunter,DNA 匹配速度持续提升。短读段时代又催生了 SOAP、MAQ、Bowtie、BWA。但长读段(Roche/454 >400bp、Illumina >100bp、PacBio >1000bp)带来不同需求:① 偏好局部比对(结构变异不打断比对,读尾端错配影响小);② 必须容忍频繁 indel(454/PacBio 错误主要来自 indel)。

Lam et al.(2008)的 BWT-SW 通过 FM-index 在后缀树上做 SW,准确率与标准 SW 相同,数千倍更快。但长查询仍慢于 BLAST。BWA-SW 在 BWT-SW 基础上引入启发式加速。

原文

Long-read alignment has different objectives from short-read alignment. In long-read alignment, we would prefer to find local matches because a long read is more fragile to structural variations and misassemblies in the reference. Lam et al. (2008) furthered this idea by implicitly representing the suffix tree with an FM-index to achieve a small memory footprint. Their new algorithm, BWT-SW, is able to deliver identical results to the standard Smith–Waterman alignment, but thousands of times faster. While BWT-SW is still slower than BLAST on long query sequences, it finds all matches without heuristics. Our BWA-SW algorithm follows this route.

BWA-SW = BWT-SW 的启发式加速版

BWT-SW 保证找到所有局部匹配,但对长查询不够快。BWA-SW 引入两条启发式(Z-best + 种子过滤)来再次加速。这是典型的"用精度换速度"权衡——在长读段含有大量冗余信息的情况下非常合理。

2. 方法

2.1 算法概览

BWA-SW 为参考和查询分别建 FM-index,在两者间执行 DP。参考 → 前缀树 𝒯(X);查询 → 前缀 DAWG 𝒢(W)(合并相同 SA 区间节点的 DAG)。

两条启发式:① Z-best——每查询节点只保留 Z 个最高分参考匹配;② 非重叠报告——查询上不重叠的比对优先。

原文

BWA-SW builds FM-indices for both the reference and query sequence. It implicitly represents the reference sequence in a prefix trie and represents the query sequence in a prefix directed acyclic word graph (prefix DAWG). A dynamic programming can be applied between the trie and the DAWG.

2.2 前缀树与前缀 DAWG

前缀树中不同节点可能有相同 SA 区间。DAWG 合并这些节点,使节点与 SA 区间一一对应——一个 DAWG 节点可代表多个前缀关系的子串。

图 1:前缀树 vs 前缀 DAWG("GOOGOL") (A) 前缀树 (B) 前缀 DAWG G[1,3] O[4,6] GO[1,2] GL[3,3] GOO[1,1] GOL[1,1] OG[4,4] OO[5,5] OL[6,6] OGO[4,4] OGOL[4,4] [1,3] [4,6] [1,2] [3,3] [4,4] [5,5] [6,6] [1,1] DAWG 节点 [4,4] 代表 {OG, OGO, OGOL} [1,1] 被 [1,2] 和 [4,4] 同时指向(两条入边)
图 1 (A) 前缀树与 (B) 前缀 DAWG 的对比。DAWG 合并相同 SA 区间的节点,[4,4] 代表三个子串 OG/OGO/OGOL。
原文

The prefix DAWG of X is transformed from the prefix trie by collapsing nodes having an identical interval. Thus in the prefix DAWG, nodes and SA intervals have an one-to-one relationship, and a node may represent multiple substrings of X, falling in a sequence where each is a prefix of the next.

为什么用 DAWG?

标准 SW 遍历全参考序列(3Gb × 1kb 不可行)。前缀 DAWG 将查询中相同子串合并为单节点。DP 时若 u 不匹配 v,则 u 也不匹配 v 的任何后代→只需访问靠近 𝒯 根的节点,省去全基因组扫描。

2.3 前缀树与 DAWG 的 DP

Guv=max(0,Gu′v′+S(u′,u;v′,v),Iuv,Duv) (1)
Iuv=max(Gu′v−q,Iu′v−r) (2) Duv=max(Guv′−q,Duv′−r) (3)

u′ 为 u 父节点,v′ 为 v 父节点,q=gap open,r=gap extension。

原文

The DP is performed by traversing both 𝒢(W) and 𝒯(X) in reverse post-order. Noting that once u does not match v, u does not match any nodes descending from v, we only need to visit the nodes close to the root.

2.4 SW 加速与种子截断

DP 常数大,用 DP 找种子区间对(Guv 够高且 SA 区间 ≤ 3),再由标准 SW 延伸。1kb 比对约 10–20 种子。与 BWT-SW 区别:BWT-SW 分数够高即触发 SW,不管 SA 区间大小。

原文

When Guv is good enough and the SA interval size of v is below a certain threshold (3 by default), we save the (u,v) pair, called a seed interval pair. Typically for 1 kb alignments there will be 10–20 seeds.

SA 区间截断的直觉

SA 区间 ≤ 3 → 独特 → 值得 SW。SA 区间大 → 重复 → 不触发 SW。BWT-SW 不检查 SA 区间大小,重复序列触发大量无意义 SW——这是 BWA-SW 快于 BWT-SW 的主因。

2.5 启发式加速

Z-best 策略

即使 500bp 唯一查询,也可能数百万 𝒯(X) 节点正分匹配。每 u 只保留前 Z 个最高分。Z=1 对 <5% 错误率的 200bp 读段已够。辅以反向-反向比对降低漏检率。

原文

At each node u in 𝒢(W) we only keep the top Z best scoring nodes in 𝒯(X). We find even Z=1 works well with high-quality 200 bp reads (<5% sequencing error rate). We also align the reverse query sequence to the reverse reference sequence.

Z-best 风险控制

真实比对通常含多个种子。P(全部被剪枝) ≈ (1/Z)^种子数。Z=10、种子数=10 → 10⁻¹⁰。长比对比短比对更安全。

种子过滤(链化)

带宽 50bp 区域内的种子组成链。短链被长链完全覆盖且种子数 < 长链的 1/10 → 丢弃。10kb 模拟数据时间减半,准确率无损失。

原文

We chain seeds that are contained in a band (default band width 50 bp). If a short chain is fully contained in a long chain and the number of seeds in the short chain is below one-tenth of the number of seeds in the long chain, we discard all the seeds in the short chain.

2.6 比对质量

MAPQ = 250 · c₁ · c₂ · (S₁−S₂)/S₁ (4)

c₁=覆盖>4种子?1:0.5;c₂=正+反都找到?1:0.2。

原文

To estimate the mapping quality of a BWA-SW alignment, we fit an empirical formula: 250·c1·c2·(S1−S2)/S1.

3. 结果

3.1 实现

BWA 组件,输入 FASTA/Q + BWA 索引,输出 SAM。FM-index 约 3.7GB,典型 <4GB,1Mb 查询峰值 6.4GB。支持多线程,自动调参。

原文

Memory usage is dominated by the FM-index, about 3.7 GB for the human genome. On typical sequencing reads, the total memory is <4 GB.

3.2 模拟数据

~10M bp,100bp–10kb,错误率 2%–10%。

程序100bp200bp500bp1000bp10000bp
2%5%10%2%5%10%2%5%10%2%5%10%2%5%10%
CPU 时间 (s)
BLAT6855775598195384861078699512131586259926281742710
BWA-SW16512584222168118249172152234168150158134120
SSAHA2487279629345193222365252331182136863155415833113
Q20% / Err%
BLAT Q20%68.725.53.092.052.97.897.186.321.497.796.439.098.499.094.0
BLAT Err%0.992.485.470.551.724.550.171.124.410.010.523.980.000.001.28
BWA-SW Q20%85.162.219.893.888.749.796.195.585.196.996.595.098.498.598.1
BWA-SW Err%0.010.050.170.000.020.130.000.000.040.000.000.000.000.000.00
SSAHA2 Q20%85.583.878.293.493.191.996.696.596.197.797.697.4
SSAHA2 Err%0.000.010.190.010.000.010.000.010.040.000.000.00
表 1 2.5 GHz Xeon E5420 单核。BWA-SW 速度是 BLAT 的 2–10×、SSAHA2 的 10–60×。
原文

Table 1 shows the CPU time, fraction of confidently aligned reads and alignment error rates for BLAT, BWA-SW and SSAHA2 given different read lengths and error rates. BWA-SW is clearly the fastest, several times faster than BLAT and SSAHA2 on all inputs, and its speed is not sensitive to the read length or error rates.

表 1 解读
  • 速度:BWA-SW 是 BLAT 的 2–10×,SSAHA2 的 10–60×。100bp→10kb 时间几乎不变(165s→120s)。
  • 准确率:≥500bp 时 BWA-SW Q20% 85–98%,错误率接近 0。BLAT 在 10% 错误率时不可用(Q20% 3–39%)。
  • BWA-SW 速度对读段长度和错误率几乎不敏感——这正是设计目标。

3.3 嵌合体与真实数据

两个测试嵌合体(1000bp+300bp):BWA-SW 全部正确;SSAHA2 找不到 300bp 片段。454 数据(355bp 平均):BWA-SW 错过 993 条 SSAHA2 正确比对的短读段;SSAHA2 错过 1,188 条 BWA-SW 正确比对的长读段。

原文

When we align the two chimeric reads, BWA-SW reports four alignments. The latest SSAHA2 fails to find the alignment of the 300 bp pieces. On real data (SRR003161), BWA-SW misses 993 alignments which SSAHA2 aligns well, while SSAHA2 misses 1,188.

4. 讨论

种子策略: BLAT/SSAHA2 用固定 11–12bp 精确 k-mer(~10⁵ 种子/次),BWA-SW 用 DP 生成变长种子(~10–20 种子/kb)。

与 BWT-SW 的区别: ① 完整性(BWT-SW 保证全部找到);② 对齐两个 FM-index;③ 内层遍历参考前缀树使 Z-best 可行;④ 只报告不重叠比对。

原文

BWA-SW, BLAT and SSAHA2 all follow the seed-and-extend paradigm. The major difference comes from the seeding strategy. BLAT and SSAHA2 identify short exact matches as seeds, typically of length 11 or 12 bp. BWA-SW resolves this issue by using a few long gapped seeds in unique regions.

BWA-SW 的历史定位

BWA-SW 是 BWA 系列第一个 seed-and-extend 算法。它证明了两件事:① BWT 索引上做启发式 DP 可以比 BLAT/SSAHA2 更快更准;② SA 区间截断能有效抑制重复区域的 SW 调用。

BWA-MEM 后来沿用了这一框架,但将种子搜索替换为 FMD-index 上的 SMEM 搜索——更快更简洁。