比对算法总结(三)——BWA-MEM 算法

2021-08-29 本文已影响0人生信小书童

今天要介绍的这篇文献是生物信息领域的大神李恒2013年发表于《BIOINFORMATICS》杂志上的 BWA-MEM 算法的文章，在此之前李恒已经于2009年发表了短序列比对的算法工具 BWA《Fast and accurate short read alignment with Burrows–Wheeler transform》，于2010年发表了长序列比对的算法工具BWA-SW《Fast and accurate long-read alignment with Burrows–Wheeler transform》。BWA-MEM 算法在BWA、BWA-SW基础上又做了诸多改进，同时适用于短reads和长reads的比对分析。自此之后 BWA-MEM 和 Bowite2 一起成为序列比对领域，特别是人基因组重测序领域最常使用的比对工具，它们都属于BWT索引结构的比对算法。

摘要

BWA-MEM 是一种新的比对算法，用于将测序 reads 或者组装后 contigs 比对至大型参考基因组，例如人参考基因组。它会自动选择局部比对和 end-to-end 比对模式，支持PE reads 比对和嵌合体比对。该算法对测序错误有良好的稳定性，适用的reads 长度范围较广，从70bp至几Mb。对于100bp的序列，BWA-MEM算法的表现是目前所有比对软件中最好的。BWA-MEM算法被嵌套入了 BWA 软件，具体代码和说明见 http://github.com/lh3/bwa.

1 引言

NGS测序刚兴起时测试reads长度在36bp左右，在这个阶段有许多比对软件被开发出来。对于36bp的reads ,采用 end-to-end 比对的方法，将reads的所有碱基都比对到参考基因，并给出比对结果和相应的编辑距离，这对于短 reads 的比对是合理的。但随着测序技术的进步，测序的reads变长，对于100bp及以上长度的reads, 这类比对软件不再适合。由于reads 中结构变异的存在以及基因组序列错误组装，需要新的比对算法在空位罚分条件下允许长gap的存在。在此背景下已有一些适用于长reads比对的软件算法被开发出来，包括BWA-SW (Li and Durbin, 2010)、 Bowtie2(Langmead and Salzberg, 2012)、 Cushaw2 (Liu and Schmidt, 2012)、GEM(Marco-Sola et al., 2012), 但这些软件都有相应的弱点：在100bp的reads比对时BWA-SW速度比Bowite2要慢；Bowite2 和 Cushaw2 在600bp的reads 比对时速度较慢；GEM采用了end-to-end 的比对模式并不适用于空位罚分比对，限制了它在长reads比对上的应用。在100bp-1000bp的reads长度范围，没有比对软件在整个区域都有较好的效果。

2 方法

2.1 单端read的比对

2.1.1 种子序列定位及再定位

BWA-MEM 沿用了经典的 种子序列定位及延伸算法（seed and extend ）。初始种子序列的采用了超精确匹配序列（SMEMs）算法 (2012年文章-算法5)，它会在每个匹配位置上找到覆盖该位置的最长精确匹配。但是，有时真正的匹配位置并不含任何的SMEMs。为了减少错误种子序列定位造成的错配，这里采用了种子序列再定位的方法。假如有一个SMEM的长度为l，它在参考基因组中出现的次数是k。如果默认的 l 长度（28bp）太长,我们会重新选取跨越SMEM中心的序列作为种子进行再次定位，前提是它至少在基因组上有k+1次匹配。这些种子序列可以通过SMEM算法中减少匹配次数实现。

2.1.2 链形成与链过滤

作者把一组共线且彼此相邻的种子序列称为链。再构建链的过程中，我们会过滤掉较短的链。过滤短链的标准是它被包含在长链中，并且长度比长链短50%或者38bp。链过滤的目的是为了减少下一步延伸过程的错误。由于形成的链并不一定和最终的匹配相关，因此链的构建过程不需要十分精确。

2.1.3 种子序列延伸

首先根据链的长度进行排序，然后再根据链中种子序列的长度进行排序。排好的种子序列中，去除那些之前有比对位置的种子序列，并用空位罚分的动态规划算法扩展新的序列比对。
BWA-MEM 算法的种子序列延伸与标准的延伸有两方面的不同。第一，假设在延伸步骤，在参考基因组的x位置及查询序列的y位置获得了最高的延伸分数。如果在（x,y）位置的分值与最高的延伸分值的差异大于 Z+|x−y|×PgapExt 的值（Z是设置的阈值，PgapExt 是空位延伸罚分），就终止延伸。这种启发式的算法可以避免造成过度延伸，比对时跨越糟糕的区域。
第二，当进行种子序列延伸时，BWA-MEM会持续记录查询序列延伸至最后一个碱基过程中的最高延伸值。如果延伸至最终碱基的分值与最优局部比对的分值差异低于设定阈值时，拒绝最优局部比对延伸，选择end-to-end延伸。BWA-MEM 使用这种策略自动选择局部比对和end-to-end比对模式。这样就可以保证含有突变的reads的末尾碱基可以比对至参考基因组，减少比对的偏倚，同时可以避免在进行嵌合体序列比对至末尾碱基时引入大量的错配碱基和空位序列。

2.2 双端read的比对

2.2.1 错配比对的找回

与BWA最早版本类似，BWA-MEM 也是同时处理一批reads。在进行每个批次处理时，软件会记录可靠比对比对的插入序列的长度分布，并计算插入序列的平均值 μ 和标准差 σ。对于top 100个候选比对，如果一条reads的配对reads在 [μ - 4σ,μ + 4σ] 范围内是未比对的状态时，BWA-MEM 会采用 SSE2-based Smith-Waterman 算法重新进行比对。这样做的目的是为了找回长重复区域的错误比对。常规的单read 比对以及SSE2-SW 算法找回的比对结果将会用于下一步的配对分析。

2.2.2 配对分析

在考虑正向 read 的第 i 条比对和反向 read 的第 j 条比对时，如果配对reads的方向正确，BWA-MEM 软件会计算配对reads比对的距离 $d_ij$ 。同时记录第 i 条比对和第 j 条比对的得分: $S_ij$ = $S_i$ + $S_j$ - min{ -a $log_4$ P( $d_ij$ ),U }，其中 $S_i$ 和 $S_j$ 是SSE2-SW 算法的比对得分；a 是匹配分值；P(d) 是假设插入序列长度程正态分布时，观测到的插入序列大于比对距离d的概率；U 是控制是否进行配对分析的阈值：如果 $d_ij$ 足够小时， -a $log_4$ P( $d_ij$ ) < U，BWA-MEM 倾向于认为正反向read是比对配对的；否者倾向认为正反向read是比对不配对的。BWA-MEM会选择 $S_ij$ 的最高分值作为最终的配对比对。这种配对联合分析的策略既考虑了比对的分值、也考虑了插入序列的大小和read是嵌合体的可能性。

3 结果

1、用模拟数据评估了NovoAlign,GEM, Bowtie2,Cushaw2, SeqAlto, BWA-SW, BWA 和 BWA-MEM 的表现。在PE和SE reads的比对上，在准确性方面，NovoAlign是最好的，BWA-MEM略低于它，但优于其他软件。在速度方面，reads长度100bp时，BWA-MEM 速度与GEM、Bowite2 相似。在reads长度达到650bp时，BWA-MEM 速度是Bowite2 的6倍。
2、作者认为 BWA-MEM 算法在较长序列比对上有更优的表现，第一基于它先进的种子序列定位方法，第二基于结合了动态规划算法，保证了查询序列长度的线性时间复杂度。
3、BWA-MEM 是一款快速准确的比对软件，它适用于从70bp至几Mb的序列比对。

image.png

4 参考文献

[1] Li H. (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2 [q-bio.GN].