Bi-level error correction for Pa
最新的测序技术,如太平洋生物科学公司(PacBio)和牛津纳米孔机器可以产生数千个碱基长度的长读,这比Illumina机器产生的数百个碱基长度的长得多。
然而,这些长时间的读取容易导致更高的错误率,例如15%,使得下游的分析和应用程序非常困难。误差校正是一个提高测序数据质量的过程。最近提出了混合校正策略,将低错误率的Illumina读序列结合起来,解决长读序列中噪声较大且性能较好的排序问题。为了进一步提高PacBio长读的质量,本文提出了一种新的双色混合纠错方法——双色混合纠错框架。在第一个层次上,我们的方法使用基于de Bruijn图的纠错思想,迭代地在k-mer长度不断增加的情况下对实体k-mer进行路径搜索。在第二层中,我们结合了与第一级不同参数下的处理结果。特别地,使用多序列对齐算法对这些类似的长读取进行对齐,然后使用投票算法确定读取的每个位置的最终碱基。我们比较了三种最先进的方法在三个实际数据集上的性能。结果表明,双色始终能达到最高的单位比。与当前的方法相比,双色还可以获得更高的比对比(>:3%)和更高的比对读取次数。在第三个数据集上,我们的方法在比对读数和基因组覆盖方面与现有方法具有很强的竞争力。我们算法的c++源代码可以在https://github.com/yuansliu/Bicolor免费获得。
Index Terms—error correction, PacBio long reads, de Bruijn graph, multiple sequence alignment
Introduction
第二代测序技术是一种高通量、低成本、高质量的测序技术,先后应用于重测序、从头测序、转录组分析和宏基因组学等领域。然而,它产生的读数相对较短——Illumina产生的读数的中位数长度为100bp。短读大大降低了连续性,为处理重复序列[4]提供的信息较少,装配难度较大。新一代测序(NGS)技术,例如太平洋生物科学和牛津纳米孔平台,可以产生长达50000的长读数。长读比短读提供更多的信息来解决复杂重复的问题。在Pacific Biosciences 实时测序中,早期化学反应的总体错误率比Illumina平台的[6]大约高两个数量级,这导致长时间读取的错误率更高(至少15%)。
极高错误率的缺点给下游分析和应用带来了挑战
虽然已有许多算法用于校正短读,但这些算法并不直接适用于校正长读。这是因为长读主要由插入和删除(indels)错误控制——indels比substitution常见约15倍,而短读的主要错误类型是替换。近年来,针对长读纠错问题提出了几种算法。根据是否使用短读,这些算法可以分为两类。第一类是自校正方法,它只使用有噪声的长读,包括HGAP、Canu和LoRMA方法。自校正方法存在许多局限性,如所需的高覆盖率和较大的计算成本。因此,为了提高长读纠错的性能,开发了第二类混合纠错。
混合校正方法利用短读来校正长读中的误差。由于短读比长读的错误率低(约1%),所以短读为长读纠正提供了一个很好的模板。混合校正方法有两个主要思想。第一种方法是在短读和长读之间构建映射,然后通过映射纠正长读。例如,PacBioToCA使用映射信息来选择重叠部分,这些重叠部分在每次长读时转换为短读序列的平铺。然后,通过平铺短读序列的多重对齐,为每次长读生成一个新的一致序列。LSC在映射之前使用了均聚物压缩(HC)转换。然后,它发现四种类型的纠正点:HC点、不匹配、删除和插入。这些点被它们的短读一致序列所取代。proovread方法通过映射信息和投票策略来计算一致性。proovread的新颖之处在于迭代校正步骤,它由三个预校正和一个整理周期组成。CoLoRMap基于映射信息构建加权对齐图。然后,应用经典的最短路径算法构造编辑分数最小的修正区域。对于长读的某些区域不包括在短读中,使用单端锚(OEA)扩展校正区域
但是,这些方法分别映射短读,并且不利用短读发生的上下文。其他方法,如LoRDEC和Jabba,从短读构造一个de Bruijn图(DBG),然后使用序列对齐算法将长读与DBG对齐。
LoRDEC通过寻找最佳路径(例最小化长读的两个实体k-mers之间的编辑距离)将长读对齐到DBG。Jabba使用种子和扩展策略将长读与DBG对齐。这些方法有一个共同的局限性,即长读校正的质量在很大程度上取决于k-mer的长度。如果用户设置了一个较大的k-mer,那么只有几个DBG可以映射到长读。因此,许多错误的碱基对无法纠正。另一方面,如果用户设置一个小的k-mer,那么很多DBGs都可以映射到长read这使得选择最终结果变得很困难。
本文提出了一种提高长阅读质量的新方法——双色法。我们的方法有两个层次的处理。在第一个层次,我们对固体k-mers的选择设定了严格的条件。该选择准则克服了k-mer长度影响长读到DBG映射质量的限制。然后使用不同长度的k-mers迭代地纠正长读。因此,在k-mer初始长度不同的情况下,我们可以得到几个预校正的长读数。在第二层,我们使用多重序列对齐(MSA)算法来对齐这些类似的预校正长读[21],然后使用投票算法来获得,最后更正长读。该方法的核心思想是利用不同长度的k-聚类(k-mers),将预先校正的长读集合组合起来。实验结果表明,该方法比现有的纠错方法具有更好的纠错性能。
2.METHODS
我们的双色算法是一种双级结构的长读纠错算法。Fig1为双色的原理图。第一级由n个迭代校正器组成,每个校正器使用不同长度的k-mer。迭代校正器在初始k- mer下迭代地校正了长读m次的噪声。该迭代校正器的初始k-mer在后续迭代中增大了k的大小。因此,我们可以在第一级得到n个预先校正的长读数。然后利用MSA和第二层的投票算法对这些预校正的长读进行处理。第二级的输出是最终纠正的长读。 image image.gif
2.1****第一级:长读预校正
迭代校正是第一级计算的核心。类似的迭代方法已用于短读程序集[22]、[23]、短读校正[24]和自校正[14]。将LoRDEC[20]修改为迭代版本(称为iLoREDC)来执行计算。LoRDEC有三个主要步骤:(1)使用短读构造DBG;(2)长读时确定固体/弱k-mers;(3)在两个实体k-mers之间编辑距离最小的DBG中搜索路径。DBG是大多数第二代汇编器(如Velvet[25]、Minia[26])的核心。DBG将短读连接到图中。然后,通过找到可靠的k-mers,长读可以与DBG对齐。在这里,长读中的实k-mers被保留为正确的子字符串,假定没有错误。我们假设误差只存在于弱k-mers中。因此,弱k-mers可以通过搜索实k-mers之间的路径来纠正。
设L为噪声长读,奇数k为k-mers的初始长度,m为迭代次数。iLoRDEC的迭代校正过程如下:
步骤1:使用短读来构建DBG,如果对应的k-mers与(k - 1)基重叠,则边连接两个节点。这些短读中出现少于s次的k-mers被过滤掉了。
步骤2:在长读L中查找实k-mer。给定长读L的所有k-mer,如果L的i-th k-mer和(i+ 1)-th k-mer都在DBG中,则L的i-th k-mer为实k-mer,否则为弱k-mer。一个或多个连续的实体k-mers构成一个实体区域,一个或多个连续的弱k-mers构成一个弱区域。特别地,长读的开头和结尾的弱区分别称为头区和尾区。
步骤3:纠正L的薄弱部位。在DBG中找到L的实体区域之间的路径来修正弱区域。如果找到多个路径,则选择编辑距离最小的路径作为校正序列。
步骤4:通过搜索到这些区域的最小编辑距离的路径来纠正这些头部和尾部区域。
步骤5:使用Dijkstra算法找出第一个和最后一个实体k-mers之间的最短路径。
步骤6:更新m = m - 1和k = k + 2。
步骤7:如果m >为0,转到步骤1,使用修正后的序列作为输入。否则,输出这个修正后的序列。校正序列的输出称为预校正长读。
步骤3、4和5的详细信息可以在[20]中看到。
通过与LoRDEC的比较,iLoRDEC做了一些修改:
1)在步骤2中,我们加强了固体k- mers的选择。通过LoRDEC,如果L的i-th k-mer在DBG中,它就被视为固体。如果我们使用一个大k,长读可能不包含一个固体k-mer。因此,长读中的错误碱基对不会被纠正。如果我们设置一个小k,长读可以有很多固定的k-mers。长读过程中,固体k-mers序列的重复和假阳性现象时有发生,可能会导致长读过程的过校正。为了克服这个问题,如果DBG中存在两个后续的长读k-mer,我们只选择第一个k-mer作为实值。该选择准则可提高固体k-mers的可靠性。
2)iLoRDEC在步骤3、4和5中只执行一次传递,而LoRDEC在两个方向上执行两次传递。LoRDEC纠正了长读的反向互补,并在第一遍中输出一个经过纠正的长读。在第二遍中,LoRDEC将校正后的长读转换为其反向互补序列,并对该序列进行校正。以下两个原因促使Salmela和竞争对手[20]进行了两次传球:
(1)下一遍使用新的实体k-mers作为起始节点;
(2)不同区域的终点指向不同的路径。实际上,iLoRDEC是一种迭代算法,在后续的迭代过程中,使用新的实体k-mers作为开始节点或结束节点。因此,我们不考虑反向互补的长读。
3)加入步骤6和7,对m轮不同长度的k-mers进行迭代,每轮k增加2。
第一级有n个迭代校正器。每个校正器使用不同的k-mer初始长度迭代地校正长读。因此,我们可以在这个级别上获得n个预先校正的长读。
2.2****二级:基于MSA的校正
MSA在当前分子生物学中得到了广泛的应用,如推断序列同源性[27]、改进蛋白质二级结构预测[28]和进行系统发育分析[29]等。在我们的校正框架的第二级,MSA用于校正从第一级导出的那些预先校正的长读数。在我们的实现中使用了tool MUSCLE[30]。然后使用一个简单的投票算法生成最终的修正序列。这个简单的投票算法选择最频繁的碱基作为每个位置的最终结果。
为了说明这一点,Fig 2中描述了一个包含4个序列的示例,其中4个序列是4个预先校正的长读。我们用肌肉来校正这些预先校正过的长读数。由于S1、S2、S4的二碱基为C, S3的二碱基为A,这些预校正长读的第二位置最频繁的碱基为C。然后,最后纠正读数的第二个基数是C。
image image.gif
** 3****结果与分析**
本节给出了修正结果和一些分析。通过与现有的三种算法:LoRDEC[20]、CoLoRMap+OEA[15]和CoLoRMap [15]。在[15]、[20]、CoLoRMap和LoRDEC中均有报道,与pacBioToCA、LSC和proovread相比,取得了相当的性能。我们没有直接与pacBioToCA、LSC或proovread比较我们的性能。All 实验 进行 了 计算 集群 上 运行 Red Hat Enterprise Linux6.7(64 bit)2x2.3 GHz Intel Xeon E5-2695 v3(14Cores)和128 GB RAM.
3.1****数据集
该算法在三个数据集上进行了测试:大肠杆菌(E.coli)的一个细菌基因组、酿酒酵母(Saccharomyces cerevisiae)和果蝇(果蝇)的两个真核基因组。它们是[15]中使用的基准数据集。这些数据集的更多细节见Table 1。
image image.gif
3.2****与LoRDEC、CoLoRMap、CoLoRMap+ OEA比较
在将Bicolor算法与LoRDEC[20]、CoLoRMap[15]和CoLoRMap+OEA[15]进行性能比较时,使用了默认的参数设置(见Table 2)。
image image.gif
为了测量这些校正方法的性能,我们使用BLASR[31]将长读序列对齐到参考基因组。对于每个读取,我们在选项“-noSplitSubreads -bestn 1”下存储一个最佳对齐。然后我们计算了以下统计数据作为度量:
对齐读数:与参考基因组对齐的长读数。
比对比率:比对碱基的数目与长读碱基总数的比率。
同源比:在参考基因组中,匹配碱基的数量与对齐区域的长度之比。
基因组覆盖范围:通过长时间的读取,基因组排列区域所占的比例。
对齐读取的数量衡量校正算法的通量。与参考基因组对齐的长读的数量越多,表示更正的长读的噪音越大。比对比率和恒等比表示数量和对齐底座的质量。他们一起测量校正的准确性。基因组的覆盖范围定义了参考基因组被正确读取所覆盖的程度。这种评价方法已被目前最先进的[20]、[15]、[18]、[14]方法广泛采用。
image image.gif
比较结果如Table3所示。在大肠杆菌的数据集上,所有的方法都能达到接近性能的方法是最高的。通过双色比对回参考基因组的读数至少比其他方法多471个。与LoRDEC和CoLoRMap相比,我们的比对率分别提高了3.2%和1.7%。而LoRDEC与CoLoRMap的比对比对比甚至比原始长读时的噪声还要小,且没有任何校正。
在酵母菌数据集上,经双色校正后的读数为246122个返回到参考基因组。这个数字至少比其他方法多4548.比LoRDEC的比对率80.672%,比CoLoRMap的比对率82.072分别高出2.7%和1.3%。双色也达到了最高的身份比97.969%,高于LoRDEC的身份比97.810%,比CoLoRMap的身份比96.515%高1.4%。
在第三个数据集果蝇中,双色校正后的长读序列比对参考基因组的读序列比对比对LoRDEC的少一些。它 可以看到Bicolor 实现 更高的比对比率 (0.2%) and 更高的一致性比率(1%).与CoLoRMap相比,Bicolor具有4413个更多的对齐读数,并且可以实现更高的标识率。我们注意到,与双色相比,CoLoRMap的比对比率可以提高2.1%比 Bicolor(37.544%的一致性比率)。可以看出,这个数据集有很多错误的基础,因为只有313989在901530个读断可以与参考基因组进行比对,原始数据的比对率相对较低(仅为37.079%)。这导致长期稳定的k-mers读数非常不可靠。此外,DBG中搜索的路径与预期的路径相去甚远。另一方面,CoLoRMap可以将短读和长读对齐,不依赖于k-mers。经双色校正后,与参考基因组比对的读序列更多,比对率低于彩色图谱。值得注意的是,双色达到了最高的同一性比率。
在典型的基因组覆盖情况下,所有方法的性能都非常接近(Table4)。仍然可以理解,CoLoRMap的性能最好。在大肠杆菌数据集上,所有方法均能达到100%的覆盖率。在数据集酵母菌上,CoLoRMap的表现略好于LoRDEC和双色酵母菌。具体来说,Bicolor和LoRDEC的覆盖率仅比CoLoRMap低0.012%和0.03%。此外,双色得到最低的覆盖率(即,93.915%)在果蝇的数据集上,与覆盖率分别为93.978%和94.735%的LoRDEC和CoLoRMap相比,覆盖率分别少了0.063%和0.82%。
image image.gif
3.3****从LoRDEC到iLoRDEC的性能改进
在第一层中,我们将LoRDEC打磨为迭代版本。这里我们比较了iLoRDEC和LoRDEC的性能。对不同参数的大肠杆菌数据集进行了实验研究。用四种不同k-mer初始长度(即,13;15;17和19),不同迭代轮数(从1到5)如Table 5所示。在[20]中,Salmela和竞争对手声称LoRDEC在默认参数下获得了最佳结果(见Table 3的第二行)。通过与LoRDEC的最佳结果比较,我们发现iLoRDEC在三种统计量下,在六组参数下的性能优于LoRDEC。值得注意的是,iLoRDEC总是能够获得更高的比对比。这些优越的性能验证了从LoRDEC到iLoRDEC的显著性能改进。
image image.gif
3.4****基于MSA的校正效果
在第一级校正后,我们得到n个预先校正的长读数。然后,在第二层使用工具肌肉来对齐这些类似的长读数。为了验证基于msa的校正的有效性,我们将iLoRDEC校正的结果与四种不同k-mer初始长度(即,n = 4)和5个不同的迭代轮数对数据集e.coli,以获得最终校正的长读。最终校正结果的比对统计如Table6所示。通过对比Table6和Table5中的对齐统计量,我们可以看出MSA校正后的对齐次数和比对比对比iLoRDEC的结果好很多。此外,恒等比非常接近Table5中最高的恒等比。结果表明,使用多组预校正长读来获得最终的校正长读可以提高性能。验证了MSA算法和投票算法的有效性。
3.5****优化时间成本的参数设置
该方法中k-mer的初始长度、迭代校正器的个数n和第一级圆数m是最重要的参数。其他四个参数,即, LoRDEC将实体k-mers的阈值、最大错误率和分支限制以及从LoRDEC继承的目标k-mer的数量设置为默认值(见Table2)。如果ki很大,许多长读就无法纠正,因为它们可能不包含任何可靠的k-mers。我们建议iLoRDEC使用的k-mer的初始长度应该小于LoRDEC使用的默认值。但是,较小的ki会导致更高的DBG复杂度,导致iLoRDEC的运行时间更长。按照LoRDEC的指示,k-mer的初始长度建议在集合内{13,15,17,19}用于细菌和真核生物的小基因组。对于大基因组物种,我们建议ki∈{13,15,17,19,21}。结果表明,当n =1时,双色膜的性能随着iLoRDEC的性能下降而下降。考虑到投票算法和运行时间, 我们 建议n>=3.
选择一个好的迭代轮数是很棘手的。Fig3显示了四种不同初始k-mers和五种不同迭代轮数下的比对比对比对的趋势(Table6)。从图中可以看出,当迭代轮数为2时,对直比可以达到最高水平。此外,随着迭代轮数从2轮增加到5轮,对准率降低。然而,当迭代轮数小于4时,对准率仍然相对较高(大于90%)。因此,当最佳迭代轮数小于4时,可以得到较好的对准比。该图还表明,单位比与迭代轮数成正比。这是因为迭代轮数越高,纠正的错误就越多。但是,当迭代次数设置大于3时,单位比并没有显著提高。因此,当迭代轮数大于2时,我们可以得到一个更好的恒等比。此外,当迭代轮数增加时,运行时间可以显著延长。因此,我们建议迭代轮数应小于5。预计校正结果应具有较高的比对率、较高的识别率和较低的时间消耗。在这项工作中,我们建议迭代轮数为3或4。
在第二个层次,我们在实验中设置了最快的选项—肌肉’-maxiters 1—diags’,因为计算MSA的时间复杂度很高。
image image.gif
image image.gif
3.6****运行时间对比
为了比较这些方法的运行时间,我们使用Linux/Unix time命令来记录实时时间。在我们的实验中,所有的核心都被用来运行程序。这些方法的运行时间见Table7。LoRDEC是最快的方法。由于CoLoRMap是基于映射的方法,因此它比LoRDEC慢。特别是OEA的过程非常耗时。双色包含两个计算阶段。第一阶段有许多iLoRDEC。尽管我们做了一些改进,但是预计运行时间比LoRDEC长很多倍。另一个原因是MSA的复杂性非常高。我们使用了最快的肌肉,但仍然花费了很多时间。双色跑得比CoLoRMap+OEA快。
image image.gif
4.****结论
本文介绍了一种用于PacBio长读纠错的双层框架。
在第一个层次,它利用不同长度的k-mers和迭代算法来确定多组初步校正的读数。
然后我们的方法结合这些初步的结果,基于MSA的校正在第二级。
通过对三个基准数据集的性能评估,证明了该方法与三种最先进的算法相比,能够达到最高的识别率。
通过对大肠杆菌和酵母数据集的比对,提高了比对率。
我们的方法也有一些缺点。
首先,在酵母菌和果蝇的数据集上有一小部分基因组覆盖丢失。
其次,运行时间比除OEA方法外的其他方法都长。
我们未来的工作将集中在这些方面,以提高速度。
参考文献:
Liu, Y, Lan, C, Blumenstein, M & Li, J 2019, 'Bi-level error correction for PacBio long reads', IEEE/ACM Transactions on Computational Biology and Bioinformatics.