动机:已有短读段比对器对 >200bp 的读段效率低下。长读段场景下仅 BLAT 和 SSAHA2 可用,但按碱基/时间计算远慢于短读段比对器。
结果:提出 BWA-SW(Burrows-Wheeler Aligner's Smith-Waterman Alignment),将长达 1Mb 的序列比对到大型参考数据库。准确率与 SSAHA2 相当,比 BLAT 更准,速度数倍到数十倍更快。
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.
自 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.
BWT-SW 保证找到所有局部匹配,但对长查询不够快。BWA-SW 引入两条启发式(Z-best + 种子过滤)来再次加速。这是典型的"用精度换速度"权衡——在长读段含有大量冗余信息的情况下非常合理。
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.
前缀树中不同节点可能有相同 SA 区间。DAWG 合并这些节点,使节点与 SA 区间一一对应——一个 DAWG 节点可代表多个前缀关系的子串。
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.
标准 SW 遍历全参考序列(3Gb × 1kb 不可行)。前缀 DAWG 将查询中相同子串合并为单节点。DP 时若 u 不匹配 v,则 u 也不匹配 v 的任何后代→只需访问靠近 𝒯 根的节点,省去全基因组扫描。
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.
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 区间 ≤ 3 → 独特 → 值得 SW。SA 区间大 → 重复 → 不触发 SW。BWT-SW 不检查 SA 区间大小,重复序列触发大量无意义 SW——这是 BWA-SW 快于 BWT-SW 的主因。
即使 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.
真实比对通常含多个种子。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.
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.
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.
~10M bp,100bp–10kb,错误率 2%–10%。
| 程序 | 100bp | 200bp | 500bp | 1000bp | 10000bp | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2% | 5% | 10% | 2% | 5% | 10% | 2% | 5% | 10% | 2% | 5% | 10% | 2% | 5% | 10% | |
| CPU 时间 (s) | |||||||||||||||
| BLAT | 685 | 577 | 559 | 819 | 538 | 486 | 1078 | 699 | 512 | 1315 | 862 | 599 | 2628 | 1742 | 710 |
| BWA-SW | 165 | 125 | 84 | 222 | 168 | 118 | 249 | 172 | 152 | 234 | 168 | 150 | 158 | 134 | 120 |
| SSAHA2 | 4872 | 7962 | 9345 | 1932 | 2236 | 5252 | 3311 | 8213 | 6863 | 1554 | 1583 | 3113 | – | – | – |
| Q20% / Err% | |||||||||||||||
| BLAT Q20% | 68.7 | 25.5 | 3.0 | 92.0 | 52.9 | 7.8 | 97.1 | 86.3 | 21.4 | 97.7 | 96.4 | 39.0 | 98.4 | 99.0 | 94.0 |
| BLAT Err% | 0.99 | 2.48 | 5.47 | 0.55 | 1.72 | 4.55 | 0.17 | 1.12 | 4.41 | 0.01 | 0.52 | 3.98 | 0.00 | 0.00 | 1.28 |
| BWA-SW Q20% | 85.1 | 62.2 | 19.8 | 93.8 | 88.7 | 49.7 | 96.1 | 95.5 | 85.1 | 96.9 | 96.5 | 95.0 | 98.4 | 98.5 | 98.1 |
| BWA-SW Err% | 0.01 | 0.05 | 0.17 | 0.00 | 0.02 | 0.13 | 0.00 | 0.00 | 0.04 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| SSAHA2 Q20% | 85.5 | 83.8 | 78.2 | 93.4 | 93.1 | 91.9 | 96.6 | 96.5 | 96.1 | 97.7 | 97.6 | 97.4 | – | – | – |
| SSAHA2 Err% | 0.00 | 0.01 | 0.19 | 0.01 | 0.00 | 0.01 | 0.00 | 0.01 | 0.04 | 0.00 | 0.00 | 0.00 | – | – | – |
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.
两个测试嵌合体(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.
种子策略: 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 系列第一个 seed-and-extend 算法。它证明了两件事:① BWT 索引上做启发式 DP 可以比 BLAT/SSAHA2 更快更准;② SA 区间截断能有效抑制重复区域的 SW 调用。
BWA-MEM 后来沿用了这一框架,但将种子搜索替换为 FMD-index 上的 SMEM 搜索——更快更简洁。