基于长读数和多序列比对的间隙填充方法

毋 东，魏亚伟，罗军伟，敖 山

（河南理工大学计算机科学与技术学院，河南焦作 454000）

0 概述

随着人类基因组计划的完成，基因组数据成为后基因组时代的重要标志，这些海量基因组数据的分析和挖掘深入到生命科学研究的各个方面，是众多生命科学研究的基础，更是计算机科学特别是生物信息学研究所面临的重大挑战。《国家“十三五”科技创新规划》和《国家自然科学基金委“十三五”发展规划》都明确将基因组数据分析以及生物大数据分析作为以“重大计划”形式倡导的研究内容。

目前，随着全基因组测序技术的发展，人们利用不同的序列组装方法获得越来越连续和完整的基因组序列。但是由于测序错误［1］、基因组重复区［2］、多态性等问题的存在，现有的序列组装方法往往输出一个包含间隙（gap）的序列集合，即scaffold［3］集合，其中每一个scaffold 是一条碱基序列，但是其包含一些gap 区域，通常由连续的“N”表示。而这些gap 区域可能对应基因编码或调控区域，这给下游分析中的基因表达和调控、结构变异分析、物种进化等研究增加困难。

为减少scaffold 中gap 的数量，增加序列的连续性和完整性，研究人员设计了gap 填充方法。gap 填充方法利用由测序技术获得的读数（read/序列片段）集合，确定scaffold 集合中gap 区域的序列。当前流行的测序技术主要包括第二代和第三代测序技术。第二代测序技术具有高速度、高通量、低成本等优点，但是第二代测序技术产生的读数长度较短，在一百到五百碱基左右，越来越多的研究表明，短读数对于复杂重复区的gap 区域的填充能力非常有限。而以太平洋生物科学公司的单分子实时测序技术［4］和牛津纳米技术公司的纳米孔单分子技术［5］为代表的第三代测序技术，能够产生长度达到数万碱基的长读数。这些长读数可以跨过基因组序列中大部分的重复区，有效弥补第二代测序技术的不足。但是第三代测序技术的测序错误率达到15%左右，这对gap 填充的质量和效率造成较大的影响。

本文提出一种基于长读数和多序列比对的gap填充方法gapLM。将scaffold 集合转化为不含gap 的序列（contig）集合，并利用比对工具把长读数集合比对到contig 集合上，然后对获得的比对结果进行修正。针对一个具体的gap，识别和提取覆盖该gap 区域的长读数序列，并将这些提取出来的序列分为3 个不同集合。最后利用多序列比对方法和新的策略对3 个不同的序列集合进行一致化，并对该gap 区域进行填充。

1 相关工作

近年来，随着测序技术的发展，研究人员提出许多gap 填充方法。根据gap 填充方法使用的读数类型，本文将gap 填充方法分为两类：一类是基于双端短读数的gap 填充方法；另一类是基于长读数的gap填充方法。gap 填充示例如图1 所示。

图1 gap 填充示例Fig.1 Example of gap filling

1）基于双端短读数的gap 填充方法

基于双端短读数的gap 填充方法往往利用短读数正确率高和组装工具成熟等特点对gap 进行填充，其主要思路是基于读数集合和scaffold 集合的比对信息确定和gap 关联的读数集合，然后再采用不同的策略对gap 进行填充，如图1（a）所示。一些gap 填充方法利用读数或者k-mer（长度为k的子序列）之间的重叠关系，从gap 区域紧邻两端进行扩展找到覆盖gap 区域的序列。文献［6］提出一种IMAGE 的gap填充方法，该方法使用比对工具将双端读数比对到scaffold 集合上，收集一端读数能够比对到gap 紧邻区域的双端读数集合，然后通过Velvet［7］对读数集合进行组装，构建contig（更长的序列）并填充gap 区域。文献［8］提出一种gapCloser 方法，首先确定和gap 相邻区域有重叠关系的读数集合，并在局部组装时加入容错机制［9］，从而提高填充gap 的准确性。文献［10］提出一种gapFiller 方法，该方法将对应到gap区域的读数转换成k-mer 集合，然后从gap 左右两端选择具有重叠关系的k-mer 进行合并，反复扩展逐渐缩小gap 区域。

此外，一些启发式方法也被用于gap 填充的过程。此类方法一般先构建读数重叠图或者De Bruijn图，然后从图中选择一条路径代表gap 区域的序列。文献［11］提出一种FinIS 方法，该方法针对每个gap区域，建立一个读数重叠图，然后采用混合整数规划方法和路径评分函数确定最佳路径。文献［12］提出一种Sealer 方法，该方法针对一个gap，确定一个读数集合和相对应的k-mer 集合，使用布隆过滤器数据结构存储这些k-mer 集合，并构建De Bruijn 图，在De Bruijn 图中利用广度优先搜索寻找连接gap 区域两侧序列的路径。文献［13］提出一种gap2Seq 方法，该方法是在构建的De Bruijn 图中确定最接近gap 长度的路径，由于该过程是NP-hard 问题，因此gap2Seq方法基于一种简单的动态规划算法寻找近似解。文献［14］提出一种gapReduce 方法，该方法通过迭代改变k-mer 频次和长度阈值构建De Bruijn 图，并基于双端读数的insert size 分布设计评分函数选择路径。

2）基于长读数的gap 填充方法

基于长读数的gap 填充方法的基本思想是利用比对方法将长读数与scaffold 集合进行比对，然后选择gap 对应的长读数去填充，如图1（b）所示。文献［15］提出一种PBJelly 方法，该方法采用BLASR［16］将长读 数比对 到scaffold 集合上，识别覆盖gap 区域的长读数集合，并基于OLC［17-18］组装算法找到填充gap 区域的一致序列。文献［19］提出一种GMcloser方法，该方法首先将scaffold 中的gap 区域删除，使其变成不包含gap 的contig 集合，然后将这些contig 与长读数集合进行比对。该方法基于比对概率估算模型，选择最大可能比对到gap 区域的长读数序列进行填充。文献［20］提出一种LR_gapcloser方法，该方法将每个长读数分成长度相同且有序的短序列片段，然后使用BWA-MEM［21］算法将所有的短序列片段比对到scaffold 上。该方法基于碱基覆盖度评分函数以及短序列片段的比对方向和顺序，过滤掉一些错误比对，进而提高填充的准确性。

2 本文方法

本文提出一种基于gapLM 的gap 填充方法，该方法以长读数集合和scaffold 集合作为输入，通过分析长读数集合和scaffold 集合之间的比对信息，对gap 区域进行填充。该方法主要包括以下4 个步骤：1）删除scaffold 中的gap 区域，获得一个新的contig集合，然后利用比对工具BWA 将长读数比对到contig 集合上；2）通过分析contig 和长读数的之间比对位置的差异，对比对结果进行修正；3）根据上一个步骤获得的比对信息，提取对应gap 区域的序列，并将这些序列分为3 种不同的类型；4）基于多序列比对方法，采取一种新的策略将一个gap 区域关联的3 种不同类型的序列进行一致化，获得最终的填充序列。

2.1 初始比对结果

由于scaffold 中的gap 区域是由“N”组成的，并且可能比较长，因此直接把长读数比对到scaffold上，可能会使一些长读数无法比对到gap 的临近区域。因此，首先删除每条scaffold 中的gap 区域，这样一条scaffold 变成多条contig，每条contig 是不包含gap 的长序列。然后利用比对方法BWA 将长读数比对到contig 集合上。

2.2 比对结果修正

经过上一步骤，本文方法可以获得一个比对结果，但是由于长读数的测序错误率较高，并且基因组中往往存在一些重复区，因此现有的比对工具容易造成一些错误的比对信息，或者给出的比对区间位置和真实情况相比具有一些偏差。为了获得更加准确的gap 区域对应的序列，本文方法对每一条比对信息进行判断和修正。

针对一条长读数和一个contig 的比对结果，本文方法抽取它们之间的比对位置、比对方向和比对质量信息，并将一条比对信息定义为一个六元组（Rs，Re，Cs，Ce，Ori，Q）。其中，长读数上的［Rs，Re］区间和contig 上的［Cs，Ce］区间能够比对上，比对方向为Ori（0 代表反向比对，1 代表正向比对），Q代表比对质量分数。如果Q低于一个阈值（默认为30），本文方法认为该比对信息不可信，它不会用于后续的步骤。对于Q高于阈值的比对信息，本文方法计算长读数和contig 之间比对区间的4 个偏移大小，分别是［Rs，Re］区间的左右两端偏移大小DRs、DRe和［Cs，Ce］区间的左右两端偏移大小DCs、DCe。如果4 个偏移量中有1 个大于α（默认500），则认为该比对是错误的比对信息，本文方法忽略该比对信息。如果4 个偏移量均小于α，则本文方法对初始的2 个比对区间进行修正，得到更新后的比对区间，具体过程如算法1 所示。

算法1RevisePosition（Rs，Re，Cs，Ce，Ori，Len（R），Len（C），α）

比对区间信息示例如图2 所示。一条长读数能分别比对到2 条contig 上，其中contig 1 上的区间［Cs1，Ce1］和长读数上的区间［Rs1，Re1］能够正向比对上，contig 2 上的区间［Cs2，Ce2］和长读数上的区间［Rs2，Re2］能够正向比对上，如图2（a）所示；利用算法1 对这2 个比对信息进行修正，得到最终的比对区间信息，如图2（b）所示。

图2 区间信息示例Fig.2 Sample of interval information

2.3 gap 区域对应的3 类序列集合

当对所有的比对信息进行修正后，获得一个新的比对信息集合。针对1 个gap 和其左右相邻的2 个contig，根据比对到这2 个contig 的长读数位置以及方向信息，抽取长读数中对应于该gap 的序列区域。有些长读数可以跨过该gap 区域，并且可以同时比对到这2 个contig 上。而有些长读数只能比对这2 个contig 中的1 个，且只能覆盖gap 区域的左右两侧部分区域。因此，本文方法针对1 个gap，把可能覆盖它的长读数区域抽取出来，并分别保存到3 个不同的集合中：跨过序列集合（SpanSet），左侧序列集合（LeftSet）和右侧序列集合（RightSet）。

针对1 个gap，假设其左右相邻的2 个contig 是Ci和Cj。本文方法首先确定能够比对到Ci和Cj的所有长读数，如果其中1 个长读数能够同时比对到Ci和Cj，这2 个比对 信息分别是Alignment1=（Rsi，Rei，Csi，Cei，Orii，Qi）和Alignment2=（Rsj，Rej，Csj，Cej，Orij，Qj），则本文方法用算法2 判断该长读数是否能够连接Ci和Cj，并且判断该长读数在这2 个contig 之间的序列是否对应该gap 区域。如果算法2 返回false，表明该长读数无法连接Ci和Cj。如果算法2 返回true，则表明该区间对应gap 区域，并把对应gap 区域的序列保存到SpanSet。其中抽取序列的长度不超过该gap 长度加β（默认200）。ReverseComplement是对一个序列进行反向互补操作。

算法2DetermineSpanSet（Alignment 1，Alignment 2，SpanSet，β）

当一个长读数能比对到Ci，但不能比对到Cj时，则该长读数可能覆盖该gap 区域的左侧区域，本文方法用算法3 判断该长读数是否能覆盖gap 区域的左侧，如果能，则抽取对应gap 区域的序列，并保存到LeftSet 集合中。

算法3DetermineLeftSet（Alignment 1，LeftSet，β）

当一个长读数能比对到Cj，但不能比对到Ci时，则该长读数可能覆盖该gap 区域的右侧区域，本文方法用算法4 抽取对应gap 区域的序列，并保存到RightSet 集合中。

算法4DetermineRightSet（Alignment 2，RightSet，β）

根据上述步骤，本文方法获得对应1 个gap 的3 个不同类型的序列集合，注意其中一些序列集合可能为空。

2.4 gap 区域填充

在填充1 个gap 时，关键是对该gap 关联的3 个集合中的序列进行一致化，获得一致序列。racon［22］是一个快速的序列一致算法，在该方法中，首先要确定一个目标序列（targe sequence），然后利用minimap2［23］将一个序列集合比对到该目标序列，并获得它们之间的重叠信息（overlap），最后利用该重叠信息对目标序列进行纠错达到序列一致化的目的。但是当目标序列和序列集合中的序列长度较短时，minimap2 往往无法获得这些序列的重叠信息，从而无法获得一致序列。针对racon 的缺陷，本文方法在racon 无输出结果时，直接利用多序列比对方法mafft［24］获得多条序列的比对信息，然后设计一种一致化策略获得一致序列。本文方法基于序列一致算法racon 和多序列比对方法mafft，采用一种新的策略对一个gap 进行填充。

针对一 个gap 和其对应的SpanSet、LeftSet 和RightSet 集合进行填充，具体的填充步骤如下：

1）如果SpanSet 中的序列个数大于等于3，利用racon 尝试求出一个一致序列。在此过程中，首先利用minimap2 对这些SpanSet 中的序列进行重叠检测，统计每条序列有多少条其他序列能够比对到它本身，并选择获得最多比对的序列作为目标序列，如果同时有多条，则选择和gap 长度最接近的序列作为目标序列。然后以该目标序列、重叠信息和其他序列作为racon 的输入，racon 的输出作为一致序列。如果racon 没有输出任何结果，则利用mafft 对SpanSet 中的所有序列进行比对，并获得这些序列的布局（layout），最后基于少数服从多数的原则，选择每个碱基位置获得最多支持的碱基，并获得一条一致序列。通过上述步骤获得的一致序列用SpanConsensusSequence 表示。

2）如果SpanSet 中的序列个数小于3 大于0，则首先从SpanSet 中找到一个目标序列（方法和上一步找目标序列的方法一样），然后利用minimap2 把3 个序列集合中剩余的序列比对到目标序列上，接着用racon 进行一致化。如果racon 没有输出，则利用mafft 对SpanSet 和LeftSet 中的所有序列进行多序列比对，求出一个一致序列。然后再用mafft 对SpanSet 和RightSet 中的所有序列进行多序列比对，求出一个一致序列，最后对这2 条一致序列进行重叠检测，并合并这2 个序列求出一条新的一致序列。如果LeftSet 或者RightSet 为空集，则直接把目标序列当作一致序列。通过上面步骤获得的一致序列用SpanConsensusSequence 表示。

3）如果SpanSet 中的序列个数为0，即为空集，则利用racon 分别对SpanSet 和RightSet 中的序列进行一致化求解，分别求出一个左侧一致序列和一个右侧一致序列（求解过程和上述的方法相同）。如果利用racon 无法求出其左侧或者右侧一致序列，则利用mafft 分别对LeftSet 和RightSet 中的序列进行多序列比对，并分别求出左侧和右侧一致序列，这2 条一致序列分别用LeftConsensusSequence 和RightConsensus Sequence表示，注意如果LeftSet或者RightSet 为空集，则对应的一致序列为空。

4）如果对于一个gap，通过上述步骤可以获得SpanConsensusSequence，则直接利用该序列填充这个gap。如果对于一个gap只有LeftConsensus Sequence，则根据gap 的长度，从左往右填充这个gap，最多填充gap 长度。如果对于一个gap 只有RightConsensusSequence 序列，则根据gap 的长度，从右往左填充这个gap，最多填充gap 长度。如果对于一个gap 既有LeftConsensusSequence，又有Right ConsensusSequence，则同时从左往右，和从右往左填充gap，最多填充gap 长度。

通过上述步骤，获得一个gap 的最终填充序列。当依次处理完所有的gap 后，则完成整个gap 填充过程。

3 时间复杂度

本文方法包括比对阶段、比对修正阶段、确定gap 对应的3 个序列集合阶段和gap 填充阶段4 个阶段。

1）比对阶段。本文方法调用BWA-MEM 算法将长读数比对到scaffold 集合上，该阶段的时间复杂度为O（w）+O（s），其中，w是所有长读数长度之和，s是所有scaffold 长度之和。

2）对比对修正阶段。本文方法依次调用算法1处理所有的比对信息，其时间复杂度是O（m），m是比对信息个数。

3）确定gap 对应的3 个序列集合阶段。针对每一个gap，本文方法确定所有能够比对其左右两侧contig 的比对信息。如果一个长读数能跨过该gap，则调用算法2；如果一个长读数只比对到该gap 的左侧contig，则调用算法3；如果一个长读数只比对到该gap 的右侧contig，则调用算法4。这个步骤的时间复杂度是O（t×m），其中t是gap 的个数。

4）填充gap 阶段。针对一个gap 区域及其关联的3 个序列集合，利用多序列比对算法racon 和mafft产生一致序列。racon 时间复杂度是O（p×L），mafft时间复杂度是O（p2×L）+O（p×L2），p是进行多序列比对的序列个数，L是序列的长度。

4 实验与结果分析

4.1 实验数据

为验证gapLM 的有效性，将gapLM 在2 个真实Pacibio 测序数据集上和其他3 种常见的gap 填充方法（GMcloser，PBJelly，LR_Gapcloser）进行比较。这2 个数据集分别来自2 个物种：酿酒酵母S288C（S.cerevisiae）和秀丽隐杆线虫（C.elegans）。本文方法采用2 个由NGS 组装产生的scaffold 集合作为原始输入序列，使用不同的gap 填充方法对这些原始scaffold 集合进行填充。这2 个长读数集合和scaffold 集合的详细信息如表1 所示。数据集都来自于文献［20］。

表1 数据集信息Table 1 Datasets information

4.2 评价指标

本文使用QUAST［25］对最终的填充结果进行评估。QUAST 将参考基因组与gap 填充后的scaffold集合进行比对，并获得多个评价指标。组装结果的好坏主要取决于长度和准确度这2 个方面。本文利用QUAST 进行相关评价时，主要列出以下4 项评测指标：

1）N50。contig 或scaffold 按照长度大小进行排序，并依次将其长度相加，当累加值超过所有contig（或者scaffold）长度之和的1/2 时，最后一个contig 或scaffold 的长度就是N50［26］的值。

2）NGA50。在组装结果发生错误的位置断开contig 或scaffold，然后重新计算N50 值。

3）基因组覆盖比例（GF）。将组装好的scaffold或者contig 与参考基因组进行比对后，计算覆盖到参考基因组的碱基百分比。

4）组装错误（MA）。在contig 或scaffold 与参考基因组进行比对后，发生在contig 或scaffold 中错误的个数。

4.3 结果分析

当gapLM、PBJelly、GMcloser 和LR_gapcloser分别对scaffold 集合进行填充后，为了更好地评价填充结果，本文将填充后的scaffold 在gap 区域处打断，产生不包含gap 区域的contig 集合。这样填充越多越正确的gap 填充方法产生的contig 数目也就越少，并且准确性也高。然后利用评价工具QUAST 将这些contig 集合和基因组参考序列进行比对，并获得相应的评价指标。

4 种gap 填充方法最终结果的评价如表2 所示，对于S.cerevisiae 数据集，其具有相对较小的基因组。与参考数据集相比，原始的输入数据集包含29 个错误组装（MA）。gapLM 对可能覆盖gap 区域的序列进行划分，并分别抽取有关序列进行多序列比对产生一致序列。在第1 组数据集上，与其他方法相比，gapLM 获得最大的NGA50值。虽然MA值比GMcloser 要高，但是GMcloser 的GF 值是最低的，填充效果不明显。而对较大的C.elegans 数据集，本文方法获得的gap 填充结果的N50 和NGA50 值都是最优的，并且本文方法的MA 值在4 种方法中排名第2，这说明本文方法能够获得较好的填充结果。

表2 gap 填充评价结果Table 2 Evaluation results of gap filling

在运行时间上，所用gap 填充方法主要消耗在比对阶段，这和选择的比对工具相关。PBjelly 使用BLASR 进行比对，GMcloser 使用MUMcer 进行比对，LR_Gapcloser 和gapLM 都使用BWA-MEM进行比对。但是，LR_gapcloser 把长读数打断成一些小片段，然后将这些小片段比对到scaffold 集合上，这样减少耗费在比对上的时间。gap 填充评价结果如表2 所示，其中，contigs≥1 000 bp。从表2 可以看出，LR_Gapcloser 在2 个数据集上都运行最快，并且消耗的内存最小，而GapLM 在运行时间和最大内存消耗上比GMcloser 和PBjelly 更优。

与其他3 个常见的gap 填充方法相比，gapLM 的N50 和NGA50 值相对较高，并且gapLM 的错误率也较低，这表明填充区域序列的准确度比其他方法具有一定优势。因此，gapLM 产生更加准确和有效的填充结果。虽然GMcloser 方法利用长读数进行填充，但是它填充后的scaffold 集合的覆盖度和原始scaffold 集合的覆盖度一样，说明其没有进行有效填充。PBjelly 方法能够利用长读数对gap 进行填充，但是其错误率显著增加。LR_gapcloser 方法在填充结果上取得不错的效果并且其运行效率最优，但是由于其直接利用长读数上的序列去填充对应的gap区域，因此其结果仍然存在一定的错误率，特别是在大基因组数据上，其填充结果的N50 值较低。

5 结束语

本文提出一种基于长读数和多序列比对的gap填充方法GapLM。将scaffold 集合中的每条contig在gap 区域分割成不包含gap 的contig 集合，然后把长读数比对到contig 上，并对比对结果进行修正以获得更加准确的比对信息。基于比对位置和方向信息，获得可能覆盖一个gap 的3 个序列集合。利用多序列比对方法和一种新的策略，对这3 个序列集合进行一致化，获得一致序列以填充gap 区域。将gapLM 与目前常见的3 种其他gap 填充方法在2 个真实数据集上进行比较，实验结果表明，gapLM 能够产生更加准确的填充结果。但该方法仍然存在运行时间稍长等不足，因此设计一种快速的针对gap 填充的比对方法是下一步的工作重点。另外，为提高gap填充的完整性和准确性，充分挖掘长读数和双端短读数的互补优势，同时利用长读数和双端短读数进行gap 填充也是未来研究的方向。