开放科学(资源服务)标志码(OSID)：

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填充方法

基于双端短读数的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的比对结果，本文方法抽取它们之间的比对位置、比对方向和比对质量信息，并将一条比对信息定义为一个六元组(R_s，R_e，C_s，C_e，Ori，Q)。其中，长读数上的[R_s，R_e]区间和contig上的[C_s，C_e]区间能够比对上，比对方向为Ori(0代表反向比对，1代表正向比对)，Q代表比对质量分数。如果Q低于一个阈值(默认为30)，本文方法认为该比对信息不可信，它不会用于后续的步骤。对于Q高于阈值的比对信息，本文方法计算长读数和contig之间比对区间的4个偏移大小，分别是[R_s，R_e]区间的左右两端偏移大小DR_s、DR_e和[C_s，C_e]区间的左右两端偏移大小DC_s、DC_e。如果4个偏移量中有1个大于α(默认500)，则认为该比对是错误的比对信息，本文方法忽略该比对信息。如果4个偏移量均小于α，则本文方法对初始的2个比对区间进行修正，得到更新后的比对区间，具体过程如算法1所示。

算法1    RevisePosition(R_s，R_e，C_s，C_e，Ori，Len(R)，Len(C)，α)

输入   一条比对信息(R_s，R_e，C_s，C_e，Ori)，读数长度Len(R)，contig长度Len(C)和α

输出   如果比对正确则修正该比对信息，并返回true，否则返回false

1.if(Ori == 1)

2.if(R_s < C_s)

3.DC_s←R_s；DR_s←R_s；

4.else

5.DC_s←C_s；DR_s←C_s；

6.if(Len(R)–R_e > Len(C)–Ce)

7.DCe←Len(C)–C_e；DRe←Len(C)–C_e；

8.else

9.DR_e←Len(R)–R_e；DC_e←Len(R)–R_e

10.else

11.if(R_s < Len(C)–C_e)

12.DC_e ← R_s；DR_s ←R_s；

13.else

14.DR_s ← Len(C)–C_e；DC_s ← Len(C)–C_e；

15.if(Len(R)–R_e > C_s)

16.DR_e ← C_s；DC_s ← C_s；

17.else

18.DR_e ← Len(R)–R_e；DC_e ← Len(R)–R_e；

19.if(DR_s < α & & DR_e < α & & DC_s < α & & DC_e < α)

20.R_s ← R_s–DR_s；

21.R_e ← R_e+DR_e；

22.C_s ← C_s–DC_s；

23.C_e ← C_e+DC_e；

24.return true；

25.else

26.return false；

比对区间信息示例如图 2所示。一条长读数能分别比对到2条contig上，其中contig 1上的区间[C_s1，C_e1]和长读数上的区间[R_s1，R_e1]能够正向比对上，contig 2上的区间[C_s2，C_e2]和长读数上的区间[R_s2，R_e2]能够正向比对上，如图 2(a)所示；利用算法1对这2个比对信息进行修正，得到最终的比对区间信息，如图 2(b)所示。

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是C_i和C_j。本文方法首先确定能够比对到C_i和C_j的所有长读数，如果其中1个长读数能够同时比对到C_i和C_j，这2个比对信息分别是Alignment1=(R_si，R_ei，C_si，C_ei，Ori_i，Q_i)和Alignment2=(R_sj，R_ej，C_sj，C_ej，Ori_j，Q_j)，则本文方法用算法2判断该长读数是否能够连接C_i和C_j，并且判断该长读数在这2个contig之间的序列是否对应该gap区域。如果算法2返回false，表明该长读数无法连接C_i和C_j。如果算法2返回true，则表明该区间对应gap区域，并把对应gap区域的序列保存到SpanSet。其中抽取序列的长度不超过该gap长度加β(默认200)。ReverseComplement是对一个序列进行反向互补操作。

算法2    DetermineSpanSet(Alignment 1，Alignment 2，SpanSet，β)

输入   比对信息Alignment 1和Alignment 2，gap长度gapDistance，序列集合SpanSet和一个参数β

输出   true或者false

1.if(Ori_i == Ori_j & & Ori_i == 1)

2.if(R_ei > R_sj & & | R_sj-R_ei| < gapDistance + β)

3.sequence ← [R_ei，R_sj]；

4.保存sequence到集合SpanSet；

5.return true；

6.if(Ori_i == Ori_j & & Ori_i == 0)

7.if(R_si > R_ej & & | R_si-R_ej| < gapDistance + β)

8.sequence ← [R_ei，R_sj]；

9.ReverseComplement(sequence)，并保存到SpanSet；

10.return true；

11.return false；

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

算法3    DetermineLeftSet(Alignment 1，LeftSet，β)

输入   一个比对信息Alignment 1，gap长度gap Distance，序列集合LeftSet和一个参数β

输出    true或者false

1.if(Ori_i == 1)

2.d ← Min(Len(C_i) - R_ei，gapDistance + β)

3.保存[R_ei，R_ei+d]到LeftSet；

4.return true；

5.if(Ori_i == 0)

6.d ← Min(R_si，gapDistance + β)

7.保存ReverseComplement([R_si-d，R_si])到LeftSet；

8.return true；

9.return false；

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

算法4    DetermineRightSet(Alignment 2，RightSet，β)

输入   比对信息Alignment 2，gap长度gapDistance，序列集合RightSet和1个参数β

输出    true或者false

1.if(Ori_j == 1)

2.d ← Min(R_sj，gapDistance + β)

3.保存[R_sj-d，R_sj]到RightSet；

4.return true；

5.if(Ori_j == 0)

6.d ← Min(Len(C_j)-R_ej，gapDistance + β)

7.保存ReverseComplement([R_ei，R_ej + d])到RightSet；

8.return true；

9.return false；

根据上述步骤，本文方法获得对应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(p²×L)+O(p×L²)，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]。

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，这说明本文方法能够获得较好的填充结果。

在运行时间上，所用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填充也是未来研究的方向。