第三代测序中基于德布鲁因图的长读错误纠正算法

第三代测序中基于德布鲁因图的长读错误纠正算法
摘要——PacBio单分子实时测序平台可以产生大量的长读序列，这对基因组的从头组装非常重要。尽管这些长读取具有15%的高错误率，但是由于它们的高错误率而放弃它们是不明智的。Illumina测序平台产生了长度在100 bp左右的短读，错误率低，成本低。但是组装形成的分支很多，不利于后续对基因组的分析。本文提出了一种新的混合纠错方法LecdB，它使用精确的短读来纠正具有较高错误率的长读。首先，从参考序列的短读构建两个固定长度和可变长度的de Bruijn图。然后，遍历要校正的长读数，以找到与固定长度的德布鲁因图一致的强K-mer。使用最大精确匹配算法将没有强K-mer的长读与可变长度de Bruijn图中的节点进行比对。实验表明，与其他基于de Bruijn图的长读校正算法相比，可以获得更好的结果。
关键词—德布鲁因图，纠错，通量，最大精确匹配
I.介绍
DNA测序确定四种含氮碱基(腺嘌呤[A]胞嘧啶[C]鸟嘌呤[G]胸腺嘧啶[T])的序列，是生命科学中的一项基础实验。作为第一代测序平台，Sanger [1]测序基于在体外DNA复制期间通过DNA聚合酶选择性掺入链终止双脱氧核苷酸。这种方法产生的读数具有接近99.999%的极高准确度，但长度相对较短，需要较长的测序时间和较高的成本。以Illumina为代表的第二代测序平台是基于Sanger测序技术开发的。高通量、短运行时间、低成本和低错误率是第二代测序技术的主要特点。第二代测序技术的缺点之一是产生的读长度很短，大约100-200个碱基。由结果组装而成的de Bruijn图有很多分支，无法确定哪一个是正确的路径，不利于后续的基因组分析和从头组装，无法有效利用。与第二代测序技术相比，以PacBio单分子实时(SMRT)测序平台和Oxford Nanopore Technologies的单分子测序平台为代表的第三代测序平台产生的平均读长度可以达到10kb-15kb，可以解决基因组组装、转录组重建、宏基因组学等现代生物学和医学的计算问题。但是第三代测序技术最明显的缺点是错误率可能高达15% [2]。
由于长读包含丰富的遗传信息，长读的纠错是目前最重要的任务。目前，纠错算法可以分为三种类型:长读自纠错；
短读取直接与长读取对齐，然后使用对齐信息进行纠错；
短读被组装形成与长读码对齐的重叠群。
对准信息用于误差校正。长读取自校正算法由LoRMA代表[3]。基于迭代的de Bruijn图，通过逐渐增加K-mer的长度来多次校正长读数。由于不使用参考短读数，需要更深的测序深度，这导致数据集过大，对计算机硬件的要求相对较高。PacBioToCA [4]是长读取误差校正算法的代表，其将短读取与长读取对齐。在短读与长读比对之后，在长读的重复区域中，仅保留具有最高相似性的比对结果用于纠错。但占用大量内存和硬盘空间，运行时间长，不符合当前对测序速度的要求。
基于将短读组装成重叠群的长读纠错算法可以由de Bruijn图算法来表示。de Bruijn图算法最早由Leena Salmela提出，并开发了纠错工具LoRDEC [5]。短读用于构建德布鲁因图，然后映射到长读进行纠错，大大提高了纠错速度。Jabba [6]在此基础上打破了种子长度的限制，使用伪匹配的方法，通过使用最大精确匹配算法来寻找种子，这样种子可以更长，在将长读数与de Bruijn图对齐时对齐精度更高，纠错结果更准确。FMLRC [7]在图形创建过程中使用FM-index来动态选择K-mer大小，并构造两个德布鲁因图，以便可以进行两次误差校正。
我们提出了一种基于德布鲁因图的混合纠错方法。构造了两个德布鲁因图。首先，在长读中搜索强K-mer，然后使用具有固定长度K-mer的de Bruijn图来纠正错误。第二，具有可变长度K-mer的de Bruijn图用于那些没有找到强K-mer的长读取的错误校正。这种方法可以保持高通量和高比对一致性，并尽可能减少运行时间。
II.德布鲁因图的基本原理
在生物医学领域，长度为K的弦称为K-mer。K-mer的具体解释如下:设r是一串碱基，其字母表为∑ 𝜖 {𝐴，𝐶，𝐺，𝑇}.从r中取一个长度为k的子串α，α= r[I]r[I+1]⋯r[I+k1]，0≤I ≤| r | k+1，我们称α为r中的K-mer，K-mer的长度有固定值和可变值两种形式。当K为固定值时，长度为r的字符串可以产生r-k+1个K-mers。当K的值可变时，K-mer是不同长度的字符串。de Bruijn图是表示序列之间重叠关系的有向图。图中的每个节点都是一个K-mer。基因组序列中所有相同的K-mer记录为图上的一个节点，相同K-mer的重复次数记录为节点上的丰度。序列中不同的K-mers记录在不同的节点中。如果前一个K-mer的后缀与下一个K-mer的前缀重叠一定长度，则在两个节点之间画一条有向线以形成一条边。根据K-mer的长度，de Bruijn图有两种形式:(1)K-mer的长度是固定长度。图1(a)显示了K=3的德布鲁因图的一个例子，其中K-mers的长度相同。由于K-mers的长度是固定的，两个K-mers重叠K-1个碱基就足够了。(2)K-mers的长度是可变的。图1(b)显示了一个德布鲁因图的例子，其中K-mers的长度是不同的。由于K-mer的长度是可变的，所以对于de Bruijn图的构建只设置了最小值K，并且每个K-mer的长度需要大于K-1，并且两个相邻节点中的K-mer需要重叠K-1个碱基。
de Bruijn图可用于生物信息学中基因组的从头组装。基因组的从头组装不使用任何参考序列，直接通过reads构建德布鲁因图，然后对德布鲁因图进行处理，例如消除短的死端和其他结构，或者基于

(a)(二)
图一。(a)具有固定长度的K-mer和K = 3的De Bruijn图(b)具有可变长度的K-mer和K = 5的De Bruijn图。

关于图的拓扑结构。这是在没有参考序列的情况下对物种基因组进行测序的唯一方法。
II.基于德布鲁因图的混合纠错算法
A.纠错的第一步
首先，从输入参考序列的短读构建固定长度的de Bruijn图。对于要校正的序列的每个长读，从长读的第一个碱基中提取具有固定长度K的K-mer。值得注意的是，从长读数中提取的K-mer的K值应该等于在构建德布鲁因图中使用的K值。例如，当K=19时，我们从长读数中提取前19个碱基，然后在de Bruijn图中搜索相同的K-mer。如果有相同的K- mer，将其标记为强K-mer，并假设强K-mer是正确的。如果没有找到相同的K-mer，则该K-mer被标记为弱的。然后向前移动一个碱基，继续提取K-mers进行匹配，直到遍历完长read，标记出所有强K-mers。如果在长读取中没有发现强K-mer，则该长读取将在第二个纠错步骤中被处理。如果找到了强K-mer，则从长读的第一个碱基到第一个强K-mer的区域被定义为头部区域，从长读的最后一个强K-mer到最后一个碱基的区域被定义为尾部区域，并且头部/尾部中间的区域被定义为内部区域。根据位置的不同，不同地区的处理步骤也不同。在图2中，矩形中的黑点是强K-mers，左右区域是头部和尾部区域，中间区域是内部区域。
1)内部区域的误差校正
对于内部区域错误校正，将找到的第一个强K-mer作为源节点，将下一个强K-mer作为目标节点，并将源/目标节点和最大分支限制作为输入参数，以在de Bruijn图中找到从源节点到目标节点的最短路径。这是将长读序列转化为德布鲁因图中相似序列所需碱基修饰最少的途径。在字符串计算中，编辑距离是将字符串A转换为B所需的编辑次数，用于衡量两个字符串之间的相似性。因此，最短路径算法也被称为最小编辑距离算法。使用动态编程算法可以找到最小编辑距离的路径。
本文根据de Bruijn图上的深度优先策略，从源节点到目标节点进行搜索，计算两个节点之间区域的最小编辑距离，并存储在动态规划矩阵中。这

图二。矩形内的黑点代表强K-mer，矩形外的区域代表弱区。薄弱区域可分为头部、尾部和内部区域。

当达到允许的最大编辑距离或到达图表的终端时，路径搜索将立即停止。当达到最大分支数量时，它也将停止探索路径。如果成功找到最小编辑距离对应的路径，记录下来，用路径上的碱基替换长读上对应的内部区域，实现纠错。如果没有找到路径，则在源/目标节点之间添加长度与源/目标节点之间的长度相同的虚拟边，从而确保从第一个强K-mer到最后一个强K-mer存在路径。然而，在三种特殊情况下，我们不搜索源/目标节点之间的路径。
a)当源目标节点对应于相同的强K-mer时:它们之间的区域被认为是正确的，不执行路径搜索。
b)当它们具有相同的强K-mer并且相隔一定距离时:这可能是由于读中的重复或错误K-mer造成的，我们选择跳过这种情况。
c)在长读中，源节点和目标节点之间的距离过大:动态规划矩阵需要过多的内存来计算最小编辑距离，消耗了过多的计算资源，我们也选择跳过这种情况。
2)头/尾区域的误差校正
对于头部区域纠错，首先在长读中提取从第一个碱基到第一个强K-mer的子序列，然后反向整个序列，使强K-mer在起始位置。所以纠错方法和尾区的子序列是一样的，我们就只设计尾区的纠错算法。我们使用尾部区域的强K-mer、该区域中的子序列和最大分支极限作为算法输入参数。这里把长read中最后一个强K-mer作为源节点，没有目标节点，只能通过深度优先搜索找到编辑距离最小的路径，把尾部区域的长度作为可能的延伸长度。这个过程允许找到从强K-mer开始的任何路径，并且可以通过优化强K-mer的前缀长度来找到最佳扩展，直到它到达长读取的末尾或者到达最大分支限制，或者到达de Bruijn图的死胡同。
B.纠错的第二步
对于所有需要纠正的长读，在第一步之后，有两种情况:(1)我们已经成功地找到了强K-mers，并完成了长读的纠错。我们认为这些长读是正确的，不会继续处理。
(2)在长读中没有发现强K-mer。由于我们没有强K-mers作为校正基础，我们无法在第一步中完成校正过程。对于校正过程的第一步，我们寻找由长读数和具有固定长度K-mer的de Bruijn图共享的K-mer。因此，通过减少K-mer的长度，可以构建其节点中较短K-mer的de Bruijn图，以找到更小的强K-mer，并且可以重复校正过程的第一步。但是如果K-mer的长度太短，重复的K-mer会增加很多，并且边缘
节点之间也会增加。在寻找最小编辑距离时，很难找到正确的路径。此外，每次改变K-mer的长度，都必须重建de Bruijn图，这将消耗大量的内存和时间，因此，我们不得不考虑另一种方法来代替。
首先，参考序列的短读数用于构建可变长度的de Bruijn图[8]。构建图时，只设置最小K-mer长度，构建最小K-mer长度de Bruijn图。然后，处理德布鲁因图以找到相邻节点之间的唯一路径。如果发现这一点，则合并这两个节点，并且这两个节点中的K-mer将被组合以形成更长的K-mer。所有节点都需要相同的处理，直到任何一对相邻节点之间都没有唯一的路径。因此，构建了具有不同K-mer长度的de Bruijn图。然后，使用最大精确匹配算法在de Bruijn图的节点中寻找长读取和K-mers之间的匹配子序列。这些匹配的序列被用作种子，并且基于种子找到de Bruijn图中种子之间的最佳路径。然后根据种子的顺序将长读对齐到这个最优路径，完成长读的纠错。

图3。(a)图中的黑点代表德布鲁因图中的每个节点，下面的长线段代表长读数，短线段代表德布鲁因图和长读数之间的最大精确匹配。(b)长read和de Bruijn图的最终比对路径。

1)使用最大精确匹配(MEM)算法寻找种子
最大精确匹配表示两个序列之间完全匹配的字符串。最大精确匹配称为种子。将对应于de Bruijn图中每个节点的碱基序列和它们的反向互补序列连接以形成重叠群，然后可以找到重叠群和要校正的长读之间所有可能的最大精确匹配。为了找到所有的种子，我们通过使用重叠群构建了一个增强的稀疏后缀数组(ESSA[9])。ESSA构造的索引字符串在字典中排序，因此在长读中匹配它们的子序列可以被快速找到。构建ESSA后，逐一输入第一步中未校正的长读数。对于每个要纠正的长读取，遍历它以在ESSA中找到具有最小种子长度的匹配前缀。如果在ESSA中找到前缀，则将其记录为种子，然后逐渐增加前缀长度，直到达到长读取的最大长度。详细的搜索过程可以在[9]中找到。遍历长读数后，记录找到的所有种子以供后续分析。图3(a)显示了德布鲁因图和长读数之间的最大精确匹配。德布鲁因图的节点中的种子可以分为四种情况:(1)节点对应的子序列是种子。(2)对应于一个节点的子序列的一部分是A种子。(3)在长读中，节点中的子序列包含多个种子，并且它们可以相隔一定的距离。(4)长读中的种子对应于de Bruijn图中两个或多个节点中的子序列。
1)种子和对齐
经过最大精确匹配算法的处理后，我们记录下所有的种子，然后根据种子在长读中的顺序，找出在de Bruijn图中是否存在种子顺序相同的序列。如果有这样的序列，我们认为长读与序列对齐。首先，根据局部对齐来考虑和分析de Bruijn图的每个节点中的种子。如果在de Bruijn图的一个节点中有多个种子，并且它们的顺序与长读的一段中的顺序一致，则该节点中的种子被连接以形成与长读的段的不准确匹配，并且该节点被认为是正确的。认为长读数中的部分与de Bruijn图对齐。其次，我们考虑德布鲁因图中不同节点之间的种子。根据全局对齐，我们首先检查是否所有的种子只存在于它们对应的节点中，即其他节点没有相同的种子。然后我们检查种子是否跨越两个或更多的节点。如果种子只存在于它们对应的节点中，并且不跨越两个或更多个节点，则认为种子是可靠的。然后按照种子在长读中出现的顺序，连接这些种子对应的节点。如果存在唯一的路径，则该路径被认为是正确的，并且长读数中的部分已经与de Bruijn图中的路径对齐。如果有多条路径，则考虑每条路径的边长，即一条路径中从开始节点到结束节点的序列的总长度。我们选择具有与长读数对齐的最短边长的路径。因为边的长度越短，相应节点中的种子就越可靠。如果没有找到合适的路径，则认为德布鲁因图是错误的或者长读数是错误的

图4。LecdB工作流。

纠错结果和参考基因组之间的相似性，使用索引标识。为了评估纠错结果中的碱基数量与原始长读取中的碱基数量的比率，使用了索引通量。其他评价指标还有平均长度、实际运行时间、CPU时间。
对于Jabba，在M_zebra(表1)和parrot(表2)的实验中，aligned和identity两个指标都高于0.89，是所有方法中性能最好的。因为Jabba使用最大精确匹配来寻找种子。对于长读取，找到更长的种子意味着de Bruijn图中的对齐路径更可靠，并且纠错结果的质量会更高。在M_zebra(表1)和parrot(表2)的实验中，LecdB的比对和同一性都高于LoRDEC，但差别不大。这是因为LoRDEC和LecdB在校正过程的第一步都采用深度优先搜索算法，只能达到局部最优解。因此，在探索过程中，可能由于正确分支得分低而探索另一条路径，导致纠错结果质量差。在纠错过程的第一步之后，LecdB再次尝试纠正未纠正的长读，所以评价指标会略高于LoRDEC。
然后我们分析通量。从表中可以看出，LoRDEC在M_zebra(表1)和parrot(表2)实验中的通量高于其他工具。LecdB的通量略低于LoRDEC，但相差不大。值得注意的是，Jabba在M_zebra(表1)和parrot(表2)实验中的通量是最小的，这是因为给定的K值太高，导致图的分支较少。找到种子后不容易扩展路径或者因为路径太长而放弃纠错，所以通量最低。
平均长度方面，在M_zebra(表1)和parrot(表2)实验中，LecdB的平均长度最高，LoRDEC的平均长度相对低于LecdB，Jabba的平均长度最短。
就程序运行时间而言。在M_zebra(表1)和parrot(表2)实验中，LoRDEC的运行时间最短，LecdB的运行时间是LoRDEC的两倍，Jabba的运行时间最长。由于Jabba是多线程的，由于大量的增强稀疏后缀树，要找到精确匹配需要更多的时间。LecdB在第二个纠错步骤中有更少的数据，因此找到精确匹配的时间更短，它比Jabba执行得更好。
III.结论
提出了一种基于德布鲁因图的混合纠错方法LecdB。首先，利用精确的短读构造固定长度和可变长度的两个de Bruijn图。然后，通过遍历长读取，将固定长度的de Bruijn图中的相同K-mer搜索为强K-mer。基于强K-mer，采用最小编辑距离算法进行第一步纠错。接下来，使用最大精确匹配算法将没有强K-mer的长read与变长de Bruijn图进行精确匹配，并将精确匹配作为种子。最后，根据用于纠错的种子顺序，将长读数与德布鲁因图中的节点对齐。实验结果表明:LecdB改善了通量不足的问题，并且索引对齐和相似性也比较高。