HG-CoLoR用一个变阶de Bruijn graph混合校正高噪声长读数

HG-CoLoR用一个变阶de Bruijn graph混合校正高噪声长读数

1 简介

几年来，长读长测序技术不断发展，可以解决大型复杂基因组的组装问题，在此之前，仅使用短读长测序技术很难解决这些问题。这些长读长测序技术的两个主要参与者是 Pacific Biosciences 和 Oxford Nanopore。后者随着 MinION 设备的发布，可以在一台简单的笔记本电脑上运行，实现了低成本和简单的长读取测序。

尽管长读长可以达到数十 kbps 的长度，但它们也达到了非常高的错误率，太平洋生物科学公司约为 10-15%，牛津纳米孔的错误率高达 30%。由于这种高错误率，在下游问题中使用它们之前纠正长读取是强制性的。许多方法可用于短读纠错，但它们不适用于长读，主要有两个原因。首先，这些方法不能扩展到长读取的更高错误率。其次，这些方法中的大多数都关注替换错误，这是 Illumina 数据中的主要错误类型，而插入和删除在长读中更为频繁。

1.1 相关作品

最近，已经开发了几种用于长读纠错的方法。这些方法可以分为两大类：要么长读长通过相互对齐来进行自我校正 [PBDAG-Con ( Chin et al. , 2013 ), PBcR ( Berlin et al. , 2015 )]，要么要么采用混合策略，使用互补的短读。在这种情况下，短读可以与长读对齐 [Nanocorr ( Goodwin et al. , 2015 ), CoLoRMap ( Haghshenas et al. , 2016 )]，也可以组装成与长读对齐的 contig [ HALC（鲍和兰，2017 ）]。 De Bruijn 基于图的方法，通过遍历图的路径来纠正长读取，最近也开始发展，在混合情况下 [LoRDEC ( Salmela and Rivals, 2014 ), Jabba ( Miclotte et al. , 2016 )] ，以及在非混合情况下 [LoRMA ( Salmela et al. , 2017 ), Daccord (Tischler and Myers, 2017, unpublished)]。 NaS ( Madoui et al. , 2015 ) 不使用短读来纠正长读，而是使用长读作为模板来招募短读并将它们组装成重叠群，用作更正序列。这种方法需要将短读取与长读取对齐，以便找到种子，这些种子是与长读取对齐的短读取。然后将种子与所有其他短读取进行比较，以招募新的短读取，对应于长读取的低质量区域。

1.2 贡献

我们介绍了 HG-CoLoR，这是一种新的长读混合纠错方法，它结合了最先进的两种思想：短读与长读的对齐（如 CoLoRMAp），以及使用 de Bruijn图，从短读（如在 LoRDEC 和 Jabba 中）构建。然而，与这些方法不同，HG-CoLoR 使用可变阶 de Bruijn 图，而不是经典图。因此，HG-CoLoR 专注于种子和扩展方法，其中通过将短读取与长读取对齐找到的种子用作可变阶 de Bruijn 图上的锚。然后遍历该图，以便将种子链接在一起，从而纠正短读取未覆盖的长读取区域。我们的实验表明，与最先进的混合和非混合长读纠错方法相比，HG-CoLoR 在运行时间和结果质量之间提供了最佳平衡，并且是唯一能够有效地扩展到真核基因组。他们还表明，HG-CoLoR 纠错的效率意味着令人满意的装配结果。

变阶de Bruijn图

2.1 de Bruijn 图

de Bruijn 图是一种在组装工具中广泛使用的数据结构。它的节点被定义为读段的 k -mers，它的边表示由节点表示的 k-mers 之间长度为 k-1 的前缀-后缀重叠。然而，尽管它很有用，但众所周知，de Bruijn 图面临困难，因为 k -mer 的大小在构建时是固定的。一方面，选择较大的 k 可以让图更好地处理重复区域，但会导致覆盖不足区域的边丢失。另一方面，选择较小的 k 将允许在覆盖不足的区域中正确检索图的边，但在重复区域中会导致更多的分支，从而带来更多的困难。

为了克服这些问题，现代组装器通常构建多个不同阶的 de Bruijn 图。虽然这种方法可以提高组装集的质量，但它也大大增加了运行时间和内存消耗，因为需要构建和存储多个图。

最近，一些方法是制定允许代表权的所有de Bruijn图，最多为K，在一个单一的数据结构。在单个数据结构中例如，流形 de Bruijn 图 ( Lin and Pevzner, 2014 ) 将任意子串与节点相关联，而不是关联 k -mers。然而，这种结构主要具有理论意义，因为它尚未实施。提出了可变阶 de Bruijn 图的另一种实现 Boucher 等人（2015 年）。它依赖于 Bowe 等人对 de Bruijn 图的简洁表示。 (2012) ，并支持允许动态更改图形顺序的其他操作。但是，当前的实现只支持构建最多 64 个，这太严格了，因为我们不想限制最大的顺序。

因此，我们引入了可变阶 de Bruijn 图的新实现。它依赖于 PgSA ( Kowalski et al. , 2015 )，一种允许回答关于一组读取的各种查询的索引结构。

2.2 PgSA 概述

PgSA 是一种数据结构，它允许对一组读取进行索引，以回答给定字符串 f ：

在哪些读取中 f ？

在多少次读取中 f ？

的出现位置是 f 什么？

出现的次数是 f 多少？

在哪些读取中 f 只出现一次？

在多少次读取中 f 只出现一次？

的出现位置是什么 f 在 reads 中只出现一次

在这些查询中， f 可以作为 DNA 符号序列给出，也可以作为一对数字给出，分别表示读取 ID 和 f 在该读取中的起始位置。

如前所述，为了回答这些查询，必须建立读取索引。 PgSA 如下构建它。首先，将所有具有重叠的读数与这些重叠连接起来，以获得假基因组。如果在创建假基因组后留下了一些没有发现重叠的读数，它们会在其末尾简单地连接起来。然后，计算伪基因组的稀疏后缀数组，以及允许从伪基因组中的原始集合检索读数的辅助数组。该辅助数组的每条记录都将原始读取集合中的读取 ID 与伪基因组中的读取偏移相关联，并且还包含标志数据，这些数据带来关于读取的补充信息，用于处理查询。查询是通过对后缀数组进行简单的二进制搜索来处理的，并使用此补充信息。

由于在假基因组计算期间读数重叠，并且由于 PgSA 不记录任何有关其长度的信息，它只允许索引和查询一组恒定长度的读数。但是，查询字符串的长度不是在编译时设置的，因此 PgSA 支持 f 对可变长度

2.3 变阶de Bruijn图表示

最大阶 K ，并且读取的 K -mers 用 PgSA 索引，以便能够表示所有 de Bruijn 图的节点，直到这个最大阶。阶的 de Bruijn 图，给定节点的边 k ≤ K 是通过查询索引来检索的，使用第三个查询（即 f 什么？），后缀长度为 k - 1 的 k -mer。查询返回一组数字对，每对代表一个 K -mer ID 和查询字符串在该 K -mer 中的出现位置。然后处理这些对，并且仅保留那些位置分量不代表 k -1 的 K 那些（以便可以将出现向右扩展为 k -mer）。这些扩展的出现表示 k 的前缀-后缀重叠的 k - 1 。 mer 与由当前考虑的节点表示的 k-mer，因此定义了该节点的边缘

由于通过查询索引来检索边，因此向后遍历图也很容易。对于给定的订单 k ，不是使用节点表示的 k -mer 的后缀来查询，而是简单地使用它们的前缀来查询索引。然后以与前向遍历相同的方式处理返回的对集合，除了仅保留其位置分量不代表 k -1 的 K 来定义边。的 de Bruijn 图中检索任何给定节点（向前或向后）的边的 k ≤ K 中给出了允许在补充算法 S1 。

3.1 概述

如前所述，HG-CoLoR 结合了最先进的两个想法：短读取与长读取的对齐，以及使用具有可变顺序的特殊性的 de Bruijn 图。为此，它专注于种子和扩展方法，通过将短读取与长读取对齐来发现种子。种子被定义为一个 5 元组（ id 、 pos 、 len 、 score 、 seq ），其中： id 是与种子关联的长读的 id ， pos 是长读上对齐的开始位置， len 是比对的长度， score 是比对的分数， seq 是在这个位置与长读比对的短读的实际共有序列。一旦种子被检索到，它们通过扩展它们的序列连接在一起，在前面描述的可变阶 de Bruijn 图的帮助下。该图是从短读中构建的，通过选择最大顺序 K 并用 PgSA 索引它们的 K -mers，并通过查询索引进行遍历，如前所述。对于每次长读，遍历图以将相关种子链接在一起，这些种子用作锚点。因此，将两个种子链接在一起所遵循的图形路径指示了长读的未覆盖区域的校正序列。最后，一旦所有的种子都链接起来，通过进一步遍历图来扩展获得的序列的尖端，以到达原始长读的末端。 HG-CoLoR 的工作流程总结在图 1 及其四个主要步骤如下所述。

图。1。

在新标签中打开下载幻灯片

HG-CoLoR 的工作流程。首先，使用 QuorUM 纠正短读取，以尽可能多地消除测序错误。然后， K 为该图选择 K KMC3 获得来自校正的短读数为了进一步降低错误率，对校正的短读数应用过滤步骤，并 K 去除出于同样的原因，只有来自校正的短读数的固体 K -mer 用 PgSA 进行索引，以表示可变阶 de Bruijn 图。然后在 BLASR 的帮助下将先前过滤的更正短读取与长读取对齐，以找到种子。然后独立处理每个长读取。对于它们中的每一个，遍历该图以将相关种子链接在一起，用作锚点，以便检索长读取的未覆盖区域的校正序列。然后，将所有种子链接在一起后获得的序列的尖端通过遍历图在两个方向上扩展，以到达初始长读的末端。最后输出校正后的长读

是 LoRDEC）具有高度相似性，但使用源自种子的序列作为图上的锚点与尤其（使用来自长读。事实上，在高度错误的长读取的情况下，即使是短的、可靠 k -mers 也很有可能出错。这种错误 k -mer 将导致使用错误的锚点，从而导致不满意的校正结果。然而，由于短读取是准确的，它们产生的种子可以用作可靠的锚，几乎没有错误的机会。此外，使用这些种子作为锚点还允许直接构建具有较大 k ，而无需执行多轮校正，并且在每一步都增加 k ，与 LoRMA 的方式相同。

3.2 短读校正和图构建

尽管在进行任何校正之前短读数已经准确，但它们仍然包含一小部分错误。由于 HG-CoLoR 试图从短读取中构建高最大阶的可变阶 de Bruijn 图，因此必须从该数据中删除尽可能多的错误，以避免图中的错误路径。为此，在 QuorUM (Marçais et al., 2015) 的帮助下纠正了短读，在我们测试的所有短读纠错工具中，它提供了运行时间和纠正质量之间的最佳权衡。

最大阶数 K 然后为该图选择 Kokot KMC3 从校正的短读数中提取 K-mers ( et al. , 2017 )。为了进一步降低short reads的错误率，从而避免图上的错误种子和嵌合路径，包含弱 K -mers（即用于出现小于某个阈值的K-mers）的short reads被过滤掉，不以下步骤，并且仅使用实心 K -mers 来构建图形。

3.3 种子检索和合并

为了检索种子，在 BLASR 的帮助下，短读取与长读取对齐（ Chaisson 和 Tesler，2012 年）。尽管它最初旨在将由插入和删除错误主导的长读取与参考基因组对齐，但在将短读取与长读取对齐时，BLASR 也设法很好地处理了此类错误。然后独立处理每个长读取，并对相关种子应用两个分析阶段。

首先，种子是成对处理的，按照它们在长读上的对齐起始位置的升序排列。如果给定的一对连续种子的对齐位置和对齐长度表明它们重叠了足够长的长度，则比较它们假定的重叠序列，并相应地合并两个种子。如果假设的重叠序列不重合，或者如果比对位置和比对长度表明两个种子确实重叠，但长度不足，则仅保留具有最佳比对分数的种子。第一个合并阶段的算法在补充算法 S2 。

其次，一旦所有具有重叠比对位置的种子都被合并或过滤掉，这些种子再次被成对处理，按照它们在长读上的比对起始位置的升序，计算具有接近比对位置的连续种子之间的序列重叠. 与上一步一样，如果给定的种子与下一个种子完全重叠足够长的长度，则将两个种子合并。此步骤允许考虑在对齐步骤期间未检测到的长读取中的小插入错误，这可能导致下一步中的链接困难。第二个合并阶段的算法在补充算法 S3 。

我们选择始终对齐所有短读取，而不是仅对齐其中的一个子集，因为这种方法可以发现更多种子。因此，可以纠正更多的长读，以及这些长读的更大区域，并且需要探索图的较小部分。此外，对齐步骤所花费的时间比例平均仅占整个方法运行时间的 20-40%，并且随着数据集大小的增加而趋于减少。因此，仅对齐短读取的一个子集而不是对齐它们都不会显着减少该方法的整体运行时间。

3.4 种子链接

一旦为所有长读找到并合并了种子，HG-CoLoR 就会独立处理每个长读，并尝试通过将它们视为对并遍历图来将它们相关联的种子链接在一起。对于给定的对，具有最左侧对齐位置的种子称为源，具有最右侧对齐位置的种子称为目标。为了将一对种子链接在一起，源的最右边 K -mer 和目标的最左边的 K -mer 被用作图上的锚点。然后遍历该图，以找到两个锚点之间的路径。当找到这样的路径时，它所指示的序列被用作对长读的未覆盖区域的校正。搜索两个种子之间的路径首先从源到目标执行，如果找不到路径，则再次执行搜索，从目标到源。搜索是双向进行的，因为根据遍历的起点，可能会探索到图的不同部分，从而导致不同的遍历。

HG-CoLoR 从最高阶开始遍历可变阶 de Bruijn 图。仅当该节点没有当前订单的任何边时，或者如果其当前订单的所有边都已被探索并且不允许到达目标，则该订单在给定节点处减少。当图的阶数减少时，来自源和来自目标的 k -mers 的大小也相应减小，因此它们仍然可以用作锚点。还设置了最小顺序，因此 HG-CoLoR 不会遍历表示短且可能无意义的重叠的 de Bruijn 图。

当面对给定阶 数 k ，HG-CoLoR 执行贪心选择。的节点的边 k 总是首先探索通向代表具有最高出现次数太多节点 k ，尽管有校正和过滤步骤，但可能包含测序错误。阶的 de Bruijn 图探索边之后 k < K ，由于没有找到较大阶的边，因此 K 在继续遍历之前，为了避免探索太多的分支路径。

当找到两个种子之间的路径时，由于贪心选择和图的阶数仅局部降低的事实，它被认为是最优的。因此选择它作为长读的未覆盖区域的校正。我们自愿以这种方式选择最佳路径，而不是探索多个路径并选择与长读取对齐的最佳路径作为校正，以避免过高的运行时间。

我们还设置了一个不匹配阈值 t 当将两个种子链接在一起时，，则源和目标可以链接在一起 K 源（或目标） K 少于 t 不匹配 K -mer （分别来源）。尽管有校正和过滤步骤，这样的阈值允许克服仍然存在于种子上的少数错配错误。然而，不能正确处理两个种子包含错误的情况。

此外，HG-CoLoR 可能会尝试将出现在图中未连接区域的两个种子链接在一起。例如，如果源自另一条染色体的短读段错误地与长读段对齐，则可能会出现这种情况。在这种情况下，将这些种子连接在一起是不可能的，因为它们之间不存在任何路径。因此，为了避免在这种情况下对图进行昂贵的探索，设置了对分支探索的最大数量的限制。如果达到此限制，并且没有找到链接源和目标的路径，则放弃当前的链接迭代，并且 HG-CoLoR 尝试跳过无法达到的目标。换句话说，源保持不变，无法到达的目标被忽略，目标被重新定义为以下种子，并执行新的链接迭代。中给出了该过程的说明补充图 S1 。

然而，由于跳过种子可能导致大量失败的链接尝试，即不可能在允许的最大分支探索数内将两个种子链接在一起，因此还设置了对可以跳过的最大种子数的限制。一旦达到此限制，由于没有任何连接尝试成功，HG-CoLoR 会在源和跳过的第一个种子之间填充长读取的未覆盖区域，并使用原始长读取的碱基。然后再次遍历该图，以便将剩余的种子链接在一起，从第一个被跳过的种子开始，并纠正长读的剩余部分。

3.5 提示扩展

最后，很明显，种子并不总是在长读的开头和结尾对齐。因此，为了尽可能接近其原始长度，一旦给定长读的所有种子都已链接，HG-CoLoR 继续遍历图以扩展产生的校正长读的尖端。以与上一步相同的方式，遍历从可变阶 de Bruijn 图的最高阶开始，并且仅当给定节点没有当前阶的任何边时，该阶才在给定节点处减少。校正的长读数的尖端因此被延伸，直到到达原始长读数的末端或分支路径。实际上，在提示扩展的情况下，当面对分支路径时，HG-CoLoR 不知道选择哪条路径并继续扩展，也没有任何锚点，这与它试图将两个种子链接在一起时不同。因此，多个分支的贪婪选择和探索是没有用的，出现这种情况时就干脆停止扩展。在这种情况下，校正后的长读物的末端从原始长读物用碱基进一步延伸，直到到达原始长读物的末端。此外，我们选择在到达原始长读的末端时始终停止扩展，即使在可以进一步遍历图的情况下，因为我们寻求执行纠错而不是组装。进一步扩展长读取确实会使该方法更接近重叠群组装过程，而不是实际的纠错管道。

3.6 输出

填充种子未覆盖的区域，并用原始长读数的碱基扩展校正的长读数的尖端，可以保持长读数的连通性，这对于脚手架等应用非常有用。对于校正碱基的重要性占主导地位的其他应用，还提出了经典的修剪和/或拆分输出。实际上，在校正过程结束时，如果一个已校正长读长的碱基来自种子或可变阶 de Bruijn 图，则它可以被认为是已校正的，如果它来自原始长读长，则可以认为它是未校正的。 HG-CoLoR 以大写字母输出校正的碱基，以小写字母输出未校正的碱基。因此，HG-CoLoR 可以输出每个更正的长读取的三个版本。第一个版本，所有未校正的碱基都被保留，一个修剪版本，其中从长读取末端的未校正碱基被移除，以及一个拆分版本，其中所有未校正的碱基被移除，并且仅长读取的区域可以纠正的作为单独的序列输出。

4 结果与讨论

我们对 Pacific Biosciences 的模拟数据和真实的 Oxford Nanopore 数据进行了实验。在这两种情况下，都包括了来自 A.baylyi 、 E.coli 和 S.cerevisiae ，并使用了互补的真实短 Illumina 读数。对于真正的 Oxford Nanopore 数据，包括来自 秀丽隐杆线虫 。由于该生物体没有质量令人满意的真正 Illumina 读数，因此使用 ART 模拟了互补的 Illumina 短读数（ Huang 等人，2012 年）。的帮助下获得了模拟的 Pacific Biosciences 长读数，这是 Stöcker et al. , 2016 一种基于真实 Pacific Biosciences 数据开发的错误模型的最先进的模拟器。使用以下参数来获得大约 15% 的错误率： –prob-ins 0.22、–prob-del 0.08 和 –prob-sub 0.02 。中给出了不同数据集的所有详细信息补充表 S1 。

我们将 HG-CoLoR 与混合纠错工具 CoLoRMap、HALC、Jabba、LoRDEC、Nanocorr 和 NaS 以及三种自纠错工具 Daccord、LoRMA 和汇编器 Canu 中使用的方法进行比较（ Koren 等人， 2017 年）。由于运行时间长，NaS 只能在快速模式下运行。由于 Jabba 和 Daccord 等方法只提出了拆分输出，我们还拆分了所有其他工具的输出，标记了未更正的碱基，以进行公平比较。除非另有说明，所有实验均在配备 16 核的 32 GB RAM 机器上运行。

4.1 参数

上使用 HG-CoLoR 进行了多轮校正， 酿酒酵母 以试验参数。因此，我们发现使用最大阶 K = 100 的可变阶 de Bruijn 图在运行时间、校正的长读取数、拆分长读取的比例、平均长度和校正碱基数之间产生了最佳折衷（参见补充图 S2 ）。对于两个合并阶段，允许合并两个种子的最小重叠长度设置为 99，相应于 K。为图表选择对于第二个合并阶段，两个连续种子之间的对齐位置的最大距离设置为 10。该图的最小顺序设置为 k = 40，因为将其设置为较大的值会导致更大的拆分比例，因此由于局部覆盖率下降，因此长读取时间更短。将其设置为较小的值还会导致更多的拆分和更短的长读取，因为除了更大的运行时间之外，由于探索了无意义的边缘，尤其是在重复区域中（参见补充图S3 ）。分支探索的最大数量设置为 1250，因为减少它也会导致更多的拆分和更短的长读取，并且增加它几乎不会产生更好的结果，但会增加运行时间（参见补充图 S4 ）。出于类似的原因，种子跳过的最大数量设置为 5，错配阈值设置为 3。对于短读取与长读取的对齐，BLASR 使用默认参数，但 bestn 设置为50 而不是 10，以获得更多的种子，从而纠正更多的长读。再一次，将这个参数增加到更大的值只会影响运行时间，并没有显着改善校正结果，而减小它会导致校正的长读取数量下降。由于我们在实验中仅使用了 50 倍的短读取覆盖率，因此将 K -mer 固体阈值设置为 1（即所有 K -mer 都被视为固体）。最后两个参数高度依赖于输入短读取的特性。特别是， bestn 参数应在使用高于 50× 的短 read 覆盖率时降低，而 solid 参数应在使用高于 50× 的短 read 覆盖率时增加，或者在使用非常准确的短 read 时，显示错误率较低超过 1%。 Canu 使用参数 -correct, stopOnReadQuality=false ，由于长读数的高错误率， corOutCoverage = 300 ，以便尽可能多地纠正长读数，并且将基因设置为每个参考的确切碱基数基因组。此外， -nanopore-raw 用于校正真正的 Oxford Nanopore 长读数， -pacbio-raw 用于校正模拟的 Pacific Biosciences 长读数。其他工具使用默认或推荐参数运行。为了进行更好的比较，在运行 Jabba 之前使用 QuorUM 纠正了短读取，而不是使用作者推荐的工具 Karect ( Allam et al. , 2015 )。所有工具都使用 16 个进程运行。

4.2 模拟数据对比

为了评估不同工具的准确性，我们首先在模拟的 Pacific Biosciences 长读数上对其进行了测试。因此，我们能够知道在哪些位置引入了哪些错误，并精确测量了校正的准确性。为此，我们使用了 LRCStats ( La et al. , 2017 )，这是一款专门用于测量模拟数据准确性的软件。上的所有工具的输出统计 A.baylyi 和 E.coli 数据在表 1 ，并在下面讨论。因此 S. cerevisiae 数据集被自愿排除在比较之外。这些统计数据在补充表 S2 。

Table 1.

Statistics of the simulated Pacific Biosciences long reads after correction with the different methods, as reported by LRCStats

Method	Original	CoLoRMap	HALC	HG-CoLoR	Jabba	LoRDEC	Nanocorr	NaS	Canu	Daccord	LoRMA
A.baylyi
错误率 (%)	17.8534	0.1023	0.0388	0.0310	0.0844	0.0374	0.5777	0.2961	8.5607	0.4967	21.7703
吞吐量 (Mbp)	71.9	62.6	63.7	64.6	61.6	61.1	64.4	48.2	59.1	64.7	0.2
删除 (%)	3.8909	0.0193	0.0134	0.0121	0.0629	0.0092	0.1000	0.0360	1.5961	0.0955	2.1088
插入 (%)	13.9605	0.0286	0.0184	0.0178	0.0233	0.0200	0.4967	0.1916	6.9289	0.4481	12.8142
换人（%）	0.7186	0.0700	0.0130	0.0059	0.0111	0.0137	0.0675	0.1138	0.7860	0.0290	8.9844
拆分读取 (%)	不适用	13.85	3.61	0.01	2.97	15.19	0	0	0	1.33	42.18
运行	不适用	57 分钟	22 分钟	47 分钟	2 分钟	6 分钟	2 小时 52 分钟	24 小时 24 分钟	10 分钟	20 分钟	5分钟
大肠杆菌
错误率 (%)	17.9267	0.1036	0.0601	0.0596	0.0462	0.0669	0.3983	0.1935	8.8525	0.4498	32.9756
吞吐量 (Mbp)	93.0	78.4	81.3	83.4	78.0	78.5	83.3	62.4	75.7	83.8	0.3
删除 (%)	3.9091	0.0198	0.0215	0.0280	0.0389	0.0133	0.0680	0.0257	1.6533	0.0866	2.8487
插入 (%)	14.0186	0.0303	0.0269	0.0347	0.0109	0.0397	0.3360	0.1385	7.1604	0.4019	19.3079
换人（%）	0.7215	0.0693	0.0209	0.0063	0.0036	0.0240	0.0549	0.0602	0.7997	0.0306	13.2367
拆分读取 (%)	不适用	13.13	9.62	0.03	4.57	26.12	0	0	0	0.23	42.58
运行	不适用	1 小时 25 分钟	24 分钟	45 分钟	2 分钟	8 分钟	3 小时 17 分钟	28 小时 48 分钟	12 分钟	27 分钟	5分钟

注意：突出显示每个统计数据的最佳结果。

LoRMA 的运行时很有竞争力，但它显然是所有工具中表现最差的，并且在两个数据集上显示出最小的吞吐量和最高的错误率。实际上，使用 LoRMA 纠错后的错误率甚至高于原始错误率。尽管长读的覆盖深度较低，但其他两个自我校正工具确实设法降低了错误率，尽管 Canu 在两个数据集上仅将错误率降低到 8% 左右，并且会从更高的覆盖率中受益。 Canu 和 Daccord 还设法产生了高通量的校正碱基，Daccord 甚至在该指标上优于所有其他方法。除 NaS 外，所有混合方法的吞吐量都相当。这种较低的吞吐量并不意味着 NaS 校正的碱基比其他工具少，而是来自这样一个事实，即 LRCStats 报告的吞吐量是根据原始读数和校正读数之间的比对计算得出的。由于 NaS 将短读取组装成重叠群并将它们用作校正的长读取，这表明产生的长读取的不可忽略量不会与它们最初来自的长读取重新对齐，而是在基因组的另一个区域对齐. 所有混合方法也令人满意地降低了错误率，始终低于 1%，最高来自 A.baylyi 数据集用 Nanocorr 校正，只有 0.57%。当更具体地查看不同类型的错误时，似乎 Jabba 在删除方面遇到了最大的困难，在 CoLoRMap 方面面临着替换，在所有其他工具方面都面临着插入方面的困难，这些都是分析数据集上最常见的错误。还值得注意的是，CoLoRMap、HALC 和 LoRDEC 确实拆分了重要比例的长读，这意味着这些长读的大区域没有被纠正，因此，参考基因组的大区域有很高的机会被发现. 在查看运行时，Jabba 总是表现最好，而 Nanocorr 尤其是 NaS 显示出令人望而却步的运行时。 HG-CoLoR 没有达到这样的运行时间，尽管它比除 CoLoRMap 之外的所有其他工具都慢。 HG-CoLoR 的错误率也低于所有其他工具，除了 大肠杆菌 数据集上的 Jabba。在这两个数据集上，HG-CoLoR 的吞吐量也高于除 Daccord 之外的所有其他工具。最后，它产生的拆分长读取的比例小于所有其他可以拆分读取的工具。因此，HG-CoLoR 在这两个模拟数据集上提供了运行时间和结果质量之间的最佳平衡。

4.3 真实数据对比

根据模拟数据的结果，我们选择将 LoRMA 从与真实数据的比较中排除，因为它的总体结果并不令人满意。我们还选择排除 CoLoRMap、HALC 和 LoRDEC，因为它们产生的拆分长读取占重要比例。最后，由于它在真实数据集上的所有混合工具中表现最差，我们也排除了 Nanocorr。中给出并评论了从该比较中排除的工具的结果补充表 S3 和 S4 。

我们用两种不同的方法评估不同工具的准确性。首先，我们分析长读取与参考基因组的对齐情况，其次，我们研究可以从校正的长读取中生成的组装质量。

4.3.1 基于对齐的比较

长读在校正之前与 Last ( Kielbasa et al. , 2011 ) 对齐，因为它更好地处理原始长读。然后运行不同的校正工具，将获得的校正长读数与 BWA mem ( Li and Durbin, 2010 ) 对齐，因为它们具有很高的准确性。结果在表2 并在下面讨论。

表 2。

使用不同方法校正前后真实牛津纳米孔长读数的统计数据

方法	原来的	HG-彩色	贾巴	在里面	唱歌	好的
A.baylyi
读取次数	89 011	25 278	16 618	24 063	8122	19 623
拆分读取 (%)	不适用	1.01	4.90	0	5.47	53.02
平均长度	4284	11 157	10 260	8840	9345	3244
碱基数 (Mbp)	381	285	179	213	81	175
平均同一性 (%)	70.09	99.75	99.40	99.82	97.79	91.92
基因组覆盖率 (%)	100	100	99.82	100	99.79	100
运行	不适用	1 小时 56 分钟	2 分钟	94 小时 18 分钟	32 分钟	45 分钟
大肠杆菌
读取次数	22 270	21 970	21 005	21 818	17 154	17 478
拆分读取 (%)	不适用	0.07	4.98	0	0.38	34.40
平均长度	5999	6093	5797	7926	7080	4495
碱基数 (Mbp)	134	134	128	173	122	119
平均同一性 (%)	79.46	99.84	99.81	99.86	96.23	98.51
基因组覆盖率 (%)	100	100	99.43	100	99.99	99.99
运行	不适用	1 小时 05 分钟	3 分钟	72 小时 02 分钟	36 分钟	30分钟
酿酒酵母
读取次数	205 923	72 228	33 484	71 793	—	—
拆分读取 (%)	不适用	5.13	11.47	0	—	—
平均长度	5698	6724	6455	5938	—	—
碱基数 (Mbp)	1173	512	243	426	—	—
平均同一性 (%)	55.49	99.10	99.54	99.59	—	—
基因组覆盖率 (%)	99.90	99.40	93.32	98.70	—	—
运行	不适用	8 小时 36 分钟	12 分钟	>16 天	—	—
秀丽隐杆线虫
读取次数	363 500	278 614	219 840	—	340 826	—
拆分读取 (%)	不适用	8.85	20.46	—	0	—
平均长度	5524	5127	3905	—	5408	—
碱基数 (Mbp)	2008	1567	1058	—	1843	—
平均同一性 (%)	71.07	98.93	99.85	—	85.63	—
基因组覆盖率 (%)	99.99	99.95	95.40	—	99.89	—
运行	不适用	80 小时 34 分钟	58 分钟	—	16 小时 38 分钟	—

注意：reads 列的数量是指已纠正的长读取数量，而不是输出片段的数量。上 NaS 的精确运行时间被省略， 酿酒酵母 因为在 16 天内没有计算结果，因此停止了执行。该数据集的 NaS 校正读数来自 Genoscope 网站。给定方法的省略结果意味着该方法未运行或无法纠正长读数。突出显示每个统计数据的最佳结果。

Jabba 显然在运行时表现最好，比所有其他工具高出几个数量级。它还产生了与高同一性对齐的更正长读数。然而，尽管高度准确，但这些校正过的长读数并未完全覆盖任何参考基因组。由于 Jabba 使用固定顺序的 de Bruijn 图，因此这些未解决的区域可能来自产生的拆分长读取的重要比例，因此面临局部覆盖率下降的困难。根据作者的建议，使用 Karect 而不是 QuorUM 对短读取进行预处理，并没有显示出任何显着的改进（参见补充表 S5 ）。

除了 Jabba，这两种自校正工具在运行时间方面优于其他两种混合校正方法。上产生的错误高达 8% A.baylyi ，尽管覆盖率为 106 倍，但使用 Daccord 校正的长读取的平均长度也小于原始长读取的平均长度，这是由于拆分长读取的比例很高。在 S. cerevisiae 数据集上，这两种工具都无法进行校正。 Canu 因错误而停止，这是由于长读取的高错误率引起的，而 Daccord 由于内存限制而失败，即使在我们可以访问 256 GB RAM 的具有大量资源的集群上运行时也是如此。这强调了一个事实，即混合纠错仍然更有效地纠正高噪声的长读取，即使来自中等大小的基因组。

因此，只有 NaS 和 HG-CoLoR 能够产生校正的长读数，这些长读数覆盖了具有高同一性的整个参考基因组，除了 酿酒酵母 ，因为原始长读数和短读数都没有覆盖整个基因组。然而，HG-CoLoR 输出更正的长读数，从而更好地覆盖参考基因组。对于所有数据集，用 NaS 校正的长读数比用 HG-CoLoR 校正的长读数具有稍高的同一性。然而，尽管在快速模式下运行，NaS 在所有数据集上都比 HG-CoLoR 慢了几个数量级。

在较大的 C. elegans 数据集上，NaS 由于其运行时间太长而没有运行，并且由于内存限制，Daccord 再一次无法执行校正，即使在上述集群上也是如此。 Canu 产生了最大数量的校正长读数。他们很好地覆盖了参考基因组，但仍然显示出很高的错误率，接近 15%。正如在模拟数据的实验中所提到的，如此大比例的错误是由于使用的 long reads 的覆盖率较弱，而更好的覆盖率可以大大改善 Canu 的结果，如在真实的 A. baylyi 和 E.大肠杆菌 数据集。 HG-CoLoR 尽管速度较慢，但仍能产生正确的长读数，显示出高同一性，并很好地覆盖了参考基因组。此外，值得注意的是，HG-CoLoR 的内存峰值仅为 10GB，即使在合理的设置下也能够扩展到大型基因组。

因此，尽管它的运行时间比自校正方法大，并且与 NaS 相比，它在校正的长读取的比对同一性方面略有劣势，但 HG-CoLoR 在运行时间和结果质量之间表现出最佳平衡，并且被证明是唯一能够有效扩展到真核基因组的工具。

4.3.2 基于汇编的比较

使用不同工具进行校正后，使用 Canu 组装长读数，无需执行校正和修剪步骤。我们只使用 Canu 而不是任何其他可能是混合的组装工具，因为组装不被认为是对校正的关键评估，也不是最重要的结果。以下参数用于所有程序集： OvlMerSize 、 –assemble = 17、OvlMerDistinct 0.9925 = 和 OvlMerTotal = 0.9925 。组大小参数独立设置为每个参考基因组的确切碱基数。所有其他参数都设置为其默认值。我们承认 Canu 最初的设计目的是将一整套未校正的读取作为输入，并同时执行校正和组装。然而，在实践中，这两个任务是独立的。长读数首先被纠正，然后只有可以纠正的长读数被组装。因此，直接将 Canu 可以纠正的长读取提供给其组装模块，而不是提供整组读取并执行所有任务，不会影响组装结果。使用 MUMmer 将程序集与参考基因组进行比较（库尔茨等人。 , 2004 年）。结果在表 3 并在下面讨论。

表3。

从校正后的真实牛津纳米孔长读数生成的组件的统计数据

方法	HG-彩色	贾巴	在里面	唱歌	好的
A.baylyi
长读覆盖率	79×	50×	59×	22×	49×
重叠群的数量	1	14	1	3	1
NG50	3 634 461	216 679	3 629 508	2 887 573	3 520 381
基因组覆盖率 (%)	99.99	89.03	100	99.39	100
身份（％）	99.94	99.94	99.99	97.04	97.06
大肠杆菌
长读覆盖率	29×	28×	37 ×	26×	26×
重叠群的数量	1	41	1	3	2
NG50	4 659 731	138 730	4 635 116	3 155 369	4 558 944
基因组覆盖率 (%)	100	95.81	99.90	99.82	100
身份（％）	99.99	99.99	99.99	97.23	97.84
酿酒酵母
长读覆盖率	41 ×	20×	34×	—	—
重叠群的数量	67	138	119	—	—
NG50	297 575	47 164	146 459	—	—
基因组覆盖率 (%)	97.57	68.67	97.44	—	—
身份（％）	99.92	99.99	99.95	—	—
秀丽隐杆线虫
长读覆盖率	15×	10×	—	18×	—
重叠群的数量	352	1 369	—	1 049	—
NG50	458 250	—	—	116 510	—
基因组覆盖率 (%)	98.41	53.86	—	88.75	—
身份（％）	99.86	99.97	—	95.60	—

注意：报告的身份代表一对一的对齐方式。省略 NG50 大小意味着组装没有达到参考基因组大小的一半。突出显示每个统计数据的最佳结果。

与我们在表 2 中观察到的一致，使用 Jabba 校正的长读数无法覆盖整个参考基因组这一事实导致高度碎片化的组装，无法解析参考基因组的大区域。结果，尽管它们的准确性很高，但这些经过校正的长读数产生了最碎片化的组装，覆盖了最少的参考基因组，并且显示了最小的 NG50 大小。特别是，在 秀丽隐杆线虫 数据集上，组装没有达到参考基因组大小的一半。然而，MUMmer 报道了略大于 50% 的基因组覆盖率，这是由重复区域中非常小的重叠群（几乎没有长读长的大小）的双重映射引起的。

如前所述，与混合方法相比，这两种自校正工具产生的长读取覆盖深度较弱，显示出更高的错误率。尽管如此，这些长读数仍然组装成少量的重叠群，在 A.baylyi 和 E.coli 数据集上都显示出高 NG50 大小。然而，由于长读取的高错误率，与参考基因组相比，这些组件显示出最低的同一性。此外，在这两个数据集上，Daccord 在连续性、NG50 和组装的基因组覆盖方面优于 Canu。在 A.baylyi 数据集上，使用 Daccord 校正的长读取甚至组装成单个重叠群。

在 A.baylyi 和 E.coli 数据集上，NaS 和 HG-CoLoR 产生了正确的长读数，这些读数组装成一个单一的重叠群。因此，NG50 尺寸非常相似，两个组件仅显示出两个细微差别。一小部分 大肠杆菌 由 NaS 的校正读数产生的组装没有覆盖其次，在 A.baylyi 数据集上，与参考相比，由 HG-CoLoR 的更正读数生成的程序集显示出略低的身份。然而，在 酿酒酵母 数据集上，HG-CoLoR 在连续性、NG50 大小和组件的基因组覆盖率方面优于 NaS。

在较大的 秀丽隐杆线虫 数据集上，尽管用 Canu 校正的长读数很好地覆盖了参考基因组，但它们生成的组装未能解析 11% 的参考基因组。如前所述，与参考基因组相比，同一性也很低。因此，只有用 HG-CoLoR 校正的长读数才能导致组装显示出高同一性，并很好地覆盖参考基因组。重叠群的数量和 NG50 的大小也比其他组装的更令人满意，尽管它们与实际基因组的那些仍然相去甚远。

4.4 HG-CoLoR 结果的细节

我们进一步关注 HG-CoLoR 结果，并更深入地分析执行校正的方式，通过观察已校正的长读取的碱基来源，无论是模拟数据集还是真实数据集。我们还评论了短读的长度如何影响结果的质量。

4.4.1 碱基的起源

中提供了一个描述，表明这些碱基是来自种子、图遍历还是原始长读取表 4 。看起来，在三个模拟数据集上，错误率约为 17%，大多数已校正的长读取的碱基来自种子。由于与其他数据集相比，这些长读的错误率相对较低，短读确实可以更容易对齐，覆盖长读的更大区域，因此需要更少的图遍历。只剩下很少的未校正碱基，在 酿酒酵母 数据集上，最大比例仅略大于 1%。这些观察结果也适用于真实的 大肠杆菌 数据集，其错误率略高于 20%。

表 4。

来自种子、图遍历和原始长读的碱基比例，在用 HG-CoLoR 校正的长读中

数据集	种子基地 (%)	图基数 (%)	原料碱 (%)
A.baylyi （模拟）	68.95	30.97	0.08
大肠杆菌（模拟）	67.01	32.74	0.25
酿酒酵母（模拟）	63.97	34.82	1.21
A.baylyi （真实）	36.37	60.53	3.10
大肠杆菌（真实）	62.61	37.04	0.35
酿酒酵母（真实的）	24.41	60.22	15.37
C. elegans (真实的)	48.94	44.23	6.83

在其他三个显示出更高错误率的真实数据集上，更多的基础来自图的遍历。在 C. elegans 数据集上，错误率为 29%，来自种子的碱基几乎与来自图遍历的碱基一样多。随着错误率的上升，这种趋势进一步加剧。在错误率高于 44% 的 S. cerevisiae 数据集上，来自图遍历的碱基数量是来自种子的两倍多。无论如何，对于所有数据集，来自图遍历的不可忽略比例的碱基强调了图的有用性。

未校正碱基的比例也随着错误率和基因组复杂性的增加而上升。上可以观察到最大比例的此类碱基 酿酒酵母 。然而，在 C. elegans 错误率低于 A.baylyi ，几乎 7% 的碱基无法纠正。这可能来自基因组的复杂性，以及重复区域因此更难解析的事实，即使图中使用了较高的最大阶数。

4.4.2 短读长的影响

我们实验中使用的所有短读长都是 250-300 bp，尽管绝大多数 Illumina 读长都较短，并且仍在 100-150 bp 范围内。然而，HG-CoLoR 并不要求短读取要那么长才能执行有效的纠错。 HG-CoLoR 在模拟和真实数据集上的结果，当使用 125 bp 短读取的 50 倍覆盖率时，在补充表 S6-S8 ，并表明 HG-CoLoR 在如此短的读取下实际上表现得更好。

4.5 变阶 de Bruijn 图的影响

为了进一步验证使用可变阶 de Bruijn 图而不是经典图的有用性，我们在真实的 酿酒酵母 数据集上进行了另一项实验。我们使用与之前的实验相同的参数，除了图的最小阶数，我们将其设置为与最大阶数相同的值，100，以便仅遍历经典的 de Bruijn 图。中给出了该实验的结果，以及与使用可变阶 de Bruijn 图时获得的结果的比较表 5 。

表 5。

可变阶 de Bruijn 图对真实 酿酒酵母 数据集

图形	变阶	古典
更正
读取次数	72 228	72 228
拆分读取 (%)	5.13	6.89
平均长度	6724	6525
碱基数 (Mbp)	512	507
平均同一性 (%)	99.10	99.06
基因组覆盖率 (%)	99.40	99.40
运行	8 小时 36 分钟	6 小时 13 分钟
集会
长读覆盖率	41×	41×
重叠群的数量	67	107
NG50	297 575	153 845
基因组覆盖率 (%)	97.57	96.57
身份（％）	99.92	99.92

这些结果表明，使用可变阶 de Bruijn 图确实允许执行更有效的纠错，特别是在对应于短读取覆盖不足的区域的子图中。当使用经典的 de Bruijn 图时，这些子图会导致更多的未校正区域，从而导致更多的分裂和更短的长读取，实际上可以使用可变阶 de Bruijn 图正确校正，而不会损失任何质量。

尽管图的类型不会影响长读的基因组覆盖率，但使用经典的 de Bruijn 图时，组装的质量会大大降低。事实上，用经典的 de Bruijn 图校正长读数会导致由 40 个以上的 contigs 组成的组装，达到 NG50 小两倍，并且覆盖参考基因组比使用可变的长读数校正后生成的组装少 1%顺序 de Bruijn 图。

5 结论

我们描述了 HG-CoLoR，这是一种用于长读取纠错的新混合方法，它将最先进的两种方法结合到种子和扩展策略中。短读首先与长读对齐，以便发现种子。然后将这些种子用作可变顺序 de Bruijn 图上的锚点，该图由短读取构建，并通过遍历图链接在一起，以纠正长读取的未覆盖区域。

我们的实验表明，与最先进的混合和非混合纠错工具相比，HG-CoLoR 在运行时间和结果质量之间提供了最佳权衡，无论是在纠错本身的质量方面，以及从更正的长读取生成的程序集的质量。进一步的实验还表明，我们的方法是唯一能够有效扩展到真核基因组的方法。

这种方法的发展和我们的实验强调了这样一个事实，尽管自我校正方法已经有用，但仍然不能令人满意地适用于高噪声长读取。事实上，他们根本无法对使用早期化学方法测序的长读取进行纠错，这显示出非常高的错误率。它们也不能扩展到真核基因组，要么完全无法进行校正，要么几乎没有降低错误率。因此，无论是在真核基因组的情况下，还是在错误率非常高的情况下，混合方法在纠正长读数方面仍然更有效，因为重新测序并不总是负担得起的，而且即使是最近的数据仍然经常显示出高错误率。

我们的实验还强调了使用可变阶 de Bruijn 图而不是经典图的有用性。实际上，可变阶 de Bruijn 图允许有效地纠正长读取，特别是在短读取覆盖率不足的区域，以及经典 de Bruijn 图由于缺少边缘而无法执行校正的区域。因此，可变阶 de Bruijn 图的使用允许更正长读取的更大区域，而不会损失任何质量。这不仅改善了纠错本身，而且有利于下游组件的连续性和整体质量。

作为进一步的工作，我们计划专注于 PgSA 的新实现，因为当前不支持索引的并行查询，因此迫使 HG-CoLoR 使用互斥锁。摆脱这种需求将减少该方法的运行时间。另一个允许查询一组具有可变长度字符串的 k -mer 并支持并行查询的索引结构也可以替代 PgSA。另一个方向是尝试其他对齐器用于短读取到长读取的对齐步骤，以便可能更快地发现种子，或者纠正更多的长读取。