摘要

动机:在过去两年里,已经开发了许多将短测序读数与参考基因组对齐的程序。大多数算法对于短读非常有效,但对于大于200 bp的读来说效率很低或不适用,因为这些算法针对排序错误率低的短查询进行了大量专门调整。然而,一些测序平台已经产生了更长的读数,其他平台预计很快就会上市。对于较长时间的读取,基于散列的软件(如BLAT和SSAHA2)仍然是唯一的选择。尽管如此,就每单位时间对齐的碱基而言,这些方法要比短读校准器慢得多。

结果:我们设计并实现了一种新的算法,即Burrows-Wheeler Aligner的Smith-Waterman Alignment(BWA-SW)算法,该算法可以将长至1Mb的序列与具有数GB内存的大型序列数据库(例如人类基因组)对齐。该算法与SSAHA2一样精确,比BLAT更精确,并且比两者都快几到几十倍。

可用性: http://bio-bwa.sourceforge.net

联系人: rd@sanger.ac.uk

1简介

在开发敏感的局部对准软件后,如FASTA(Pearson和Lipman,1988)和BLAST(Altschul等。,1997)大约在1990年,自2000年以来,开发了新一代更快的方法来寻找DNA序列匹配,包括MegaBLAST(Morgulis等。,2008; 等。,2000)、SSAHA2(宁等。,2001),BLAT(肯特,2002)和PatternHunter(马)等。,2002),大大加快了毛细管测序与大参考基因组的匹配速度。当新的测序技术问世,产生了数百万次短(<100 bp)读取时,各种新算法被开发出来,速度快了10到1000倍,其中包括SOAP(Li,R。等。,2008),MAQ(Li,H。等。,2008)、鲍蒂(兰米德等。,2009)和BWA(Li和Durbin,2009). 然而,Roche/454测序技术已经在生产中产生了>400bp的读数,Illumina逐渐增加了>100bp的读数长度,太平洋生物科学公司在早期测试中产生了1000bp的读数(Eid等。,2009). 来自新测序技术的读数不再短,这有效地排除了许多专为读数不超过100 bp的新对准器。有效地将长读数与人类基因组等长参考序列进行比对,对比对工具的开发提出了新的挑战。

长读对齐与短读对齐的目标不同。首先,在短读对齐中,我们通常会对齐全长读取,以减少由于接近读取末尾的不匹配而导致的参考偏差。考虑到这一要求,我们可以设计间隔种子模板(Ma等。,2002)跨越整个阅读范围(江和黄,2008; 等。,2008; 史密斯等。,2008),或快速筛选出较差的匹配项,例如,通过应用q-gram过滤(隆隆声等。,2009; 杂草等。,2009)或者通过限定搜索过程(Li和Durbin,2009),从而加速对齐。然而,在长读对齐中,我们更喜欢找到本地匹配,因为长读更容易受到引用中结构变化和错误组装的影响,但受接近读取末尾的不匹配影响较小。其次,许多短读对准器只有在进行未映射对准或允许有限间隙(例如最大一个间隙)时才有效。他们找不到更多的差距,或者当他们为此任务进行调整时,性能会迅速下降。然而,长读比对仪必须允许存在比对间隙,因为长读中indels出现的频率更高,并且可能是某些技术(如454和Pacific Bioscience)测序错误的主要来源。

当考虑使用算法来加速长读对齐时,大多数当前软件中使用的散列表索引并不是唯一的选择。米克等。(2003)发现了一种类似Smith–Waterman的动态编程,可以应用于查询序列和引用的后缀树之间,通过自顶向下的遍历有效地将查询与后缀树中采样的每个子序列对齐。在后缀树上,相同的序列折叠在一条路径上,通过避免相同子序列的重复对齐来节省时间。等。(2008)通过用FM-index隐式表示后缀树(Ferragina和Manzini,2000)基于Burrows–Wheeler变换(BWT;BurrowsandWheeler,1994),以实现较小的内存占用。他们的新算法BWT-SW能够提供与标准Smith–Waterman比对相同的结果,但与人类基因组序列比对时速度快数千倍。虽然BWT-SW在长查询序列上仍比BLAST慢,但它可以在没有启发式的情况下查找所有匹配项。可以想象,引入启发式将进一步加速BWT-SW。我们的BWA-SW算法遵循这一路线。

在某种程度上,BWA-SW和BWT-SW也遵循种子和扩展范式。但与BLAT和SSAHA2不同,BWA-SW通过动态编程在两个FM诱导之间查找种子,并允许种子中存在不匹配和间隙。当种子在引用序列中很少出现时,它会扩展种子。通过减少高度重复序列的不必要扩展来获得速度。

在本文中,我们将描述用于长读取比对的新比对算法BWA-SW,并在模拟数据和实际数据上评估其实际性能以及BLAT和SSAHA2。我们还将简要介绍后缀数组和FM-index,但读者可以参考Li和Durbin(2009)了解更多详细信息。

2方法

2.1 BWA-SW算法概述

BWA-SW为引用序列和查询序列构建FM-indices。它隐式表示前缀trie中的引用序列,并表示前缀定向非循环单词图(前缀DAWG;Blumer等。,1985),它是从查询序列的前缀trie转换而来的(第2.3节). 通过分别遍历引用前缀trie和查询DAWG,可以在trie和DAWG之间应用动态编程。如果不应用启发式,这种动态编程将找到所有本地匹配,但不会比BWT-SW快。在BWA-SW中,我们应用了两个启发式规则来大大加快这一过程。首先,查询DAWG的遍历是在外循环中进行的,因此,在没有完成动态编程的情况下,我们知道引用前缀trie中与查询节点匹配的所有节点都是正分数。基于观察到真正的对齐往往具有较高的对齐分数,我们可以在每个节点处剪除低分数匹配,以限制仅围绕好匹配进行动态规划。因此,可以大大减少动态编程的规模。在此过程中,可以修剪真正的对齐,但实际上,这可以通过使用启发式进行控制,并且在给定长或高质量查询序列的情况下很少发生。其次,BWA-SW只报告查询序列上基本上不重叠的对齐,而不提供所有重要的局部对齐。它启发式地识别并丢弃较长路线中包含的种子,从而节省了不成功种子扩展的计算时间。

2.2符号和定义

2.2.1后缀数组和BWT

设∑={A类,C类,,T型}是核苷酸字母表,$是一个比∑中所有符号都小的符号。给定核苷酸序列X(X)=1n个−1具有n个−1=$,让X(X)[]=成为-第个符号,X(X)[,j个]=j个的子序列X(X)X(X)=X(X)[,n个−1]后缀X(X).后缀数组S公司属于X(X)是整数0,…,的置换,…,n个−1,以便S公司()=j个当且仅当X(X)j个-第th个按字典顺序排列的最小后缀。的BWTX(X)是的排列X(X),其中B类[]=$如果S公司()=0和B类[]=X(X)[S公司()−1]否则。

2.2.2后缀数组间隔

给定一个序列W公司,的后缀数组间隔SA间隔论坛属于W公司定义为
特别是,如果W公司是空字符串,R(右)(W公司)=1和论坛。所有事件的位置集W公司论坛.
C类()=#{0≤j个n个−2:X(X)[j个]<}和O(运行)(,)=#{0≤j个:B类[j个]<},其中{·}计算集合的基数(或大小)。费拉吉纳和曼奇尼(2000)证明了
还有那个论坛当且仅当是的子字符串X(X).

2.2.3 FM诱导

后缀数组S公司,数组C类O(运行)足以精确搜索中的模式X(X)FM-index(费拉吉纳和曼奇尼,2000)是三个数组的压缩表示,由压缩的BWT字符串组成B类,用于计算的辅助数组O(运行),和后缀数组的一部分S公司然而,BWA-SW使用了一个简化的FM索引,其中我们不压缩B类并存储引用数组的一部分O(运行)没有辅助数据结构。简化版对于字母表非常小的DNA序列更有效。有关施工的详细信息,请参阅我们之前的论文(Li和Durbin,2009).

2.2.4对准

对齐是元组(W公司1,W公司2,A类)其中W公司1W公司2是两个字符串和A类是一系列复制、替换、插入和删除操作W公司2进入之内W公司1。插入和删除是间隙.间隙和替换为差异. The编辑距离对齐的数量等于A类.

给定计分系统,可以计算路线的分数。我们说W公司1比赛W公司2如果W公司1W公司2可以与正分对齐,在这种情况下,我们还可以说(W公司1,W公司2)是一场比赛。

一场比赛(W公司1,W公司2)据说是包含英寸(W公司1,W公司2)关于第一个序列如果W公司1是的子字符串W公司1类似地,我们可以定义对齐(更强的条件)之间以及对齐和匹配之间的“包含”关系。

2.3前缀trie和前缀DAWG

这个前缀trie字符串的X(X)是一棵树,每条边都标有一个符号,这样从叶子到根的路径上的符号串联起来就有了一个唯一的前缀X(X)。从节点到根的边缘符号串联始终是X(X),调用节点表示的字符串。这个SA间隔节点的定义为节点表示的字符串的SA间隔。不同的节点可能具有相同的间隔,但回顾SA间隔的定义,我们知道这些节点表示的字符串必须是同一字符串的前缀,并且具有不同的长度。

这个前缀DAWG第,页,共页X(X)通过折叠具有相同间隔的节点,从前缀trie进行转换。因此,在前缀DAWG中,节点和SA间隔具有一对一的关系,一个节点可以表示X(X),如前一段所述,每一个都是下一个的前缀。图1给出了一个示例。

字符串“GOOGOL”的前缀trie和前缀DAWG。(A) 前缀trie。符号“∧”标记字符串的开始。节点中的两个数字表示节点的SA间隔。(B) 通过折叠具有相同SA间隔的节点来构造前缀DAWG。例如,在前缀trie中,三个节点具有SA间隔[4,4]。他们的父母分别有间隔[1,2],[1,2]和[1,1]。在前缀DAWG中,[4,4]节点因此具有父节点[1,2]和[1,1]。节点[4,4]表示三个字符串“OG”、“OGO”和“OGOL”,前两个字符串是“OGOL”的前缀。(A) 修改自Li和Durbin(2009)的图1。
图1。

字符串的前缀trie和前缀DAWGGOOGOL公司’. (A类)前缀trie。符号“∧”表示字符串的开始。节点中的两个数字表示节点的SA间隔。(B类)通过折叠具有相同SA间隔的节点来构造前缀DAWG。例如,在前缀trie中,三个节点具有SA间隔[4,4]。他们的父母分别有间隔[1,2],[1,2]和[1,1]。在前缀DAWG中,[4,4]节点因此具有父节点[1,2]和[1,1]。节点[4,4]表示三个字符串OG公司’, ‘OGO公司'和'OGOL公司'前两个字符串是'的前缀OGOL公司’. (A) 修改自图1在李和杜宾(2009).

2.4根据前缀DAWG对齐前缀trie

我们构造了一个前缀DAWG𝒢(W公司)用于查询序列W公司和前缀trie𝒯(X(X))供参考X(X).计算最佳得分的动态编程W公司X(X)如下所示。紫外线=紫外线=D类紫外线=0时u个是𝒢的根(W公司)以及v(v)\119983;的根(X(X)). 在节点处u个在𝒢中(W公司),对于其每个父节点u个',计算
哪里v(v)'是的父级v(v)在\119983年(X(X)),函数S公司(u个′,u个;v(v)’,v(v))给出边缘符号之间的分数(u个′,u个)还有那个(v(v)′,v(v))、和q个第页分别是空位罚分和空位延长罚分。紫外线,紫外线D类紫外线计算公式如下:
其中pre(u个)是的父节点集u个.紫外线等于表示的(可能是多个)子字符串之间的最佳分数u个和由v(v)。我们说节点v匹配u如果紫外线>0.

通过遍历两个𝒢来执行动态编程(W公司)和𝒯(X(X))以相反的后序(即所有父节点在子节点之前访问)嵌套方式。注意到这一次u个不匹配v(v),u个不匹配从以下降序的任何节点v(v),我们只需要访问靠近\119983;根的节点(X(X))无需遍历整个trie,与总是遍历整个参考序列的标准Smith–Waterman算法相比,这大大减少了迭代次数。

2.5标准Smith–Waterman加速

与时间复杂度为的标准Smith–Waterman比对相比O(运行)(|X(X)|·|W公司|),BWA-SW具有更好的时间复杂性,因为它不比时间复杂性的BWT-SW慢O(运行)(|X(X)|0.628|W公司|)(林等。,2008). 得出这个结论是因为对于短期次级比对,我们正在考虑使用单个匹配项进行多个可能的匹配紫外线比较。然而,由于与遍历前缀trie和前缀DAWG相关的复杂步骤,与每个迭代相关的常量要大得多,这使得当我们使用BWA-SW扩展唯一对齐时,BWA-SF效率低下。更有效的策略是使用BWA-SW查找部分匹配,并应用Smith–Waterman算法进行扩展。在动态编程中,我们知道在任何对中考虑的部分匹配数,因为这可以从SA间隔的大小计算得出。什么时候?紫外线足够好,SA间隔大小为v(v)低于某个阈值(默认为3),则我们保存(u个,v(v))一对,称为种子间隔对,不要深入v(v)\119983;中的节点(X(X)). 通过查找后缀数组X(X)W公司,我们可以推导种子匹配,或者简单地说种子,来自种子间隔对。随后,这些种子将由Smith–Waterman算法进行扩展。如果整个查询是一个高度重复的序列,那么它将完全与上一节中描述的算法对齐,而不使用Smith–Waterman扩展。

因为我们正在尽早停止动态编程以生成种子,所以全局最佳对齐可能包含多个种子,实际上,对于长对齐,情况往往是这样的。通常,对于1 kb的比对,将有10到20个种子。下面我们将利用这一观察来启发性地加快搜索。

BWT-SW在执行序列和前缀trie之间的动态编程时部署了类似的策略,以查找种子匹配,然后执行Smith–Waterman扩展。与我们算法的主要区别是,无论SA间隔大小如何,只要分数足够高,BWT-SW就会启动Smith–Waterman对齐。有时,重复序列可能与人类基因组中的数千个位置相匹配,每次扩展部分匹配可能会很慢。

2.6启发式加速

2.6.1 Z-best策略

到目前为止描述的算法是精确的,因为它能够提供与Smith–Waterman算法相同的结果。尽管在给定长参考序列的情况下,它比标准算法快得多,但它还不够快,无法对齐大规模测序数据。更仔细的调查表明,即使是一个唯一的500 bp查询序列,在𝒯中也有几百万个节点(X(X))可以使用正对齐分数匹配查询。这些节点中的大多数是随机匹配或在较短的低复杂度区域中的匹配。访问所有这些网站都是浪费。

为了加速对齐,我们穿过𝒢(W公司)在外环和𝒯中(X(X))在内部循环中,以及在每个节点u个在𝒢中(W公司)我们只保留顶部Z轴\119983;中的最佳得分节点(X(X))那场比赛u个,而不是保留所有匹配的节点。这种启发式策略称为Z最佳当然,当我们应用Z轴-最好的策略是,当一场错误的比赛得分较高时,我们可能会错过真实路线中包含的种子。但是,如果查询与引用几乎相同,则这种情况发生的频率会降低。此外,如果真正的排列很长并且包含许多种子,那么所有种子为假的可能性很小。模拟数据和实际数据(第3节),我们发现不相上下Z轴=1适用于高质量的200 bp读取(测序错误率<5%)。增加的Z轴到10或更高会略微提高精度,但会大大降低对准速度。

为了减少对齐错误,除了正向-正向对齐外,我们还将反向查询序列与反向引用序列对齐,即反向-反向对齐。理想情况下,正向-正向和反向-反向对齐应该产生相同的结果,但如果真正对齐中的种子具有低得分后缀(或前缀),则正向-正向(或反向-反向)对齐很可能会错过它,而合并两轮对齐会降低机会。此外,如果前向-前向对齐中的最佳对齐包含多个种子匹配,那么它出错的可能性也很小。在实现中,如果最佳对齐默认包含5个或更多种子,则我们不应用反向-反向对齐。

2.6.2在Smith–Waterman扩展之前过滤种子

与BLAST一样,BLAT和SSAHA2都报告了所有重要的比对,或者通常是几十个得分最高的比对,但这不是读取映射中最理想的输出。我们通常更感兴趣的是最佳对齐或最佳少数对齐,涵盖查询序列的每个区域。例如,假设一个1000 bp的查询序列由一条染色体上的900 bp片段和另一条染色体的100 bp片段组成;900 bp片段中的400 bp是一个高度重复的序列。对于BLAST,要知道这是一种嵌合阅读,我们需要要求它报告400 bp重复序列的所有比对,这既昂贵又浪费,因为一般来说,我们对包含在较长唯一序列中的短重复序列的比对不感兴趣。在这个例子中,一个有用的输出是报告900 bp和100 bp段的每一个比对,并指出这两个段是否具有良好的次优比对,这可能会导致最佳比对不可靠。这样的输出简化了下游分析,节省了重建重复序列详细比对的时间。

在BWA-SW中,我们说两条路线是不同的如果查询上重叠区域的长度小于较短查询段长度的一半。我们的目标是找到一组不同的比对,使该组中每个比对的得分之和最大化。这个问题可以通过动态编程解决,但在我们的例子中,读取通常是完全对齐的,贪婪近似很好。在实际实现中,我们根据对齐分数对局部对齐进行排序,从最好的列表中扫描排序后的列表,如果它与所有分数较大的保持对齐不同,则保持对齐;if对齐2被拒绝,因为它与1,我们认为2是次优路线1并使用此信息近似映射质量(第2.7节).

因为我们只保留在查询序列上基本上不重叠的对齐,所以我们不妨丢弃对最终对齐没有贡献的种子。在Smith–Waterman扩展和在不必要的扩展上花费的时间得以节省之前,可以使用另一种启发式方法检测此类种子。为了识别这些种子,我们将包含在带中的种子链接起来(默认带宽为50 bp)。如果在查询序列中,短链完全包含在长链中,并且短链中的种子数低于长链中种子数的十分之一,我们将丢弃短链中所有的种子,因为在这种情况下,短链很少能比长链更好地对齐。Z轴-最好的策略是,这种启发式方法对对准精度没有明显影响。在1000个10 kb的模拟数据上,它将运行时间减半,而精度没有降低。

2.7近似绘图质量

李,H。等。(2008)引入映射质量的概念来估计查询序列被放置在错误位置的概率。如果对齐算法保证找到所有局部对齐,则映射质量仅由这些局部对齐决定。然而,当BWA-SW部署启发式规则时,产生错误对齐的可能性也与启发式有关。为了估计BWA-SW路线的绘图质量,我们拟合了一个经验公式:250·c(c)1·c(c)2·(S公司1S公司2)/S公司1,其中S公司1是最佳对齐的分数,S公司2第二名最佳阵容的得分,c(c)1如果对齐覆盖四个以上种子,则等于1,否则等于0.5,以及c(c)2如果通过正向-正向和反向-反向对齐找到最佳对齐,则等于1,否则为0.2。

3结果

3.1实施

BWA-SW算法是作为BWA程序的一部分实施的(Li和Durbin,2009)根据GNU通用公共许可证(GPL)分发。该实现以BWA索引和查询FASTA或FASTQ文件作为输入,并以SAM格式输出对齐(Li等。,2009). 查询文件通常包含许多序列(读取)。我们依次处理每个查询序列,如果适用,使用多个线程。内存使用主要由FM-index控制,约3.7GB用于人类基因组。每个查询所需的内存大致与序列长度成正比。在典型的顺序读取中,总内存<4 GB;在一个具有100万个碱基对(Mbp)的查询序列上,峰值内存总计为6.4GB。

在实现中,我们尝试根据读取长度和排序错误率自动调整参数,以使默认设置适用于不同特征的输入。对于不熟悉算法的用户来说,这种行为很方便,并且在读取混合长度和错误率的情况下有助于提高性能。

3.2模拟数据评估

在模拟数据上,我们通过比对知道了正确的染色体坐标,评估也很简单。

3.2.1总体性能

表1显示了给定不同读取长度和错误率的BLAT(v34)、BWA-SW(版本0.5.3)和SSAHA2(版本2.4)的CPU时间、可靠对齐读取的分数和对齐错误率。除非必要,否则我们尝试使用每个对齐器的默认命令行选项。根据输入数据的特性微调选项可能会产生更好的性能。

表1。

模拟数据评估

程序韵律学100个基点
200个基点
500个基点
1000个基点
10 000基点
2%5%10%2%5%10%2%5%10%2%5%10%2%5%10%
BLAT(爆炸)CPU秒6855775598195384861078699512131586259926281742710
问题20%68.725.539252.97.897.186.321.497.796.43998.49994
误差%0.992.485.470.551.724.550.171.124.410.010.523.98001.28
体重-体重CPU秒16512584222168118249172152234168150158134120
问题20%85.162.219.893.888.749.796.195.585.196.996.59598.498.598.1
误差%0.010.050.1700.020.13000.04000000
SSAHA2型CPU秒487279629345193222365252331182136863155415833113
问题20%85.583.878.293.493.191.996.696.596.197.797.697.4
误差%00.010.190.0100.0100.010.04000
程序韵律学100个基点
200个基点
500个基点
1000个基点
10 000基点
2%5%10%2%5%10%2%5%10%2%5%10%2%5%10%
BLAT(爆炸)CPU秒6855775598195384861078699512131586259926281742710
问题20%68.725.539252.97.897.186.321.497.796.43998.49994
误差%0.992.485.470.551.724.550.171.124.410.010.523.98001.28
BWA-SW公司CPU秒16512584222168118249172152234168150158134120
问题20%85.162.219.893.888.749.796.195.585.196.996.59598.498.598.1
误差%0.010.050.1700.020.13000.04000000
SSAHA2型CPU秒487279629345193222365252331182136863155415833113
问题20%85.583.878.293.493.191.996.696.596.197.797.697.4
误差%00.010.190.0100.0100.010.04000

从人类基因组中模拟了大约10 000 000 bp不同读取长度和错误率的数据。百分之二十的误差是indel误差,indel长度由几何分布得出(密度:0.7·0.3−1). 这些模拟读数分别与BLAT(选项-fastMap)、BWA-SW和SSAHA2(选项-454表示100和200 bp读数)的人类基因组对齐。然后将对齐的坐标与模拟坐标进行比较,以找出对齐误差。在该表中的每个单元中,这三个数字是Intel E5420 2.5 GHz CPU单核上的CPU秒数、映射质量大于或等于20(Q20)的对齐百分比以及Q20对齐中的错误对齐百分比。SSAHA2和BWA-SW报告制图质量;BLAT映射质量估计为最佳和次最佳对齐分数之差除以最佳对齐分数的250倍(基本上与BWA-SW的计算相同)。

表1。

模拟数据评估

程序韵律学100个基点
200个基点
500个基点
1000个基点
10 000基点
2%5%10%2%5%10%2%5%10%2%5%10%2%5%10%
BLAT(爆炸)CPU秒6855775598195384861078699512131586259926281742710
问题20%68.725.539252.97.897.186.321.497.796.43998.49994
误差%0.992.485.470.551.724.550.171.124.410.010.523.98001.28
BWA-SW公司CPU秒16512584222168118249172152234168150158134120
问题20%85.162.219.893.888.749.796.195.585.196.996.59598.498.598.1
误差%0.010.050.1700.020.13000.04000000
SSAHA2型CPU秒487279629345193222365252331182136863155415833113
问题20%85.583.878.293.493.191.996.696.596.197.797.697.4
误差%00.010.190.0100.0100.010.04000
程序韵律学100个基点
200个基点
500个基点
1000个基点
10000个基点
2%5%10%2%5%10%2%5%10%2%5%10%2%5%10%
BLAT(爆炸)CPU秒6855775598195384861078699512131586259926281742710
问题20%68.725.539252.97.897.186.321.497.796.43998.49994
误差%0.992.485.470.551.724.550.171.124.410.010.523.98001.28
BWA-SW公司CPU秒16512584222168118249172152234168150158134120
问题20%85.162.219.893.888.749.796.195.585.196.996.59598.498.598.1
误差%0.010.050.1700.020.13000.04000000
SSAHA2型CPU秒487279629345193222365252331182136863155415833113
问题20%85.583.878.293.493.191.996.696.596.197.797.697.4
误差%00.010.190.0100.0100.010.04000

从人类基因组中模拟了大约10 000 000 bp不同读取长度和错误率的数据。百分之二十的误差是indel误差,indel长度由几何分布得出(密度:0.7·0.3−1). 这些模拟读数分别与BLAT(选项-fastMap)、BWA-SW和SSAHA2(选项-454表示100和200 bp读数)的人类基因组对齐。然后将对齐的坐标与模拟坐标进行比较,以找出对齐误差。在该表中的每个单元中,这三个数字是Intel E5420 2.5 GHz CPU单核上的CPU秒数、映射质量大于或等于20(Q20)的对齐百分比以及Q20对齐中的错误对齐百分比。SSAHA2和BWA-SW报告制图质量;BLAT映射质量被估计为最佳和第二最佳对准分数之差除以最佳对准分数的250倍(基本上与BWA-SW的计算相同)。

发件人表1我们可以看到,BWA-SW显然是最快的,在所有输入上比BLAT和SSAHA2快几倍,其速度对读取长度或错误率不敏感。当查询较长或错误率较低时,BWA-SW的准确性与SSAHA2相当。考虑到短且容易出错的读取,SSAHA2更准确,尽管它必须花费更多时间来对齐这些读取。SSAHA2未在10 kb读取上进行测试,因为它最初不是为该任务设计的,因此性能不佳。使用-fastMap选项的BLAT比SSAHA2更快,但精确度更低。在默认选项下,BLAT比SSAHA2慢几到几十倍。与-fastMap模式相比,精度更高,但总体上仍低于BWA-SW模式(未显示数据)。

在内存上,BWA-SW和BLAT都使用~4 GB内存。SSAHA2使用默认选项对≥500 bp的读取使用2.4 GB,使用−454选项对较短的读取使用5.3 GB,这会增加哈希表中存储的种子序列数,从而增加内存。此外,BWA-SW支持多线程,因此如果在多核计算机上运行,每个CPU核占用的内存可能更少。SSAHA2和BLAT目前不支持多线程。

3.2.2嵌合体检测

我们首先研究了给定嵌合体读取的每个对准器的行为。为此,我们制作了两个嵌合读码,它们都由一个染色体位置的1000 bp片段和另一个位置的300 bp片段组成。两次读取之间的主要区别在于,第二次读取中的1000 bp片段具有~750 bp的重复序列,而第一次读取是高度独特的。当我们将两个嵌合读数与人类基因组比对时,BWA-SW报告了四个比对,每个片段一个,这是所需的。最新的SSAHA2未能在两次读取中找到300 bp片段的对齐,但如果我们将300 bp片段作为单个读取进行对齐,则它能够找到对齐。较旧的版本(1.0.9)默认情况下可以在第一次读取时对齐300 bp的片段,但在第二次读取时,我们需要切换到更彻底但速度慢得多的配置,将所有点击数报告给750 bp的重复。使用-fastMap的BLAT找不到第二次读取时300 bp片段的对齐。在这两个例子中,只有BWA-SW有足够的功率检测嵌合体。

此外,BWA-SW很少产生高质量的假嵌合比对。例如,假设10 000个1kb的读取有10%的错误,但在模拟中没有嵌合体,BWA-SW预测18个嵌合体读取。这些读数上错误对齐的片段的映射质量只有2.4(最多11个),这意味着BWA-SW知道这些嵌合体是不可靠的。正如预期的那样,BWA-SW在基址错误较低的情况下产生的虚假嵌合体读取更少。

3.3实际数据评估

由于缺乏基本事实,对真实数据的评估很复杂。然而,仍然可以通过比较两个校准器的结果来评估相对准确度,使用的原则是,真正的校准器往往具有相当高的校准分数,因为大多数错误都是由于找不到种子而产生的。

假设我们使用两个对齐器对齐读数A类B类并得到不同的结果。如果两者都有A类B类如果映射质量较低,则对齐是模糊的,无论哪一个对齐是否错误都无关紧要。如果A类提供了高映射质量和A类比对分数低于B类,A类对齐可能是错误的;即使A类校准分数只比B类,A类校准不可靠,并且通过以下方式获得了高绘图质量A类仍有疑问。实际上,定义“稍微好一点”的对齐分数需要对分数差异设置任意阈值,因此这种评估方法是近似的。

表2给出了454个读数的摘要,这些读数仅由一个对准器映射或映射到不同的位置,并且由BWA-SW或SSAHA2分配大于或等于20的映射质量。我们可以看到,BWA-SW往往会以高错误率(其中946个)错过短对准,这与对模拟数据的评估一致。SSAHA2因不同原因未对准。在1188次读取中,SSAHA2产生了明显错误的对齐。我们知道,由于指定了低映射质量,这些对齐是错误的,但无论如何都会错过真正的对齐。

表2。

BWA-SW和SSAHA2在实际数据上不一致的路线汇总

条件计数BWA-SW公司
SSAHA2型
平均长度平均差异平均地图Q平均长度平均差异平均地图Q
BWA-SW≥20;SSAHA2未映射0
BWA-SW≥20合理;SSAHA2<201188398.21.3%178.4198.313.1%3.9
BWA-SW≥20有问题401837.8%41.2280.39.4%2.4
SSAHA2≥20;BWA-SW未映射94675.46.3%51.2
SSAHA2≥20可信;BWA-SW<20471299.3%2.5200.58.8%34.4
SSAHA2≥20可疑185400.21.7%13.4399.22.9%216.4
条件计数体重-体重
SSAHA2型
平均长度平均差异平均地图Q平均长度平均差异平均地图Q
BWA-SW≥20;SSAHA2未映射0
BWA-SW≥20合理;SSAHA2<201188398.21.3%178.4198.313.1%3.9
BWA-SW≥20可疑401837.8%41.2280.39.4%2.4
SSAHA2≥20;BWA-SW未映射94675.46.3%51.2
SSAHA2≥20可信;BWA-SW<20471299.3%2.5200.58.8%34.4
SSAHA2≥20可疑185400.21.7%13.4399.22.9%216.4

共有137 670 454个读数统一从中选择SRR003161型分别用BWA-SW和SSAHA2对人类基因组进行了定位。如果BWA-SW和SSAHA2校准的最左侧坐标相差超过平均读取长度355 bp,则表示读取不一致。每次对齐都会计算一个分数,该分数等于匹配数减去三乘以对齐区域中的差异数(不匹配和间隙)。BWA-SW线路称为似是而非的如果BWA-SW比对得出的分数减去相同读数的SSAHA2比对得出的得分大于或等于20(即BWA-SV比对足够好);否则,BWA-SW路线称为有问题的。以类似的方式定义了合理和可疑的SSAHA2定线。在表中,“BWA-SW≥20”表示测绘质量高于20的BWA-SV定线。总之,BWA-SW遗漏了993条(=946+47条)路线,SSAHA2对齐良好,而SSAHA2遗漏了1188条;40条BWA-SW Q20路线和185条SSAHA2 Q20路线可能是错误的。

表2。

BWA-SW和SSAHA2在实际数据上不一致的路线汇总

条件计数BWA-SW公司
SSAHA2型
平均长度平均差异平均地图Q平均长度平均差异平均地图Q
BWA-SW≥20;SSAHA2未映射0
BWA-SW≥20合理;SSAHA2<201188398.21.3%178.4198.313.1%3.9
BWA-SW≥20有问题401837.8%41.2280.39.4%2.4
SSAHA2≥20;BWA-SW未映射94675.46.3%51.2
SSAHA2≥20可信;BWA-SW<20471299.3%2.5200.58.8%34.4
SSAHA2≥20可疑185400.21.7%13.4399.22.9%216.4
条件计数BWA-SW公司
SSAHA2型
平均长度平均差异平均地图Q平均长度平均差异平均地图Q
BWA-SW≥20;SSAHA2未映射0
BWA-SW≥20合理;SSAHA2<201188398.21.3%178.4198.313.1%3.9
BWA-SW≥20有问题401837.8%41.2280.39.4%2.4
SSAHA2≥20;BWA-SW未映射94675.46.3%51.2
SSAHA2≥20可信;BWA-SW<20471299.3%2.5200.58.8%34.4
SSAHA2≥20可疑185400.21.7%13.4399.22.9%216.4

共有137 670 454个读数统一从中选择SRR003161型分别用BWA-SW和SSAHA2对人类基因组进行了定位。如果BWA-SW和SSAHA2校准的最左侧坐标相差超过平均读取长度355 bp,则表示读取不一致。为每个对齐计算一个分数,该分数等于匹配的数量减去三乘以对齐区域中的差异(不匹配和间隙)的数量。BWA-SW线路称为貌似有理的如果BWA-SW比对得出的分数减去相同读数的SSAHA2比对得出的得分大于或等于20(即BWA-SV比对足够好);否则,BWA-SW路线称为有问题的。以类似的方式定义了合理和可疑的SSAHA2定线。在表中,“BWA-SW≥20”表示测绘质量高于20的BWA-SV定线。总之,BWA-SW遗漏了993条(=946+47条)路线,SSAHA2对齐良好,而SSAHA2遗漏了1188条;40条BWA-SW Q20路线和185条SSAHA2 Q20路线可能是错误的。

对于两个对准器,大多数错误对准都是由于忽略了与最佳报告对准分数相似的对准而导致的。例如,SSAHA2对齐读取SRR003161.1261578号文件与X染色体的比对质量为244,BWA-SW以相同的比对长度和编辑距离将其与2号染色体进行比对。存在两个最佳评分比对意味着读数不能唯一放置,高达244的绘图质量不准确。SSAHA2提供这种高映射质量可能是因为它忽略了2号染色体上的匹配。在这个特定的示例中,BWA-SW正确地给出了映射质量零,尽管它可能会忽略其他示例中的可选匹配。

在模拟的100和200 bp读数中,带有−454选项的SSAHA2比BWA-SW提供更好的比对。在这个真实的数据集上,BWA-SW可能更准确,因为平均读取长度相对较长(355 bp)。为了证实这一推测,我们比较了这两个校准器在99958次运行读数上的差异SRR002644号文件平均读取长度为206 bp。这一次BWA-SW错过了1092条SSAHA2 Q20路线,并产生了39条有问题的路线;SSAHA2未命中325个,产生了10个可疑的。SSAHA2在这个较短的数据集上更准确,尽管它比BWA-SW慢9倍,并且使用的内存多40%。

4讨论

BWA-SW是一种有效的算法,用于将数百个或更多碱基对的查询序列与长参考基因组对齐。对于较长的查询或错误率较低的查询,其敏感性和特异性往往更高,在这种查询序列上,BWA-SW的准确性可与迄今为止最准确的对准器相媲美。此外,BWA-SW能够检测嵌合体,嵌合体可能是由结构变化或参考装配错误引起的,这可能会对BLAT和SSAHA2构成挑战。

BWA-SW、BLAT和SSAHA2都遵循种子和扩展范式。主要区别在于播种策略。BLAT和SSAHA2将短的精确匹配确定为种子,通常长度为11或12 bp。对于k个-两个长度序列之间的mer播种L(左)分别为,预期种子数为L(左)·/4k个,或约为105与人类基因组进行比对。使用Smith–Waterman算法扩展这些种子是很昂贵的。为了减少不必要的种子扩展,BLAT和SSAHA2在默认情况下都使用非重叠种子,并且需要多个种子匹配,这对于随机序列应该很有效,但在高度重复的区域中仍然涉及许多种子扩展。BWA-SW通过在独特的地区使用一些长间隙种子来解决这个问题。在真实的生物数据上,它节省了许多不必要的种子扩展,并带来了更好的整体性能。然而,为了减少识别长种子的时间,BWA-SW只保留动态编程矩阵的一小部分,这可能会错过真正匹配的所有种子。这种启发式方法是对齐错误的主要来源,特别是对于要对齐的序列之间只有几个有效的唯一种子的简短查询。幸运的是,在长路线上,丢失所有种子的可能性很小。我们已经证明,BWA-SW与SSAHA2同样有效。

BWA-SW与BWT-SW在几个方面不同。首先,BWT-SW保证找到所有本地匹配,而BWA-SW是一种启发式算法,可能会错过真正的点击,但速度要快得多。其次,BWA-SW对齐两个FM-index,而BWT-SW对齐一个序列和一个FM-index。为查询序列构建前缀DAWG可能有助于避免在查询中重复对齐相同的子字符串,从而提高理论时间复杂性。第三,BWA-SW在内环中遍历引用前缀trie,而BWT-SW在内部循环中遍历查询序列。如果没有启发式,BWA-SW方法将损害性能,因为在遍历引用前缀trie时,我们必须以速度换取内存,而在外部循环中遍历它会更有效。尽管如此,应用Z轴-best策略要求在不完成动态编程的情况下知道匹配查询子串的得分最高的引用节点,因此只有在内部循环中遍历引用时才有效。第四,BWA-SW只报告在查询序列上基本上不重叠的对齐,而BWT-SW与BLAST一样,报告所有具有统计意义的对齐。BWA-SW保留了路线的关键信息,并生成更小、更方便的输出。对于BWT-SW,最终用户通常需要对结果进行后处理,以筛选出他们不感兴趣的许多对齐。总之,鉴于大规模实际数据,BWA-SW已调整为实用性。

BWA-SW的高速很大程度上来自两种策略:使用FM索引和抑制更好匹配中包含的短重复匹配。虽然第一种策略不适用于基于散列表的算法,如SSAHA2和BLAT,但第二种策略可以在此类程序中实现,并且可以通过节省大量重复比对的构建时间来大大加快这些程序的速度。尽管BWT的使用减少了重复中不必要的对齐,但与哈希表查找相比,每个BWT操作都有一个很大的常量。如果基于散列表的算法结合了这些特性,那么它们仍然可能比BWA-SW更快。

致谢

我们感谢Zemin Ning对SSAHA2算法的有益评论,也感谢Kimmo Palin提供有关DAWG的文献。

基金:惠康信托077192/Z/05/Z。

利益冲突:未声明。

参考文献

阿尔特舒尔
旧金山
缺口BLAST和PSI-BLAST:新一代蛋白质数据库搜索程序
核酸研究。
1997
,卷。 
25
(第
3389
-
3402
)
模糊
A类
识别文本子单词的最小自动机
西奥。计算。科学。
1985
,卷。 
40
(第
31
-
55
)
Burrows公司
M(M)
惠勒
流行音乐播音员
一种分块无损数据压缩算法
技术报告124
1994
加利福尼亚州帕洛阿尔托
数字设备公司
开斋节
J型
单聚合酶分子实时DNA测序
科学
2009
,卷。 
323
(第
133
-
138
)
费拉吉纳
P(P)
曼齐尼
应用程序的机会主义数据结构
第41届计算机科学基础研讨会论文集(FOCS 2000)
2000
美国加利福尼亚州雷东多海滩
(第
390
-
398
)
小时
Wong(王)
世界卫生组织
SeqMap:将大量寡核苷酸映射到基因组
生物信息学
2008
,卷。 
24
(第
2395
-
2396
)
肯特
WJ公司
BLAT–类似BLAST的对齐工具
基因组研究。
2002
,卷。 
12
(第
656
-
664
)
台湾
DNA的压缩索引和局部比对
生物信息学
2008
,卷。 
24
(第
791
-
797
)
朗米德
B类
超快速和记忆有效的短DNA序列与人类基因组的比对
基因组生物学。
2009
,卷。 
10
第页。 
25兰特
 
小时
杜宾
R(右)
使用Burrows-Wheeler变换快速准确地进行短读对齐
生物信息学
2009
,卷。 
25
(第
1754
-
1760
)
小时
,等人
使用映射质量分数映射短DNA测序读取和调用变体
基因组研究。
2008
,卷。 
18
(第
1851
-
1858
)
小时
序列对齐/映射格式和SAMtools
生物信息学
2009
,卷。 
25
(第
2078
-
2079
)
R(右)
SOAP:短寡核苷酸比对程序
生物信息学
2008
,卷。 
24
(第
713
-
714
)
小时
缩放!绘制了无数寡头
生物信息学
2008
,卷。 
24
(第
2431
-
2437
)
妈妈
B类
PatternHunter:更快更敏感的同源搜索
生物信息学
2002
,卷。 
18
(第
440
-
445
)
米克
C类
,等人
OASIS:一种在线准确的生物序列局部对齐搜索技术
第29届超大数据库国际会议记录(VLDB 2003)
2003
德国柏林
(第
910
-
921
)
莫古利斯
A类
用于生产超级用户搜索的数据库索引
生物信息学
2008
,卷。 
24
(第
1757
-
1764
)
Z轴
SSAHA:一种大型DNA数据库的快速搜索方法
基因组研究。
2001
,卷。 
11
(第
1725
-
1729
)
皮尔逊
WR(额定功率)
利普曼
流行音乐播音员
改进的生物序列比较工具
程序。美国国家科学院。科学。美国
1988
,卷。 
85
(第
2444
-
2448
)
隆隆声
性虐待
SHRiMP:短颜色空间读取的精确映射
公共科学图书馆计算。生物。
2009
,卷。 
5
第页。 
e1000386
 
史密斯
AD公司
使用质量分数和更长的读取时间可以提高Solexa读取映射的准确性
BMC生物信息学
2008
,体积。 
9
第页。 
128
 
杂草
D类
RazerS–具有灵敏度控制的快速读取映射
基因组研究。
2009
,卷。 
19
(第
1646
-
1654
)
Z轴
一种贪婪的DNA序列比对算法
J.计算。生物。
2000
,卷。 
7
(第
203
-
214
)

作者注释

副主编:Dmitrij Frishman

这是一篇根据知识共享署名非商业许可条款发布的开放存取文章(http://creativecommons.org/licenses/by-nc/2.5)它允许在任何媒体上无限制地进行非商业性使用、分发和复制,前提是正确引用了原始作品。