VeChat：使用变化图纠正长读取中的错误
摘要
纠错是长读序列数据分析的第一步。目前的标准是使用共识序列作为模板。然而，在混合样本中，例如亚基因组或更高倍性的生物体中，一致性诱导的偏见可以掩盖影响较低频率单倍型的真实变体，因为它们被误认为错误。
这里介绍的新颖之处是使用基于图形，而不是基于序列的共识作为识别错误的模板。其优点是，基于图形的参考系统也能捕捉低频率的变化，因此不会错误地将其掩盖为错误。我们提出了VeChat，作为实现这一想法的一种新方法：VeChat根据变异图区分错误和单倍型特定的真变异，变异图反映了泛基因组参考系统的一种流行数据结构。在从原始输入读取中初步构建一个特殊变化图时，通过一种基于频繁项集挖掘原理的迭代过程，从该图中删除由错误导致的节点和边。终止时，图形只包含反映真实连续现象的节点和边。原始读数的最终重新校准表明需要在何处以及如何更正读数。大量的基准测试实验表明，通过VeChat校正的PacBio和ONT读数包含4到15个错误，或者分别比校正最先进方法时的错误少2到10倍。VeChat是在一个易于使用的开源工具中实现的，可在https://github.com/HaploKit/vechat.上公开获取。
关键词 纠错；长篇大论；单倍型；多倍体；超基因组
背景
第三代测序（TGS），如单分子实时测序（Pacific Biosciences，简称PacBio）或纳米孔测序（Oxford nanopore Technologies，简称ONT），在过去几年中迅速出现。除了TGS读取的明显原因，从数量级上看，TGS读取的长度往往更长，从几kbp到几Mbp不等（Logsdon等人，2020年），而下一代测序读取通常只跨越几百个碱基对之外，TGS相对便宜的事实增加了它的受欢迎程度。此外，TGS作为其方案的一部分，避免了聚合酶链反应（PCR），从而防止了相关的偏见。多亏了这些优势，TGS一直 能够在各个应用领域做出决定性的贡献。突出的例子包括单倍型阶段化（Schrinner等人，2020年）、基因组组装（Jain等人，2018b；Ruan和Li，2020年；Shafin等人，2020年；Kolmogorov等人，2020年；Miga等人，2020年）和（复杂）变体调用（Edge和Bansal，2019年；Thibodeau等人，2020年；Fujimoto等人，2021年）。
然而，另一方面，TGS的缺点是读操作受到影响的错误率显著升高。例如，作为目前最具代表性的TGS读取示例，PacBio CLR和ONT读取包含5%到15%的错误（Logsdon等人，2020年）。这与NGS短读形成明显对比，后者的错误率通常不超过1%。大多数影响长读的错误都是由插入和删除组成的，这一事实增加了困难，因为它妨碍了原则的应用和工具的直接调整，以纠正短读中的错误。无论如何，在大多数相关应用中，直接使用原始TGS读取几乎是不可能的。因此，需要新的方法和工具。
由于纠正TGS读取中的错误对于声音分析来说是必要的，并且纠正TGS读取中的错误是必须的，因此同时提出了各种TGS读取错误纠正方法。相应的方法范围可分为两大类：混合校正和自校正。混合校正解决了在过程中集成短读的问题，而自校正则寻求在没有辅助数据的情况下纠正错误。
显然，混合修正总体上反映了一种合理合理的方法（参见（Hackl等人，2014年；Salmela和竞争对手，2014年；Firtina等人，2018年；Morisse等人，2018年）。然而，混合修正存在一些实际问题。首先，虽然长读可以跨越重复区域，但短读不能；这会在分配短读到长读的过程中产生歧义（反之亦然），从而导致校正质量出现偏差。其次，短读重新引入PCR诱导的偏见。例如，由于序列内容（例如GC内容），短读没有充分覆盖基因组的某些区域。这阻碍了这些领域的错误纠正。最后但并非最不重要的一点是，采用几种不同方案对基因组样本进行测序可能是完全不可能的，这从一开始就阻碍了混合纠错的应用。作为第二类方法，
自我纠正不存在上述任何问题。然而，由于缺乏外部（例如，基于短阅读的）帮助，自我纠正面临着其他有条不紊的挑战。在人们能够从自我纠错的巨大实际优势中获益之前，克服这些挑战是关键。
就先前的相关工作而言，自我修正可进一步分为三个子类，每一个子类都有特定的算法策略和方法基础。
第一种，也是三种类型中最常见的一种，是基于多序列比对（MSA）。对于之前主要依赖计算MSA的方法和工具，请参见Racon（Vaser等人，2017年）、组装程序Canu的纠错模块（Koren等人，2017年）和FLAS（Bao等人，2019年）。
第二类原则性方法依赖于德布鲁恩图（DBG）。相应的工具在对校正过程至关重要的某个时刻使用DBG。需要考虑的普遍工具是Daccord（Tischler和Myers，2017），它基于提高局部DBG，其中局部DBG指的是反映相对较小的基因组片段的DBG。
第三类自校正方法收集在过程中同时使用MSA和DBG的方法。这种有条不紊的组合方法试图平衡两个概念的优缺点，一方面是MSA，另一方面是DBG。成功地将MSA与DBG结合使用的主要工具是LoRMA（Salmela等人，2017年）和同意（Morisse等人，2021年）。
统一所有这些方法的共同点是提出共识顺序，该顺序反映了所观察到的读数的摘要，并在纠错过程中作为指导方针。然而，基于序列的模板无法捕捉歧义，whcih解释了校正过程中的相应偏差：对于每个多态性位点，需要决定一个可能的等位基因，同时丢弃所有其他等位基因。因此，这些方法往往掩盖了混合样本（亚基因组、癌症基因组）或多倍体基因组中几乎没有覆盖的单倍型或菌株的变异。由于受影响的单倍型几乎消失，下游分析对它们仍然是盲目的
为了解决这个问题，我们开发了VeChat，这是一种自我校正方法，可以对长读取执行单体型识别错误校正。从更大的角度来看，VeChat考虑了所有可能的单倍型的全谱，这些单倍型可能会在纠错期间影响样本，而不仅仅是在纠错之后。这反映了对于大于两倍体的倍体的新颖性，因为早期的方法只处理二倍体情况（Luo等人，2021年）。从系统的角度来看，VeChat的新颖之处在于将变化图作为基本数据结构集成到错误纠正过程中。据我们所知，VeChat是实现这一目标的第一种方法。
我们对VeChat进行了测试，并在反映当前感兴趣的各种设置的数据集上将其与当前技术水平进行了广泛比较。在模拟和真实数据上的基准测试实验表明，我们的方法在所有类型的常见测序错误中都达到了最佳性能
结果
我们设计并实现了VeChat（基于[V]变异图的[ha]plo[t]型错误校正），这是一种新的单体型识别长读自错误校正方法。向量背后的关键概念是变化图。与当前的方法通常以单一一致序列为中心不同，变异图能够代表多个进化或环境一致的基因组的遗传多样性。这使得在纠错过程中，对于更高、已知或未知倍性的样本，也能保留单倍型特异性变异。在本节中，我们首先对VeChat的工作流程进行高级描述。然后，我们评估了VeChat在模拟和真实数据上的性能，并与最先进的方法进行了比较。
工作流程
请参见图1了解VeChat的工作流程。VeChat由两个循环组成。第一个周期产生预校正读取，第二个周期从预校正读取生成最终的校正读取。每个循环分6步进行。虽然这两个周期在这6个步骤上大体一致，但在影响步骤1和4的关键细节方面，它们存在分歧。
在第一个周期中，步骤1计算基于最小值的所有与所有重叠，为此我们使用Minimap2（Li，2018）。因为使用Minimap2可以避免计算基准面路线的需要，所以这一阶段进展迅速，无需进一步努力。
步骤2-6反映了图1中纠错程序的技术核心。在步骤2中，选择一个目标读取作为要纠正错误的读取，并计算一个由所有与之重叠的读取组成的读取对齐堆。随后，在步骤3中，将读取对齐堆划分为小段，其中在步骤3中，每个段产生堆的窗口状部分；在特定窗口中读取的目标部分进一步称为“目标子读取”。
随后，在步骤4中，在每个窗口中执行目标子读取的错误校正。第4步在系统上更加复杂，因为它捕获了新颖的、基于变化图的方法；关于第一个循环中使用的特定步骤版本的详细说明，请参见图2（a）。第4步涉及构建变化图，并以迭代方式修剪该图（“图2（a）中的“图修剪”和“图重新修剪”）。剪枝反映了去除虚假节点和边，并根据频繁项集模型对节点和边进行建模，该模型涉及读取覆盖率、排序错误和读取中字符的共现；详情见小节。第一个循环以串联一个目标读取的不同“目标子读取”结束，这将导致以完整原始长度进行预校正读取。然后，这些预先校正的读数将作为第二个周期的输入。
虽然（绝大多数）错误已在第一个周期内得到纠正，但少数被视为统计异常值的错误在第一个周期的迭代图修剪过程中未得到纠正，因此拒绝了纠正。第二个周期应该能发现这样的异常值。由于第一个周期比第二个周期在统计上更复杂，所以它不包括第二个周期；图1中的蓝色元素指出了第二个循环的不同路径。总的来说，如上所述，第二个循环与步骤2、3、5和6中的第一个循环相同。然而，在步骤1中，除了计算所有与所有重叠之外，还计算了基础水平对齐，这使得能够识别单倍型的读取重叠过滤成为可能。第4步，如图2（b）所示，然后以不同的方式进行，因为图形（重新）修剪不再是周期的一部分。 而不是迭代的重新修剪，这并没有导致删除我们想要删除的错误 ，在第二个周期中，通过使用动态规划算法（Lee，2003），仅从构建的变化图中导出一致性序列（在图2（b）中显示为粗黄色箭头）。在步骤5中连接一致序列后，在步骤6中连接所有目标读取生成Vechat的最终输出。参考方法中的第一个和第二个循环，详细描述了步骤1-6。
图1。VeChat的工作流程。周期1和周期2的输入和输出分别标记为紫色和蓝色。循环1和循环2都有步骤1-6，除了步骤1和4中的一些差异。目标读数以橙色突出显示。红色叉表示读取中的顺序错误。
图2。一个目标子读取的错误更正。（a） 是周期1中一个目标子读取的纠错图。假设二倍体情况（橙色路径代表最佳对齐路径，而绿色路径代表另一个真正的单倍型），说明了目标读取r的错误纠正过程。”“图形修剪”和“图形重新修剪”指的是核心错误纠正程序。这些程序依赖于一个变化图，该变化图是由读对齐堆的段构成的，读对齐堆是目标读和重叠读的多次对齐的结果，见图1。在图形修剪和重新修剪过程中，由排序错误引起的伪边（虚线箭头）将从变化图中删除。粉红色的元素表示这些过程是重复的。（b） 是周期2中一个目标子读取的错误纠正图。粗体橙色路径表示共识序列。在第1周期（见（a）中）中预校正靶亚读物中的假核苷酸“G”标记为红色，并在第2周期中进一步校正。
数据集
模拟测序数据
我们利用最新工具PBSIM2（Ono等人，2021年）分别使用内置P6C4和R103基于模型的模拟曲线模拟PacBio CLR和Oxford Nanopore读数。由于VeChat的主要应用场景是纠正来自多个基因组（如多倍体基因组和亚基因组）的长读取测序数据，因此我们模拟了这两种情况下的各种数据集。
多倍体基因组数据。我们构建了假二倍体（倍性=2，ANI:98%）、三倍体（倍性=3，ANI:96%）⇠ 98%）和四倍体（倍性=4，ANI:96%⇠ 99%）混合大肠杆菌（E.coli）菌株的基因组；请注意，平均核苷酸同一性（ANI）的定义是为了测量基因组序列的相似性，例如，FastANI（Jain等人，2018a）可以报告这一点。大肠杆菌的所有基因组序列均从NCBI数据库下载（参考基因组详情见补充表S1）。从单倍型（即菌株）独立模拟读取，并在世代混合后形成相应的多倍体基因组数据集（倍性=2,3,4）。我们模拟了这些数据集的PacBio CLR和Nanopore读取，每个单倍型的平均测序覆盖率为30倍，平均测序错误率为10%。
元基因组数据。此外，我们使用CAMISIM（Fritz et al.，2019）模拟了两个不同复杂程度的宏基因组数据集（PacBio CLR读取）。在这里，我们使用PBSIM2来模拟PacBio CLR读取，而不是CAMISIM中的内置模拟器。低复杂性数据集由10个物种（20个菌株）组成，而高复杂性数据集由30个物种（100个菌株）组成。两个数据集中使用的基因组均源自（Quince等人，2017年）（详情见补充表S1）。对于这两个数据集，菌株的平均测序覆盖率约为30倍，平均测序错误率为10%。菌株的相对丰度分别为1.9%至10.6%和0.28%至3.3%。
真实测序数据（模拟群体）
酵母伪二倍体基因组数据。我们通过混合两种ANI为98.4%的酵母菌株（N44、CBS432）构建了一个伪二倍体基因组，这两种酵母菌株来自酵母群体参考小组（见补充表S1）。相应的真实PacBio CLR读数从欧洲核苷酸档案（ENA）项目PRJEB7245下下载，我们只对每个菌株30倍测序覆盖率的长读数进行二次抽样，以进行进一步分析。
NWC元基因组数据。我们下载了两个来自天然乳清发酵剂培养物（NWC）的真实宏基因组数据集（PacBio CLR reads）（Somerville et al.，2019），并将其混合在一起（见补充表S1），然后对20%的读数进行二次抽样，从而获得了一个低复杂度的宏基因组数据集，其中包含3个物种（6个菌株）。
微生物10倍元基因组数据。我们从一个由PacBio Sequel系统Chemistry v3测序的10倍多重数据集中下载了原始长读测序数据和相应的参考基因组。https://downloads.pacbcloud.com/public/dataset/microbial_ multiplex_dataset_release_SMRT_Link_v6.0.0_with_Express_2.0/）。然后，我们随机抽取了10%的读数，从而获得了一个平均测序覆盖率约为40倍的模拟宏基因组数据集，其中包含7个物种（共9个菌株，ANI<98.5%，见补充表S1）。
基准测试：可选方法
为了实现公平和有意义的比较，我们考虑了所有流行的最先进的工具，这些工具执行TGS读取自校正。也就是说，这一选择包括Racon（v1.4.13）（Vaser等人，2017年）、同意（v2.2.2）（Morisse等人，2021年）、Canu（v2.1.1）（Koren等人，2017年）和Daccord（v0.0.18）（Tischler和Myers，2017年）。我们使用默认参数运行所有工具。
评估指标
基因组组装性能评估指标通常通过几种常用指标进行评估，如QUAST V5所述。1.0（Mikheenko等人，2018年）。具体解释请参见下文，并参见http://quast.sourceforge.net/docs/manual.html获取完整详细信息。请注意，标准QUAST标准也涵盖了纠错长读的原则性质；例如，低单倍型覆盖率反映出真实的变异被错误地识别为错误，而错误率和错误装配则反映出纠正程序忽略了错误，甚至是错误纠正导致的混淆读。和往常一样，在评估之前，从输出中过滤长度小于500bp的校正读数。请注意，我们运行QUAST时使用了选项——模糊用法一，适当地考虑到我们的数据集反映了混合样本（例如多倍体基因组或亚基因组）。
错误率（ER）。将获得的校正读数映射到参考单倍型序列时，错误率等于失配率和indel率之和。
单倍型覆盖率（HC）。单倍型覆盖率是校正读取覆盖的基本真值单倍型中对齐碱基的百分比，用于测量校正读取的完整性。
N50和NGA50。N50定义为该长度或更长的所有校正读取的集合至少覆盖给定序列的一半的长度。NGA50与N50相似，但只能在提供参考基因组时计算。NGA50只考虑对齐的区块（在错误组装事件中中断读取并修剪所有未对齐的核苷酸后），其定义为该长度或更长的所有对齐区块的总大小至少等于参考单倍型的一半的长度。N50和NGA50均用于评估校正读数的长度分布。请注意，这可能对更正的读取相对不感兴趣。尽管如此，我们还是显示了相应的结果，因为误差校正确实会对这些统计数据产生影响。
错误装配的数量（#错误装配）。校正读取中的错误组装事件表明，左右两侧序列与真正的单倍型对齐，间隙或重叠超过1kbp，或与不同的链对齐，甚至与不同的单倍型或菌株对齐。这里，我们报告给定序列数据中错误组装的总数。
基准测试结果
模拟了多倍体基因组数据集：PacBio。表1显示了从不同倍体（即2、3和4）的基因组中模拟PacBio CLR读取的误差校正基准测试结果。VeChat大约达到14- 30, 9 -26和4-11的错误率分别降低 在二倍体、三倍体和四倍体基因组。同时，它在其他方面保持了更好或可比的性能，如校正读取的数量、完整性（单倍型覆盖）、错误组装的数量和校正读取的长度（如N50/NGA50所示）。特别是，VeChat在不匹配率方面显著优于其他工具（4-69倍 比其他人低）。
模拟多倍体基因组数据集：ONT。表2显示了从不同倍性（即2、3和4）的基因组中读取模拟牛津纳米孔的误差校正基准结果。VeChat实现了大约10-20, 3 -9和2-5 降低的错误率在二倍体、三倍体和四倍体基因组，同时在所有其他方面保持更好或可比的性能。与PacBio reads一样，VeChat在失配率方面也表现出更好的性能：2-59倍比其他矫正工具，与indel rate相比，。
模拟元基因组数据集。表3显示了具有不同复杂性的宏基因组数据集的模拟PacBio CLR读取的纠错基准测试结果。VeChat实现了大约6-7和3-4倍 在低复杂度和高复杂度元基因组的错误率，同时在其他方面保持可比性能。特别是，VeChat在失配率方面明显优于其他工具6-12在低复杂度数据集上。
表1。各种多倍体基因组（倍性=2,3,4）模拟PacBio CLR读取的误差校正基准结果。每个单倍型的平均测序覆盖率为30倍，测序错误率为10%#Seq'表示已更正的读取次数。错误率等于失配率和失效率之和。结果按错误率升序排序。
表2。各种多倍体基因组（倍性=2,3,4）模拟牛津纳米孔读数的误差校正基准测试结果。每个单倍型的平均测序覆盖率为30倍，测序错误率为10%
表3。具有不同复杂性的宏基因组数据集的模拟PacBio CLR读取的纠错基准测试结果。菌株的平均测序覆盖率约为30倍，测序错误率为10%。
真实测序数据（模拟社区）。表4显示了真实测序数据的误差校正基准测试结果。表中的三个部分显示了酵母伪二倍体基因组数据集（模拟群体）的结果，NWC元基因组数据集（真实群体）的结果，以及微生物10倍元基因组数据集（模拟群体）的结果，作为表中行的第三部分？？。VeChat实现了大约2-4, 1.4 -7.8和3.3-5.6倍低错误率在 酵母、NWC和微生物10倍数据集，同时在其他方面保持可比性能。
不同的读取覆盖率
为了评估测序覆盖率如何影响纠错方法，我们将重点放在三倍体基因组上，如前所述，三倍体基因组由三个大肠杆菌菌株组成。我们模拟了不同测序覆盖率下的PacBio CLR读取，即每个单倍型10倍、20倍、30倍、40倍、50倍。补充表S2显示了误差修正的基准测试结果。总之，VeChat实现了大约2⇠47倍 降低 在所有数据集的错误率，同时在其他方面保持更好或可比的性能，如纠正读取的数量、完整性（HC）和纠正读取的长度。随着测序覆盖率的增加（从10倍增加到50倍），VeChat实现了更好的错误纠正（错误率从0.311%增加到0.017%），同时在其他方面保持了相当的性能。
不同的测序错误率
为了评估测序错误率对不同方法的影响，我们再次关注由上述三种大肠杆菌菌株组成的三倍体基因组。因此，我们模拟了不同测序错误率下的PacBio CLR读取，即5%、10%和15%。每个单倍型的平均测序覆盖率是30倍。
补充表S3显示了相应的误差校正基准测试结果。总的来说，VeChat实现了大约10-93, 9 -26和7-9倍 低错误率 在数据集 以5% 10% 15%的错误，同时在其他方面保持可比性能。在降低测序错误率（从15%降至5%）的同时，VeChat的错误纠正率无疑有所提高，错误率从0.091%降至0.009%
运行时和内存使用评估
VeChat的运行时主要由三个步骤组成：虽然读取重叠的计算（无基线对齐）非常快，但后续基于编辑距离的对齐（用于分割窗口）非常耗时。其次，驱动变异图构造的POA算法执行序列到图的对齐，其计算复杂度为O（N（2Np+1）|V |），其中N是要对齐的序列的长度，Np是图中前导的平均数，|V |是图中的顶点数（Lee et al.，2002）。第三，VeChat遵循一种迭代范式，例如在第二次迭代中执行读取重叠计算（与基准层对齐）和错误纠正（一致性生成）步骤；这也需要相当长的运行时间。
我们使用48个内核在x86 64 GNU/Linux机器上执行了所有基准测试分析。补充表S4-7中报告了不同方法的运行时和峰值内存使用评估.VeChat需要23-81个CPU小时和8-92个CPU小时 在模拟PacBio CLR和ONT从反映不同倍性（通常为2、3或4）的数据集读取时，，即1.1-6.2和0.6-6.9倍比其他方法慢。同时，它需要更高的峰值内存使用率（补充表S4、S5）。这是2.4-7.1倍慢于模拟宏基因组数据集上的其他方法，并达到更高的峰值内存使用率（补充表S6）。请注意，高度复杂的宏基因组数据集太大，Racon和Daccord无法处理，而我们的方法能够处理它。此外，VeChat为0.7⇠32倍慢 在实际数据集上，同时需要更高的峰值内存使用率（补充表S7）。
讨论
我们介绍了VeChat，作为一种对第三代测序（TGS）读取执行单体型感知错误纠正的方法。据我们所知，VeChat是第一种明确解决在校正过程中保留单倍型特异性变异的方法，它比以前的方法有所改进，尤其是在可以避免产生偏见的一致序列的情况下。结果证明了VeChat的优越性：在所有基准测试场景中，VeChat将错误率抑制至少2或3倍，如果不是，就像大多数场景一样，与领先竞争对手相比，将错误率抑制一到两个数量级。同时，VeChat保留了读取的单倍型身份，这意味着在使用VeChat进行校正后，所有读取都有助于覆盖它们产生的单倍型。对这些结果最明显的解释是，当试图从TGS读取中删除所有错误时，在纠错过程中捕获单倍型结构不仅是有益的，甚至是必要的。为了在纠错步骤中恰当地捕获单倍型特异性变体，我们直接从嘈杂的TGS读取构建变体图。请注意，从严重错误的读取直接构造变化图不是标准。事实上，乍一看，这甚至是违反直觉的，因为变异图的基本思想是由真正的单倍型特异性、足够长的序列补丁构建而成。在这里，序列补丁甚至包含15%的错误。
因此，生成的图形包含大量反映此类错误的节点和边。为了识别虚假节点和边缘，可以利用测序错误是随机分布的，而变体往往会在不同的读取中重复出现。特别是，链接虚假节点（与真实节点或其他虚假节点）的边往往很少被读取覆盖，因为与真实效果不同，读取不会共享错误。为了系统地将伪边识别为被太少读取覆盖的边，我们采用了两种来自频繁项集挖掘的度量标准，这是一种合理、有原则的方法。支持度衡量的是图中边的相对覆盖率，置信度衡量的是与它们共享的边相关的两个节点之间的关联；如果基本支持或关联太少，边缘和可能产生的孤立节点将被删除。
VeChat的一个特别效果是在错配率方面取得了显著的改善；在所有情况下，indel率的改善也很明显，但通常不一定是剧烈的。一种可能的解释是，替代事件比插入和删除事件更能主导生物体的进化过程，因此通常是菌株或单倍型的特征。VeChat似乎是第一个正确保存这些单核苷酸多态性（SNPs）的方法，因为单倍型的区别正是变异图的目的。无论如何，VeChat似乎可以防止由于产生一致序列而掩盖真实的变异。
至于未来的前景，长读测序技术的快速发展将导致测序错误率的降低。然而，由于当测序错误率下降时，VeChat的优势变得更加明显，因此在纠正错误率为5%或更低的长时间读取时，VeChat也将是一个优越的工具（见补充表S3）。
当然，未来的改进是可以想象的：尤其是VeChat需要的计算资源超过了其他方法。尤其是在大型数据集上，VeChat的运行时间更长，峰值内存使用率更高。然而，VeChat的计算效率还有改进的空间。一个关键点是，VeChat在某些地方使用现成的方法。虽然这些方法可以完成这项工作，但并不是它们提供的所有功能都是必需的，这相当于运行现成软件时的计算开销。特别是，在未来，诸如编辑基于距离的对齐或序列到图形的对齐等计算很有希望被更高效的例程所取代。
方法
步骤1：读取重叠计算
步骤1指使用（广泛流行的）基于最小值的长读取重叠计算工具Minimap2（Li，2018）计算输入读取的所有与所有重叠（第一周期：原始读取，第二周期：预校正读取）。在第一个循环中，仅执行种子链程序，而在第二个循环中，添加基准面对齐。因为这是非常快的，它可以轻松地管理我们需要处理的大量读取对。随后，根据合理的附加标准过滤不良重叠，包括删除长度不超过500 bp的重叠、自重叠或内部匹配。在这方面，我们遵循（Li，2016）中的算法5及其在（Marijon等人，2020）中的实现。此外，具有高错误率的重叠，即、 其中Lq和Lt分别是查询和目标读取的重叠长度，e是最大错误率阈值，也会被过滤掉。
在第二个循环期间，也会执行类似的程序（重叠计算和过滤），但用于预校正读取。与第一个周期不同的是，我们使用基准级对齐计算预校正的读取重叠，以便可以确定重叠的序列标识（重叠标识），然后使用另一个标准最小重叠标识（表示为 δ,δ= 0.99用于模拟测序数据和δ = 0.98实际测序数据）。值得注意的是，由于大多数测序错误已在第一个周期内得到纠正，因此来自不同单倍型的读取重叠的身份分布明显分散。因此，只需设置重叠标识≥δ，就可以直接从不同的单倍型中过滤重叠群
第2步：读取对齐堆生成
我们的工作流程然后选择一个目标读取r，并基于在第1步中计算的重叠，收集该重叠r的读取。目标读取r用作主干，对于读取r和另一个读取之间的每个重叠，执行基于编辑距离的快速对齐（Myers，1999；ˇSoˇsi'c和ˇSiki'c，2017），这将生成一个读取对齐堆。在步骤3中，只需将读取对齐堆拆分为小窗口，就可以使用基于编辑距离的对齐。与目标读取r重叠的读取的悬空端， 由图1中原始读取对齐桩中的水平虚线指示，在下文中不再进一步考虑。更多详情请参见（Vaser等人，2017年）
步骤3：窗口分割<\/br><\/br>
然后将读取对齐堆栈划分为几个长度相同的不重叠小窗口，目标读取作为参考：每个这样的窗口覆盖目标读取的500 bp。显然，分割读取对齐堆反映了一个简单的过程，因为目标读取与其重叠读取的成对对齐是读取对齐堆构造的基础，参见步骤2（Vaser等人，2017）。与一个特定窗口对应的目标读取r的部分进一步被称为“目标子读取”。这尤其意味着，除了最右边的窗口外，目标子读取的长度为500 bp，其中目标子读取的长度可以更短，将对齐桩分割为小长度窗口的原因是以下方面的计算负担大大减少：下一步4，作为我们方法的技术核心，涉及变化图，无论是在缩小原始问题的规模方面，还是在实现并行化方面，都可以从这种细分中获得巨大的利润<\/br>
第4步：目标子读取的错误纠正<\/br><\/br>
第4步反映了VeChat系统的新颖性。第4步在比较第一个循环和第二个循环时有所不同，第一个循环见图2（a）（a），第二个循环见图2（a）（b）。第一个周期的第4步要复杂得多，因为它反映了识别测序错误的关键统计因素。这些关键的统计考虑的核心思想是，真正的变异很可能在不同的读取中同时发生，而错误的发生是随机的。相应地，以下论点是有道理的。由于将读取分为子读取，相应的计算（如变异图构造和相对于读取覆盖率的边的统计评估）可以在子读取之间并行化，这大大加快了计算速度<\/br>
变异图构造<\/br>
我们称之为子读取的窗口中的子序列随后用于构造变异图G=（V，E，P）。这个变化图是一个有向无环图（DAG），其中顶点v2v代表核苷酸（a，T，C，G），边（vi，vj）2e表示由节点vi和vj代表的核苷酸在构建图的其中一个读取中以两个字母的子序列出现，相对于读取对齐堆的特定位置。因此，例如，如果vi和vj分别对应于A和G，那么相对于目标读取所暗示的坐标的读取正好在该特定位置显示AG，从而产生一条边（vi，vj）。相应地，可以将读取标识为变化图中的某些路径P=（v1，…，vl）。对于变异图的构建，我们使用了偏序对齐（POA）算法（Lee等人，2002年）及其更快的版本，并通过SIMD矢量化进行了增强，如（Vaser等人，2017年）所述。
修剪：原则
在下面，我们将使用符号v2v来表示字母表{A，C，G，T}中的字母，一个特定的节点v指的是。影响TGS读取的高错误率和POA算法可能引入的偏差（例如，因为POA取决于相对于读取的顺序），在第一个循环中构建的变化图包含许多虚假顶点和边。为了从错误的边和/或节点（顶点）中修剪图形，我们采用了频繁项集挖掘技术。其基本思想是用项目集识别边，如果相应的项目集似乎不够频繁，则从图中剪除边。删除边后，读取将与结果图重新对齐，因此必须重新计算项集计数 。
这可能会渲染更多的边，以对应不够频繁的项集。重复将边标识为不常见项集、删除它们以及重新对齐读取的循环，直到收敛，也就是说，直到没有其他边被标识为要删除。在实践中，我们确定3是实验的适当迭代次数。请注意，为了根据修改后的图重新对齐读取，我们使用POA算法，但不再重新修改图。
挖掘频繁项目集的基础模型是“市场篮子模型”。篮子对应于项目集，频繁项目集对应于出现在足够多篮子中的项目子集，反之亦然，不频繁项目集对应于没有出现在足够多篮子中的项目子集。
按照这个模型，基本项目集与变化图中的节点集V一致。然后，篮子对应于读，读被建模为路径P=（v1，…，vl），它们包含的项对应于节点v1，vl 2 V读取封面。此外，我们感兴趣的项集对应于边e=（v，w），作为项v，w的对。如果作为项的特定子集的两个连续节点v，w没有出现在足够多的篮子中，也就是说，没有以该特定顺序包含在足够多的路径P中，则相应的边e=（v，w）将从图中移除。
为了恰当地量化“足够多的篮子”，我们进一步使用支持度和置信度作为频繁项集挖掘的两个标准定义。虽然“支持”只对应于包含特定项目子集的篮子数量，“信心”对应于衡量篮子中项目的外观是否相关（或者换句话说，项目集是否相互“关联”）。这里，Support只同意覆盖特定边e=（v，w）的读取次数。置信度对应于覆盖边缘（v，w）的读取量，一方面与覆盖v的读取次数有关，另一方面与覆盖w的读取次数有关。如果既没有足够多的读取覆盖v也覆盖（v，w），也没有足够多的读取覆盖w也覆盖（v，w），我们就对e=（v，w）失去信心，因为边缘（v，w）可能反映排序错误噪声，并将其移除。
删减：定义<\/br><\/br>
为了明确上面的想法，让R（v）是所有覆盖节点v的读取。对于r∈ R（v），让进一步的Pr，v反映v在R中反映错误的概率。在处理FASTQ文件时，概率pr，v来自R的短语配置文件。对于F ASTA文件，pr，v取零。
我们现在想要确定w（v），作为节点v的权重，它反映了覆盖它的预期读取次数。注意，对于F ASTA文件，w（v）只与覆盖v的读取次数一致。对于F ASTQ文件，w（v）对应于总和1㼿 pr，v在所有读取的r∈ R（v）中，作为读取的r∈ R（v）确实反映与v相关的字母的概率之和。根据公式，我们得到
最佳对齐路径提取<\/br><\/br>
修剪算法收敛后，对应于小窗口的目标子读取将与修剪算法最后一次迭代产生的完全修剪的变化图重新对齐。然后，将图中与目标子读取的最佳对齐对应的路径作为预校正的目标子读取；请参见图2（a）中的橙色元素以获取说明。
第5步：串联
在这一步中，预校正的目标子读被串联成一个完整的预校正目标读，这与明显的、直截了当的“拼凑在一起”预校正目标子读的想法相对应；参见“5。图1中的“连接”用于说明。
步骤6：合并已纠正的目标读取
重复步骤2-5，直到所有读取都已纠正。然后将这些步骤产生的整套预校正结果作为第二个循环的输入。
第二个周期：修改步骤1和4
注意，通过周期1生成的预校正读数仍然包含少量但不可忽略的随机错误。循环2个地址以更正这些剩余错误。为此，重复步骤1-6，但在步骤1和4中进行了一些关键修改。请参见图1中的蓝色元素，以了解反映周期2的流程的工作流。具体而言，第1步的修改不仅包括基于最小值计算所有与所有重叠，还包括预校正读取的基准水平对齐（Li，2018）。这大大有助于根据序列身份相关阈值筛选来自相同单倍型的读取重叠。我们使用δ = 0.99，测序错误率为5-10%和δ =0.98 在我们的实验中，测序错误率为15%，这也是我们通常建议的。然后，对步骤4的修改涉及从每个变化图生成单个一致序列，而不是执行迭代图修剪和序列到图重新对齐。为此，使用动态规划算法（称为“最重束算法”，因为遍历算法通过最大权重选择路径）（Lee，2003），如图2（b）所示。请注意，在第二个周期中，从图形生成一致序列是合理的，因为用于纠正目标读取的重叠读取来自同一单倍型。最后，我们获得完全错误纠正的读取作为第二个周期的输出；参见图1中的“所有最终校正读数（周期2）”。