1简介
由于DNA测序的显著改进,需要评估大量片段序列(如短读)之间的序列相似性。我们解决了按编辑距离枚举大型字符串池中所有邻居对的问题,其中插入、删除和替换的成本是一。也就是说,给定一组n个等长序列ℓ,秒1,…,秒n个,任务是找到编辑距离最多的所有对天,
它通常被称为全对相似性搜索.
所有配对搜索都出现在重要的生物任务中。例如,在序列聚类所需的所有对对齐中查找种子匹配是必需的(Abouelhoda公司等。2004年). 这样的对齐可以用于检测和纠正短读中的错误(曲等。,2009年). 在的第一步从头开始基因组组装(辛普森等。, 2009;Zerbino和Birney,2008年),短读取被分解为k个-mers和后缀–长度的前缀匹配k个检测到−1。在大多数情况下,由于时间限制,会使用精确匹配。通过使用近似匹配,可以延长连续梁的长度,从而使最终装配质量更好。这个问题通过收集所有对的相似性来减少到所有对的相似性搜索k个-1个前缀和后缀放入序列池。从输出中,只报告前缀-后缀对。
基本上,大多数流行的方法都是通过以下两种方法之一或它们的组合来解决搜索问题。(i) 寻找共同点k个-mer并验证匹配(Lipman和Pearson,1985年;辛普森等。, 2009;沃伦等。, 2007;杂草等。, 2009;Zerbino和Birney,2008年). (ii)索引结构中的回溯(即后缀数组和FM-index)(朗米德等。, 2009;Li和Durbin,2009年;锂等。, 2009;拉贾塞卡兰等。, 2005;Sagot,1998年;Trapnell公司等。, 2009). 第一种类型很常见k个-字符串中的mers(即种子匹配),并验证两个字符串是否共享k个-mer确实是邻居,通过动态编程扩展匹配。当绳子足够长时,它工作得很好。然而,当字符串较短且阈值天较大,共享长度k个-mers的差距如此之小,以至于需要验证太多的候选对。第二种类型将字符串存储到索引结构中,最常见的是后缀数组。然后,通过遍历相应后缀树的节点找到相似的字符串。如果天很小,例如。天≤2,并用于最先进的短读数绘图工具,如BWA(Li和Durbin,2009年)、领结(朗米德等。,2009年)和SOAP2(锂等。, 2009). 然而,由于天变得更大,主要是因为复杂性是指数级的天并且还不知道有效的修剪。ELAND和SeqMap(江和黄,2008)将序列分解为块并使用多个索引存储所有k个-块的串联。显然,与BWA相比,它需要更多的内存,这在许多情况下都会有问题。多重排序(联合国,2008年)使用多个子串匹配有效地缩小搜索范围,但它只能根据汉明距离查找邻居。
我们的方法称为幻灯片排序通过高效的模式增长算法找到公共子串链,该算法已成功应用于项集挖掘等数据挖掘任务(汉族等。2004年). 模式对应于子字符串序列。所有模式的空间都组织为一棵树,并系统地遍历。我们的方法不依赖任何索引结构来避免存储开销。相反,在模式增长期间,使用基数排序来查找等效字符串。为了证明我们算法的正确性,首先证明了任何邻居对中都存在一个公共子串链。此外,我们故意避免通过重复检查多次报告同一对。因此,我们的方法很容易扩展到1000万个序列,并且对于短序列和大半径的序列,比种子匹配方法和后缀数组要快得多。
本文的其余部分组织如下。第2节介绍了我们的算法。在第3节给出了计算实验结果。第4节文章的结论。
2方法
两个相似的字符串共享公共子字符串串联因此,我们可以通过系统地检测普通字符串的链来检测类似字符串。在继续进行算法之前,让我们先描述一下基本属性。将间隔1除以,…,ℓ 进入之内b条任意长度的块w个1,…,w个b条, ∑我=1b条 w个我= ℓ. 每个块的起始位置定义为q个我= 1 + ∑j个=1我−1 w个j个字母表表示为∑。我们假设数据库中的每个字符串{秒我}我=1n个包括ℓ 信件,秒我∈|∑|ℓ.给定两个字符串秒,吨,秒=吨如果所有字母都相同,则保持不变。位置的子字符串我到j个被描述为秒[我,j个].
A类图案长度的k个定义为字符串和块索引的序列,
哪里x个我∈ |Σ|w个年我, 1 ≤年1<年2,…, <年k个≤b条.图案X与字符串匹配秒我带偏移量{第页= (第页1,…,第页k个),如果
的所有事件X在数据库中表示为
为了方便起见,一个索引集我(X)定义为出现在C类(X). 中的序列数我(X)定义为|我(X)|.
相邻对和模式之间的关系由以下定理来表征。
T型神灵1. —
如果是我 和j个 长度相等ℓ是邻居,即EditDist(秒我,秒j个) ≤天,我<j个,存在长度为b的图案X−d,使X与s匹配我 具有零偏移(第页1=第页2= … =第页b条−天= 0)和匹配sj个 有界偏移−⌊天/2⌋ ≤第页k个≤ ⌊天/2⌋,k个= 1,…,b条−天.
P(P)屋顶. —
有多种可能的对齐方式秒我和秒j个。对齐的特征是匹配的数量米,不匹配(f)差距克我在里面秒我和缺口克j个在里面秒j个.长度秒我等于米+(f)+克j个和的秒j个是米+(f)+克我,因为里面有任何字母秒我与中的任一字母对齐秒j个或中的间隙符号秒j个反之亦然。因此,我们得到克我=克j个≤ ⌊天/2⌋考虑到最大间隙数不超过天因此,任何字母的对齐位置都在⌊的范围内天/2⌋字母从其原始位置。
让我们分开秒我进入之内b条长度块w个1,…,w个b条由于不匹配的数量最多天,至少b条−天块与中的对应块完全匹配秒j个在任何路线上。此外,由于对齐位置绑定在⌊内天/2⌋从它们的原始位置,任何块的匹配对应秒我可以在中找到秒j个在-⌊之间的偏移范围内天/2⌋和⌊天/2⌋. ▪
演示了使用b条=5,天= 3. 这个定理意味着任何邻居对都有一个b条−天公共块和相应的块彼此靠近。它是我们后面介绍的算法的基础。
块大小5和编辑距离阈值3的示例模式。秒我与匹配X在第一个块和第三个块中没有偏移。秒j个与匹配X第一个块中没有偏移,但第三个块中有-1偏移。
2.1模式增长
在我们的算法中|我(X)|≥2由递归模式增长算法枚举。在模式增长算法中,构造了一个模式树,其中每个节点对应一个模式(). 深度节点k个包含长度模式k个.
该方法的模式生长和剪枝过程。模式是通过以深度优先的方式遍历树来枚举的。在每个节点中,通过对序列池中的子字符串进行排序(“ATA”、“TAT”、“TTA”年1= 1). 删除无用的图案(本例中为“TA”)。添加剩余元素以生成新图案。此过程通过递归调用执行,直到模式大小达到b条−天.
首先,生成长度为1的图案,如下所示。对于每个区块年1= 1,…,天+1,字符串池是通过收集{秒我}我=1n个开始于q个年1− ⌊天/2⌋,…,q个年1+⌊天/2⌋。将基数排序应用于字符串池并扫描排序结果,可以检测到等效字符串的重复(). 每个长度为1的图案,表示为X1,构造为重复字符串的组合x个1和年1,
同时,所有事件C类(X1)记录。如果秒我匹配相同的模式X1通过几个不同的偏移量,只记录最小的偏移量。它们形成与模式树深度1对应的节点。
给定一个模式X吨长度的吨,它在模式树中的子代的生成方式类似于如下所示。对于每个年吨+1=年吨+ 1,…,天+吨+1,字符串池是通过收集我(X吨)开始于q个年吨+1− ⌊天/2⌋,…,q个年吨+1+ ⌊天/2⌋. 因为字符串池仅由发生集构成,所以池的大小随着模式的增长而急剧减小。通过排序和扫描,下一个字符串x个吨+1识别,并将模式扩展为
和事件C类(X吨)更新为C类(X吨+1)还有。
为了避免生成无用的模式,一旦支持降到2以下,就会修剪模式树。此外,如果树中没有字符串,则会对其进行修剪我(X)匹配的X零偏移。随着模式以深度优先的方式增长,内存会随着模式的扩展而保留,然后在模式收缩以访问另一个分支时立即释放。这种动态内存管理使峰值内存使用量保持得相当小。
2.2从图案到成对
如定理1所示,每个相邻对()出现在索引集中我(X)具有至少一种图案。由于其中一对必须具有零偏移,因此合格对集描述为
因为不是的所有成员P(P)X对应于邻居,我们必须通过实际的编辑距离计算来验证它们是否是邻居。
这里的一个问题是同一对(我,j个)可能出现在许多不同模式的索引集中。也可能是这对(我,j个)在同一个索引集中,是从不同的偏移量导出的。在大多数应用中,最好确保没有配对报告两次。该问题的直接解决方案是通过存储所有对来检查以前是否报告了新的对,这需要大量内存。我们提出了一个替代解决方案非规范的不使用任何额外内存的配对如下。
一场比赛秒我和秒j个可以以各种方式发生,每种方式都可以描述为元组(年,第页),其中年=年1,…,年b条−天描述模式中的块,并第页是模式匹配的偏移量秒j个.我们将规范匹配定义为字典序最小的匹配年和第页,其中优先考虑年例如,考虑一下这种情况秒我=AATT公司,秒j个=ATAT公司,天=2,所有块宽度设置为1。有10种不同(年,第页)对,如所示,其中匹配的残留物以红色方块显示。在这种情况下,(1)是规范的。其中,具有重叠正方形的匹配没有正确对齐。我们不排除此类对以避免额外运行动态编程。
全部(年,第页)第页,共页秒我=AATT公司,秒j个=阿塔特。匹配残留物显示在红色方块中。由于(6)和(10)中的红色方块重叠,因此它们不符合正确的对齐方式。
判断给定匹配项(年,第页)无论是否规范,只要检查是否存在另一个字典序较小的匹配就足够了。更确切地说,由年,第页不规范,如果存在块1≤z(z)≤最大值(年),z(z)∉年和偏移-⌊天/2⌋ ≤第页≤ ⌊天/2⌋这样
此规范性检查可以在中完成O(运行)(天ℓ) 时间。
整个算法的伪代码如算法1所示。在第18行中,计算宽度为2的对角线条纹就足够了天+DP矩阵的1。因此,距离计算在O(运行)(天ℓ) 时间。
2.3备注
只需稍作修改,SlideSort就可以处理间隙打开和间隙扩展惩罚。将缺口开放和扩展成本定义为γo(o)和γe(电子)分别是。将不匹配、间隙打开和间隙扩展的数量表示为(f),克o(o)和克e(电子)分别是。然后,我们的所有对相似性搜索问题被重新定义为查找对,如下所示(f)+克o(o)γo(o)+克e(电子)γe(电子)≤天.将每个序列中的间隙数表示为克我和克j个.然后,克e(电子)=克我+克j个和克o(o)≥2(如果克e(电子)≠ 0),克o(o)=0(如果克e(电子)=0)。什么时候?克e(电子)≠0,我们有(克我+克j个) γe(电子)≤天− 2 γo(o)由于两个序列的长度相等,所以间隙的数量也相等,克我=克j个导致以下不平等,
因此,偏移量第页k个,用于k个=1,…,b条−天,以为界
什么时候?克e(电子)=0,我们可以通过零偏移找到所有对,因此偏移范围(三)也涵盖了这个案例。请注意,块大小b条必须大于max(天, ⌊天/(γo(o)+ γe(电子))⌋). 当γo(o)和γe(电子)大于替代成本。
值得注意的是,SlideSort可以处理长度稍有不同的序列,而无需进行任何必要的修改。详见补充文件1中的补充方法。
2.4复杂性
用σ=|∑|字母大小表示。SlideSort的空间复杂度为O(运行)((b条−天)dn(数字网络)日志n个+n个ℓlogσ),因为它需要一个指针数组来描述模式树,并且必须保留原始字符串。表示方式米索引集中包含的所有对的数量我(X). SlideSort的时间复杂性为O(运行)(b条天−1 天b条−天ℓn个+医学博士ℓ), 其中第一部分用于排序,第二部分用于编辑距离计算。时间复杂性取决于修剪的有效性米后一部分的最坏情况是O(运行)(n个2 天ℓ) 当所有输入短读都相同时。然而,在大多数情况下,短文阅读是多种多样的米预计规模将比O(运行)(n个2).
通过在长度为n个ℓ. 增强的后缀数组可以通过O(运行)(ℓ天+1σ天n个)时间和O(运行)(n个日志n个+n个ℓ对数σ)空间(Abouelhoda公司等。2004年). 这种时间复杂性本质上是通过在距离内生成所有变量来实现的天所有序列并找到相同的对。区别在于后缀数组的时间复杂度取决于字母表的大小,而不是SlideSort的大小。因此,SlideSort也可以应用于大型字母表(即蛋白质)。
3个实验
从NCBI序列读取存档(网址:http://www.ncbi.nlm.nih.gov/sra/),我们获得了两个数据集:人类HapMap配对测序(ERR001081)和IMR90细胞系全基因组鸟枪亚硫酸氢盐测序(SRR020262)。它们将分别称为数据集1和2。数据集1的序列长度为51,数据集2的序列长度是87。这两个数据集均由Illumina Genome Analyzer II生成。从原始fastq文件的顶部选择不包含“N”的读取。我们的算法用C++实现,用g++编译。所有实验都是在带有Intel Xeon X5570(2.93 GHz)和32 GB RAM的Linux PC上进行的。所有实验都只使用一个内核。
作为竞争对手,BWA(Li和Durbin,2009年)和SeqMap(江和黄,2008)在众多备选方案中进行选择,因为这两种方法代表两种完全不同的方法,回溯和区块组合.BWA是使用索引回溯的最佳方法之一,而SeqMap采用类似ELAND的方法,对所有块组合使用多个索引。SlideSort还与计算所有对的编辑距离的朴素方法进行了比较。BWA和SeqMap通过从所有短读创建索引并使用相同的读集进行查询,应用于所有对的相似性搜索。
请注意,BWA和SeqMap最初不是为所有对的相似性搜索而设计的,而是为读取映射而设计的。读取映射需要更大的搜索空间。尽管在不同目的的工具之间进行公平比较是困难的,但我们将映射工具作为竞争对手使用,因为据我们所知,没有任何工具可广泛用于所有对的相似性搜索。
对于我们的方法,数字b条必须确定块的数量。在以下实验中,我们设置b条相对于距离阈值天作为b条=天+k个.给,k个对应于图案大小。在以下实验中,我们尝试了k个=1,…,5,并报告了最佳结果。
两个数据集的邻居对数量如补充图S1所示。我们确认SlideSort和naive方法报告的邻居对数量完全相同,这确保了SlideSort实现的正确性。
3.1计算时间和内存使用
根据距离阈值绘制计算时间天.SlideSort在所有配置中都始终更快。随着序列数量的增加和距离阈值的增加,与BWA和SeqMap的差异变得越来越明显。由于30 GB的内存限制和30万s的时间限制,并没有获得所有结果。比较SlideSort、SeqMap和BWA的峰值内存使用情况。我们分别测量了BWA的索引步骤和搜索步骤的内存使用量,因为BWA被设计为分别执行这些步骤。在大多数配置中,BWA用于搜索步骤的峰值内存是最小的,而SlideSort的峰值内存与BWA的峰值索引内存相当或略好。10万次短读的详细结果如所示.
两个短读取数据集的计算时间。在这四种方法中,“naive”表示详尽的距离计算。
两个短读取数据集的内存使用情况。BWA的内存使用分别针对索引步骤(索引)和搜索步骤(搜索)进行评估。
表1。
| 数据集1
| 数据集2
|
---|
| 幻灯片排序 | BWA公司
| SeqMap(序列图) | 天真 | 幻灯片排序 | BWA公司
| SeqMap(序列图) | 天真 |
---|
| | 索引 | 搜索 | | | | 索引 | 搜索 | | |
---|
d=1 | 0.2 | 2.34 | 3.25 | 7.68 | 8743.46 | 0.33 | 5.07 | 6.91 | 39.59 | 15 678 |
d=3 | 0.85 | 2.37 | 562.63 | 205.26 | 23 796.1 | 1.84 | 5.09 | 1647.16 | 10 698.6 | 39 046.3 |
d=5 | 6.56 | 2.19 | 19 697.67 | 93 115.2 | 38 179.5 | 5.56 | 5.08 | 12 876.88美元 | > 300 000 | 65 244.6 |
BWA在空间复杂性方面效率最高,因为它的索引大小不依赖于距离阈值。相反,由于回溯中遍历的节点数量激增,随着编辑距离阈值的增加,BWA的时间复杂性迅速恶化。相比之下,SeqMap对所有关键块的组合进行索引和散列,这导致了巨大的内存使用量。SlideSort与SeqMap类似,它考虑了所有块组合,但内存效率更高。区别在于SlideSort是一种无索引方法,它通过深度优先遍历动态生成模式树。它允许我们在内存中只维护树的必要部分。
3.2图案尺寸的影响
研究图案尺寸的影响k个关于效率。除了天=1,最佳设置在附近k个= 2–4. 我们的方法k个=1大致对应于单种子方法,因此该结果表明使用链的有效性。总的来说,计算时间对选择k个.
不同方法的性能比较k个通过100万次短读进行评估。
3.3与单个种子的比较
我们的算法使用一系列公共子串来缩小搜索范围。与使用k个-mer派生候选对时,子字符串的总长度可以比k个-没有失去任何邻居。它产生了更高的特异性,从而减少了候选对的数量。可以检测多个而不是链k个-mers并验证那些包含多个匹配项的对(Burkhardt和Kärkkäinen,2002年). 然而,与相同总长度的链相比,这种方法的特异性更低,因为根据定理1,链中每个元素的匹配位置都是严格局部化的。
比较我们的方法和单个种子生成的候选对的数量(k个-情节中的mer)。它对应于编辑距离计算的数量。单种子法有两种变体:k个-mer/nonredundant“将以前报告的配对存储在内存中,而不包括候选中先前报告的配对。”k个-mer/redundant不使用额外内存,但对同一对进行多次计数。这里我们设置k个-mer到Ş/天这样就不会失去邻居。在图中,我们可以看到算法中候选对的数量显著减少。请注意,候选对的数量以对数刻度显示。在我们的方法中,在最大的模式大小下,候选的数量是最小的,因为子串的总长度是最大的,特异性是最优的。然而,由于模式的搜索空间被扩展,因此在这种情况下总计算时间不是最优的。
候选对数量的比较。对10万篇短文进行了评估。使用k个=1,…,5。'“邻居对”表示数据中邻居对的实际数量k个-mer/非冗余“和”k个-mer/redundant表示单种子方法的两种变体(见正文)。
3.4短文的聚类分析
SlideSort的一个主要应用是分层序列聚类,例如,它将用于纠正短读中的错误和元基因组映射的预处理。SlideSort提供无向图G公司,其中顶点表示短读取,而加权边表示相邻对的编辑距离。在层次聚类算法中,单链路是最具可扩展性的(曼宁等。, 2008). 由于单链聚类的树状图与最小生成树同构(Gower和Ross,1969年),可以通过Kruskal或Prim算法构造最小生成树(MST)来执行单链接聚类(克鲁斯卡尔,1956年;Prim,1957年).
由于存储所有边可能需要大量内存,因此我们使用了一种著名的在线算法来构建MST(Tarjan,1983年). 它从一系列边创建MST,沿途丢弃不必要的边。它基本上保持所有无循环连接的组件,如果循环由新边生成,则会从循环中删除最重的边。在我们的实验中,与SlideSort查找相似对相比,查找MST的额外计算时间微不足道().通过3D可视化工具Walrus,以编辑距离阈值3可视化数据集2的10000 000次短读中发现的最大MST(http://www.caida.org/tools/visualization/walrus/).
使用编辑距离阈值3从10000 000次短读的相邻图中可视化大型MST。左图显示了112 995个节点的360个MST,每个MST包含100多个节点。右图重点关注由6990个节点组成的最大MST。从这些MST中直接获得单链聚类的树状图。
表2。
编辑距离阈值为3的两类短读数据集搜索对和查找MST的计算时间比较
读取次数 | 数据集1
| 数据集2
|
---|
| 仅搜索对 | 查找MST | 仅搜索对 | 查找MST |
---|
10 000 | 0.08 | 0 | 0.06 | 0 |
10万 | 0.85 | 0.08 | 1.84 | 0.01 |
1 000 000 | 23.08 | 1.08 | 31.34 | 2.06 |
10 000 000 | 495.15 | 17.69 | 554.09 | 5.61 |