基于de Bruijn图和序列比对的长序列混合纠错算法

刘 刚

（广西大学计算机与电子信息学院，南宁 530004）

0 引言

基因纠错对于基因工程是很重要的。序列纠错算法旨在识别和消除（或修复）测序序列中包含的错误，从而有利于重新测序或从头测序分析。序列纠错算法应能有效地处理越来越多、越来越长的测序数据的纠错。

为了纠正第三代测序平台产生的错误率高的长序列，人们开发了混合纠错算法。所谓混合纠错算法是指利用第二代测序平台产生的高准确率的短序列来对第三代测序长序列进行纠错的算法。长序列混合纠错算法分为4类：基于短序列比对的混合纠错算法，基于重叠序列或重叠群比对的混合纠错算法，基于de Bruijn 图的混合纠错算法以及基于隐马尔可夫模型的混合纠错算法。基于短序列比对和基于重叠序列比对的混合纠错算法由于长序列的较高错误率导致比对会出现偏差，因此需要改进比对方式来提高纠错效果，且它们将长序列拆分成很多的长度很短的序列片段，这些序列片段会丢失长序列所携带的信息，这使得纠错效果不够理想。基于隐马尔可夫模型的混合纠错算法也需要将长序列与短序列比对，也需要改进比对方式来提高纠错效果。已有的基于de Bruijn 图的混合纠错算法采用固定值的k-mers（长度为的序列片段）构建，其有时很难找到长序列的有效锚点，从而难以扩展序列路径，这使得纠错质量不高。

本文研究基于值可变的de Bruijn 图和序列比对的长序列混合纠错算法，以纠正第三代测序长序列中与第二代测序短序列比对的区域和未被覆盖的区域。

1 算法

1.1 算法思想

本文的算法思想：将与长序列同物种的短序列采用短序列纠错算法其进行纠错，从纠错后的短序列中提取出k-mers（长度为的序列片段），将待纠错的第三代测序长序列与同物种的正确率高的第二代测序短序列进行比对，以生成种子；使用值可变的de Bruijn图扩展连接形成种子序列，将连接两个相邻种子的序列路径覆盖位于两个种子之间未与长序列比对的区域；遍历de Bruijn 图扩展种子序列的末端，使得种子序列能够扩展到待纠错的长序列的末端；采用短序列比对纠正长序列与短序列对准的区域，并使用种子序列路径来纠正长序列未与短序列对准的区域，从而得到经纠错之后的长序列。

1.2 算法描述

算法是通过将待纠错的长序列与同物种的短序列比对来生成种子的，种子由5 元组（id，pos，len，score，seq）组成，其中id 表示种子与长序列关联的标识符，pos 表示相同物种的短序列与长序列比对的开始位置，len 表示短序列与长序列比对区域的长度，score 表示比对得分，seq 表示在某个位置上短序列与长序列比对的共有序列。

本文给出的基于值可变de Bruijn图和序列比对的长序列混合纠错算法（hybrid de Briujin graph errors correction algorithm，简记为HdGEC算法）形式描述如下。

HdGEC

输入：相同物种的长序列long-read[0..n-1]和短序列short-reads[0..m-1]，参考基因组R[0..r-1]

输出：纠正后的长序列long-read[0...n-1]

Begin

1:采用QuorUM 算法纠正短序列short-reads[0..m-1］，以减少短序列中包含的错误；

2: 使用KMC3 工具从纠错后的短序列short-reads［0..m-1］中提取出k-mers；

3:过滤掉弱值的k-mers，将强值的k-mers 使用PgSA索引构造值可变的de Bruijn图；

4: 使用BLASR 工具将经过纠错的短序列shortreads[0..m-1］与同物种的长序列long-read[0...n-1]比对，以获得种子seeds（id,pos,len,score,seq）；

5: 利用序列比对将种子序列seeds（id, pos, len,score,seq）与参考基因组R[0..r-1]进行比对，以纠正中被短序列short-reads[0..m-1］覆盖到的区域；

6:遍历值可变的de Bruijn 图以找到任意相邻两个种子之间的序列路径，连接这相邻的两个种子得到种子序列路径，并用种子序列路径来纠正长序列long-read[0...n-1]当中未被短序列short-reads[0..m-1]覆盖到的区域；

7:遍历值可变的de Bruijn 图，以扩展得到种子序列路径的末端，使得种子序列路径覆盖到整条长序列long-read[0...n-1]；

8:输出纠正之后的长序列long-read[0...n-1]；

End.

本文算法HdGEC 使用种子作为锚点进行扩展，使用贪心策略遍历de Bruijn 图，由于可以选择1～阶的de Bruijn 图来连接种子序列，所以可以减少de Bruijn 图的遍历次数，从而缩短了序列纠错需要的时间。此外，本文算法HdGEC 结合k-mers 比対和值可变de Bruijn 图来对长序列进行纠错，其中通过k-mers 序列比対来纠正长序列中被短序列覆盖到的区域，通过de Bruijn 图纠正长序列中没有被短序列覆盖到的区域，这使得算法的覆盖率更高，能够纠正序列中更多的错误碱基，获得了更高的吞吐量（输出的纠错序列的质量）。

2 实验

2.1 实验环境与数据

实验使用的计算机是8 核Intel（R）Core i7-6700 k CPU@4.00 GHz处理器、内存容量32 GB。采用C++语言编程实现算法。实验数据集采用太平洋生物科学平台的PacBio模拟数据集和牛津纳米孔平台的ONT真实数据集。模拟数据集与真实数据集都来自于物种A.baylyi、E.coli、S.cere⁃visiae和C.elegans 中的数据，这4 种物种数据集和参考基因组来源于基因测序中心https://www.genoscope.cns.fr/externe/NaS/datasets。模拟的长序列数据集采用SimLord模拟器生成。实验使用的参考基因组、模拟数据集与真实数据集的信息如表1所示。

表1 实验数据集

2.2 模拟数据集上的长序列纠错实验

对于模拟数据集上的实验，采用纠正后的长序列的错误率（%）、吞吐量、split reads比率、运行时间这4个指标测试长序列混合纠错算法的性能，其中、和的值越小越好，值越大越好。使用LRCStats 软件测量得到算法CoLoRMap、HALC、Jabba、LoRDEC、NaS和HdGEC 在物种A.baylyi和E.coli 模拟数据集上的实验结果，分别如表2和表3所示。

从表2 与表3 的实验结果可以看到：在纠正后的长序列的错误率上，6 种混合纠错算法的值均低于0.3%，本文算法HdGEC 对物种A.baylyi 纠正后的长序列的错误率最低，对物种E.coli 纠正后的长序列的错误率是第二低的，它比Jabbaa 算法纠正后的长序列的错误率0.0462%仅高了0.0134%。在吞吐量上，除了NaS算法，其他5种混合纠错算法的指标值都较高，并且本文算法HdGEC 的吞吐量高于其他5 种混合纠错算法。在split reads 比率上，算法CoLoRMap、HALC和LoRDEC 的srr值较高，值高意味着长序列的大部分区域没有被纠正，本文算法HdGEC 的值较低，表明参考基因组与短序列比对的区域有很大的机会被发现，长序列被短序列覆盖到的更多区域将会被纠正；在运行时间方面，纠错算法Jabba 所需时间最少，而纠错算法NaS所需时间最多。

表2 长序列纠错算法在物种A.baylyi上模拟数据集的实验结果

表3 纠错算法在物种E.coli上模拟数据集的实验结果

这表明，对于物种E.coli和A.baylyi 的模拟数据集的长序列纠错，本文算法HdGEC 的吞吐量高于其他算法，并且错误率较低，产生的split reads 的比例也较小，能纠错长序列的大部分区域。

2.3 真实数据集上的长序列纠错实验

对于在牛津纳米孔真实数据集上的实验，采用纠正的长序列数量、比率（%）、纠正的长序列平均长度、纠正的碱基数量、平均同一性、覆盖率、运行时间这7个指标测试长序列混合纠错算法的性能，其中、、、和的 值 越 大越好，和值越小越好。长序列纠错算法HdGEC、 CoLoRMap、 HALC、 Jabba、LoRDEC和NaS在小型物种A.baylyi和E.coli、中型物种S.cerevisiae、大型物种C.elegans 真实数据集上运行的实验结果分别如表4～表7所示。

表7 算法在大型物种C.elegans真实数据集运行的实验结果

从表4和表5可以看到，对于小型物种真实数据集的长序列纠错，在6个算法中，本文算法HdGEC 纠正的长序列的平均长度是最长；本文算法HdGEC的split reads比率第二低，仅比算法NaS高了一点点；本文算法HdGEC、CoLoRMap、HALC、LoRDEC和NaS 的覆盖率均达到100%；算法HALC 能纠正的长序列数量最多，本文算法HdGEC 能纠正的长序列数量第二多；在纠正的碱基数量方面，本文算法HdGEC和NaS 是较多的；在运行时间方面，算法Jabba 所需时间最少，本文算法所需时间是第三少。

表4 算法在小型物种A.baylyi真实数据集的实验结果

表5 算法在小型物种E.coli真实数据集的实验结果

从表6可以看到，对于中型物种真实数据集的长序列纠错，在6个算法中，本文算法HdGEC 能纠正的长序列平均长度最长、能纠正的碱基数量最多、split reads 比率最低、覆盖率最高；在能纠正的长序列数量方面，算法HALC最多，本文算法HdGEC 第二多；在平均同一性方面，算法Jabba 最高，本文算法HdGEC 第三高；在运行时间方面，Jabba 算法所需时间最少，本文算法HdGEC 所需时间是第四少；而算法NaS运行时间超过了16天尚未计算得到结果。

表6 算法在物种S.cerevisiae真实数据集的实验结果

从表7可以看到，对于大型物种真实数据集的长序列纠错，在6个算法中，本文算法HdGEC 能纠正的长序列数量最多、能纠正的长序列平均长度最长、能纠正的碱基数量最多、split reads 比率最低、覆盖率最高；在平均同一性方面，算法Jabba 最高，本文算法HdGEC 第二高；在运行时间方面，Jabba 算法所需时间最少，本文算法HdGEC 所需时间是第三少；而算法HALC和NaS 运行时间超过了16 天尚未计算得到结果。

上述结果表明：对于中型和大型物种真实数据集的长序列纠错，本文算法HdGEC 的整体纠错质量最好；对于小型物种真实数据集的长序列纠错，本文算法HdGEC 的整体纠错质量也较好。

3 结语

为了有效纠正第三代测序平台产生的错误率高的长序列，本文给出一个基于值可变de Bruijn图和序列比对的混合纠错算法。在模拟数据集和真实数据集上的实验结果表明，与已有的长序列混合纠错算法相比，本文算法获得了整体较好的纠错质量。本文算法为了获得较高质量的长序列纠错结果，是利用de Bruijn 图来寻找种子序列路径，这需要花费较时间长。下一步工作是研究将本文的算法并行化，以在维持获得较高质量的长序列纠错结果的同时，显著减少纠错过程所需的时间。