CN104200133B - 一种基于读数和距离分布的基因组De novo序列拼接方法 - Google Patents

Abstract

本发明公开了一种基于读数和距离分布的基因组De novo序列拼接方法，采用De Bruijn图存储读数之间的重叠关系，基于读数分布提出了一种新的打分函数用在contig构建、scaffolding和填充空白区域等步骤。该打分函数充分考虑了测序深度，k‑mer频次以及在复杂重复区中insertsize的偏移。本发明简单易用，在不同的模拟和真实测序数据上表现出良好的拼接结果，较其他序列拼接方法具有更高的连续性和完整性。

Description

一种基于读数和距离分布的基因组De novo序列拼接方法

技术领域

本发明涉及生物信息学领域,特别是一种基于读数和距离分布的基因组De novo序列拼接方法。

背景技术

基因组一般是指全部编码和非编码的脱氧核糖核酸(DNA)序列。DNA序列是遗传信息的载体和蛋白质氨基酸序列合成的依据，并引导生物发育与生命机能运作。因此DNA序列是生命存在和发展的决定性因素，生命活动中发生的所有事情，都与DNA序列有着必不可分的联系。在基础生物学研究和众多应用领域中，如诊断，生物技术，法医生物学，生物系统学中，DNA序列已成为不可缺少的知识。基因组测序是指分析特定DNA片段的碱基序列，也就是腺嘌呤(A)、胸腺嘧啶(T)、胞嘧啶(C)与鸟嘌呤(G)的排列方式。测序技术是获得完整基因组序列必要途径。下一代测序技术(NGS)具有速度快，成本低，高通量等特点。De novo序列拼接是在不知道基因组参考序列的情况下，利用由测序技术产生的序列片段(读数或read)重新构建完整原始基因组序列的方法。如何设计有效地序列拼接方法构建连续和正确的原始基因组序列已经越来越引起人们的关注。

基因组的重复区是指在基因组序列中出现两次或两次以上序列片段。重复区的存在是造成序列拼接结果令人不满意的最大原因，也是序列拼接组装过程中最难解决的问题。目前通过下一代测序技术得到的读数长度一般比较短，在100碱基对(bp)左右。这种短读数可以分为两种类型：单个读数和双端读数。单个读数是在测序的时候复制一个比较短的基因组片段，然后对该片段进行测序得到一个读数。另外一种是在测序时，首先复制一个比较长的基因组片段，然后对基因组片段的左端和右端一段短区域进行测序，得到一对短读数即双端读数。双端读数中每对读数之间的间距称为insertsize，即首先复制的基因组片段长度，一般假insertsize长短服从正态分布。由于双端读数的insertsize一般比读数本身的长度长，可以跨过一些重复区，因此双端读数可以克服部分长度小于insertsize的简单重复区。因此，现有的序列拼接方法往往利用双端短读数进行序列拼接。

序列拼接的步骤一般包括四个步骤：(1)数据的预处理阶段。该阶段通过特定的纠错方法，移除或者改正测序数据中的错误碱基；(2)基因组连续片段(contigs)生成阶段。该阶段利用读数集合拼接成contigs；(3)超长序列片段(scaffold)组装阶段。该阶段使用双端读数，确定contigs之间的方向和位置关系，生成scaffolds。(4)填充空白区域(gap)阶段。利用双端读数对scaffolds中的空白区域进行填充。

现有的序列拼接方法可以分为两大类：

(1)De Bruijn图方法：基于De Bruijn的拓扑结构，删除测序错误的结点，选择合适路径确定contigs，形成scaffolds，并填充空白区域。Velvet方法首先建立De Bruijn图并对图中的长结点进行标记，然后利用双端读数标注两个长结点之间中间结点，通过这些标注的结点寻找一条正确的路径连接这两个长结点。但是Velvet针对一些短重复区却无能为力，因为只有结点的长度大于等于读数的长度时，才能够进行标注，所以，这就容易漏掉了一些长度小于读数长度的短结点。SOAPDenovo方法直接利用De Bruijn图中的长结点作为contig，然后再进行scaffolding和fillgap步骤，该方法并没有利用双端读数对contig进行左右延展，这往往会造成contig太短，影响后续的步骤。ALLPath首先填充双端读数之间的空白区域，但是当双端读数在重复区时，往往会造成一些错误的填充。IDBA采用迭代的方法，由小到大的改变k-mer长度，并不断的改变De Bruijn图，消除掉一些重复区和测序错误。但是，在一些重复区长度大于读数长度时，无法通过改变kmer大小消除重复区。Paired-De Bruijn图方法在构建De Bruijn图的时候，利用双端读数对结点之间先后顺序进行标注，但是当双端读数之间的距离变化比较大的时候，该方法的准确性会下降。

(2)Overlap-Layout-Consensus方法：根据读数之间的重叠信息对选定的种子序列进行左右扩充。SSAKE、VCAKE和SHARCGS主要利用读数之间的重叠信息，计算每个碱基位置上匹配的读数个数，一般选择匹配最多的后续碱基进行扩充。PE-Assembly方法首先利用k-mer出现的频次对种子进行筛选，在进行扩充的时候，主要还是利用双端读数找到左右端的碱基。但是，该方法并没有考虑测序深度的问题，有些区域测序深度比较高，而有些区域测序深度比较低，如果不考虑测序深度问题，必定会对序列拼接造成影响。TeleScoper提出了一个打分函数，在进行左右扩展时，构建了一个read-overlap图，并对可能正确的扩展路径进行打分，选择得分高的路径继续进行扩充。但是，该方法在设计打分函数的时候，没有考虑测序深度。

目前，虽然已有的方法利用双端读数进行序列拼接，但是，仍然有两个问题需要进一步进行解决：(1)现有的方法利用由双端读数产生的局部读数集合能够消除一部分重复区带来的影响，但是无法克服一些复杂重复区的影响。(2)现有的方法往往假设测序深度是均衡的，而实际上测序的深度往往是不均衡的。IDBA-UD方法虽然提出了测序深度问题，但是该方法仅仅从相邻结点之间的k-mer频率比值来判断结点本身的正确性和删除错误结点，而没有在contig扩充，scaffolding和填充gap三个步骤中考虑测序深度。

发明内容

本发明所要解决的技术问题是，针对上述现有技术的不足，提供一种基于读数和距离分布的基因组De novo序列拼接方法。

为解决上述技术问题，本发明所采用的技术方案是：一种基于读数和距离分布的基因组De novo序列拼接方法，包括以下步骤：

1)输入双端读数文库，构建初始De Bruijn图，并对初始De Bruijn图进行优化；

2)以De Bruijn图为基础，选择种子序列，并利用打分函数对候选扩展序列进行打分，选择得分最高的候选扩展序列与种子序列合并，并继续进行扩充，直到结束条件为止。扩展后的每个种子序列即为一条contig，所有的contig构成一个contig集合；

3)建立scaffold图，每个结点代表一个contig，边代表两条contig在真实基因组序列上的位置紧邻；

4)填充scaffold图中有边相连的两个结点之间的空白区域，通过匹配上的双端读数长生一个局部读数集合，在局部读数集合上构建新的De Bruijn图，并在该De Bruijn图上寻找能够连接两个结点的路径，如果存在这样的路径，则以该路径填充空白区域。

所述步骤1)中，构建初始De Bruijn图的过程如下：依次读取读数文库中的每条读数，每条读数分解成r-k+1个k-mer，在De Bruijn图中查找每个k-mer对应的结点，如果某k-mer对应的结点不存在，则在De Bruijn图中创造一个新的结点并对应于该k-mer；如果一个k-mer最后k-1个碱基和另一个k-mer前k-1个碱基相同，则在这两个k-mer对应的结点之间加一条有向边；当遍历完所有的读数时，初始De Bruijn图建立完成；其中r为每条读数的长度；k为k-mer的长度；且k为奇整数。

对所述初始De Bruijn图进行优化的方法为：首先De Bruijn图中的简单路径合并为一个结点；简单路径是指一个路径中起始结点的出度为1，最终结点的入度为1，并且所有中间结点的出度和入度均为1；删除初始De Bruijn图中长度小于2*k，并且出度或入度为0的结点。

所述步骤2)中，contig集合构建的具体步骤为：

1)将初始De Bruijn图中结点长度大于双端读数文库insertsize的结点为种子序列，以所述种子序列的前驱结点和后继结点为候选扩展序列；其中，insertsize为双端读数文库中每对读数之间的间距；

2)如果候选扩展序列的长度小于minLen，则以该候选扩展序列为起点，采取深度优先遍历方法，抽取出所有长度大于minLen的子路径，每个子路径对应于一个二级候选扩展序列，对所有的二级候选扩展序列打分，并以最高分作为候选扩展序列的得分值；minLen取值为2*k；

3)针对一个种子序列ss和一个候选扩展序列ec，采用公式F(ss,ec)＝RM(ss,ec)^MC(ss,ec)

计算ec的得分值，其中RM是匹配上ss和ec的双端读数的相对值；由匹配上的双端读数得到一个距离集合，MC用于衡量该距离分布和整体insertsize分布的偏移大小；如果在ec中的实际存在读数个数不符合二项式分布B(minLen，p)，则将所述ec的得分值减去一个罚分值其中a是罚分系数，b为ec中的k-mer平均频次与读数文库中的k-mer平均频次的比值；re是在ec中实际存在的读数个数，p是读数出现在读数文库中的概率；

4)如果一个种子序列只有一个候选扩展序列，则将该候选扩展序列直接合并到种子序列中；如果一个种子序列存在多个候选扩展序列，并且得分最高的两个候选扩展序列得分值都大于0.7或者相差小于0.2，或者没有候选扩展序列，则扩展终止，否则把得分最高的候选扩展序列合并到种子序列中去；迭代地对种子序列进行左右扩展，直到达到上述的扩展终止条件，最终的种子序列作为一条contig加入到contig集合中；

5)当扩展完所有的种子序列，则contig集合构建完成。

所述步骤3)中，确定contig集合中contig位置关系的过程为：

1)建立scaffold图，每个contig作为一个结点；

2)对于contig集合中的任两个contig：contig₁和contig₂，先由与contig₁和contig₂分别匹配的双端读数估计出gap的距离大小，即gapDistance；以contig₁最后insertsize-gapDistance个碱基为测试区域，把该测试区域分成若干个长度为minLen+r-1的子区域，以每一个子区域为候选扩展序列，并以contig₂为种子序列，利用公式对每个子区域打分，其中RN为子区域中存在双端读数的读数个数，SN为双端读数匹配上种子序列的读数个数；如果任何一个子区域的RM小于0.2，则认为contig₁和contig₂之间不存在相邻关系，contig₂作为contig₁右邻居的得分值为0；如果所有子区域的RM均大于0.2，则contig₂作为contig₁右邻居的得分值为所有子区域的平均RM值；当计算完contig₁和所有其它contig之间的得分值之后，选择得分最大的contig作为contig₁的右邻居，并在scaffold图中contig₁和contig₁的右邻居这两个结点之间增加一条边，同理计算contig₁左邻居得分值；

3)当完成所有contig的左右邻居的计算，并得到初始的scaffold图后，则删除错误的边。若contig₁的得分值最大的右邻居是contig₂，contig₂得分值最大的左邻居是contig₁，则认为contig₁和contig₂紧邻，两者之间的边保留；否则，删除它们之间的边。

所述步骤4)中，空白区域填充方法如下：

1)针对两个contig₁和contig₂，收集与所述contig₁和contig₂末端区域分别匹配的所有双端读数，构成一个局部读数集合；

2)以该局部读数集合和k_new构建一个De Bruijn图，并在该De Bruijn图中找到对应于contig₁右端的结点node₁，和对应于contig₂左端的结点node₂；k_new＝k+5，如果k+5大于r，则k_new取r-1；

3)从node₁出发，深度优先遍历该De Bruijn图，选择连接node₁和node₂的路径；在深度优先遍历中，以contig₁为种子序列，用公式F(ss,ec)＝RM(ss,ec)^MC(ss,ec)对候选扩展序列打分，利用公式进行罚分；并以contig₂为种子序列，以公式RM对候选序列打分，如果所有的得分值都大于0.3，则一条路径继续遍历，否则终止该路径；每次都以得分值最高的结点优先遍历，如果候选扩展序列为node₂，则遍历的这条路径作为结果填充node₁和node₂之间的区域；整个填充空白区域过程结束；如果没有找到这样的路径，则令k_new＝k_new ^*，转到步骤2)，直到k_new小于0或者k_new小于k-5；其中k_new ^*＝k_new-2。

与现有技术相比，本发明所具有的有益效果为：本发明的方法不仅利用paired-end读数和kmer频次来消除重复区，并且引入读数分布构建正确的contig，该方法能够有效的解决序列拼接数据中普遍存在的测序深度不均衡和复杂重复区的问题。

附图说明

图1为本发明方法流程图；

图2为本发明一实施例读数分布图；

图3为本发明一实施例contig扩展示意图；

图4为本发明一实施例填充空白区域图。

具体实施方式

如图1所示，本发明具体实现过程如下：

一、构建De Bruijn图

读入fasta格式的读数文库，读数文库中的每个读数的长度相同，并且双端读数的左右读数依次出现在文库中，并且分别对应于正反链。在所有文库中只存在四种碱基{A,T,G,C}。每个读数都是一个长度一定的字符串。每个k-mer是长度为k的字符串。长度为r的读数共存在r-k+1个k-mer。在初始De Bruijn图中每个结点对应于一个kmer。

依次读取每条读数，针对每条读数的r-k+1个k-mer分别在De Bruijn图中找到该k-mer的位置，如果该k-mer不存在，则在De Bruijn图中添加该结点。每两个在该读数中相邻的k-mer，即一个k-mer的最后k-1个碱基和另一个k-mer的前k-1个碱基相同，则这两个k-mer连接一条有向边。当遍历完所有的读数时，初始的De Bruijn图建立完成。

De Bruijn图中的简单路径合并为一个结点，即一个路径中起始结点的出度为1，最终结点的入度为1，并且所有中间结点的出度和入度均为1，这样的路径为简单路径。当处理完所有的这些简单路径后，De Bruijn图中结点个数会减少并且结点长度会增加。针对长度小于2*k长度的结点，并且该结点的出度或者入度为0，则这些结点被认为是由测序错误造成错误结点。这些错误结点全部被删除。k-mer的长度本方法选择为读数长度的2/3左右，并且为奇数。

二、构建contig集合

在De Bruijn图中结点长度越长，则该结点是非重复区和正确性的可能性越大。本方法选择长度大于insertsize的结点为种子序列进行扩充，只有大于该距离才能利用双端读数信息。

构建contig就是以种子序列为起点，基于De Bruijn图寻找左右后续结点，并进行合并，迭代的进行扩充，直到结束条件为止。

针对一个长度为len的序列，则该序列中存在len-r+1个读数，每个读数对应一个起点，即len-r+1个起点碱基位置，如果一个读数出现在文库中，那么这个读数对应的起点位置被标记，那么这len-r+1个位置被标记的个数符合二项式分布。假设每个碱基位置在文库中被标记的概率为p，则该序列区间中被标记的碱基个数的均值为(len-r+1)*p，方差为(len-r+1)*p*(1-p)。p值可以利用De Bruijn图中长结点进行估计，收集De Bruijn图中每个长度大于1000的结点，然后计算被标记的碱基位置个数和所有的碱基位置个数，利用这个比值计算估计p值。

针对当前种子序列，为了能够利用双端读数进行判断后续结点的正确性。每个候选扩展序列必须大于minLen长度。如果候选扩展序列的长度小于minLen，则利用深度优先遍历的方法，找到所有以该候选扩展序列为起点，长度大于minLen的结点序列。每个结点序列都对应于一个二级候选扩展序列。对所有的二级候选扩展序列打分，并以最高分作为候选扩展序列的得分值。minLen默认取值为2*k。

针对当前种子序列ss和一个候选扩展序列ec，统计候选扩展序列中被标记的碱基位置个数(MN)，在被标记的碱基位置所对应的读数存在左读数的个数(RN)，以及这些左读数能够匹配上种子序列的被标记的碱基位置个数(SN)，利用这些信息计算相对匹配得分(RM)。

当候选扩展序列在复杂重复区时，通过RM很难消除这些重复区的影响。我们通过MC对RM的可信度进行评价：

insertsize_candidate是在候选扩展序列中所有匹配上的双端读数之间的平均距离。insertsize是读数文库中的双端读数的距离。SD是读数文库中的双端读数之间距离insertsize的标准方差。

相对于当前种子序列，候选扩展序列的得分值为：

F(ss,ec)＝RM(ss,ec)^MC(ss,ec) (3)

如果候选结点中被标记的碱基位置个数不在[avg-3*sd,avg+3*sd]，avg为(minLen-r+1)*p，sd为(minLen-r+1)*p*(1-p)，则认为该候选结点是一个错误的扩展，需要对其进行罚分：

F(ss,ec)＝RM(ss,ed)^MC-P(ss,ec) (6)

a是一个罚分系数，取0.2，AKF₁是在候选扩展序列中存在的k-mer的平均频次。AKF₂是在整个文库中所有k-mer的平均频次。由于在测序深度比较低的区域，被标注的碱基位置个数本身就会比较小，所以，当对其进行罚分时，应该不进行罚分或者尽量小的罚分。本方法结合标记碱基位置个数和k-mer频次判断候选扩展序列是一个测序深度低的区域还是一个错误的候选扩展序列。

当存在多个读数文库时，只要当前种子序列的总长度大于该文库的insertsize，即可以利用该文库对候选结点进行打分。特别是insertsize比较大的时候，对于消除重复区非常有帮助。最后计算候选扩展序列得分值是每个文库得分值的平均值。

当前种子序列只有一个候选扩展序列时，那该候选扩展序列直接合并到当前种子序列中。如果当前种子序列存在多个候选结点时，根据上述打分函数对每个候选扩展序列进行打分，按照得分值进行由大到小进行排序，如果前两个得分值都大于0.7或者它们相差小于0.2，则contig扩充停止。否则选择得分最大的候选扩展序列作为当前种子序列的正确候选结点，合并到当前路径，继续进行扩充。

当左右均进行扩充完毕，则选择下一个种子序列进行左右扩充。当所有的种子序列扩充完毕，则初步contig集合构建完毕。

针对生成的初步contig集合，计算每个contig之间是否有包含关系以及重叠关系。如果一个contig包含另一个contig，则删除短的contig。如果两个contig₁和contig₂之间的首尾重叠长度len大于读数长度，则contig₁的最后insertsize个碱基去掉重叠部分len个碱基，构成一个长度为insertsize-len的区域，该区域连续划分为一些长度为r+k的子区域，计算每个子区域和contig₂的RM得分值。同理，contig₂的前insertsize个碱基去掉重叠部分len个碱基，也构成一个长度为insertsize-len的区域，并划分成长度为r+k的子区域。计算每个子区域和contig₁的RM得分值。如果所有子区域的RM值均大于0.2，则两个contig合并。形成最终的contig集合。

三、Scaffolding

在该步骤中，建立scaffold图，确定所有contig之间的先后顺序。

首先针对每个contig建立一个结点。对于contig集合中的任两个contig：contig₁和contig₂，先由能够分别匹配上它们的双端读数估计出gap的距离大小即gapDistance。以contig₁最后insertsize-gapDistance个碱基为测试区域，把该测试区域分成若干个长度为minLen+r-1的子区域，以每一个子区域为候选扩展序列，并以contig₂为种子序列，利用公式对每个子区域打分，其中RN为子区域中存在双端读数的读数个数，SN为双端读数能够匹配上种子序列的读数个数。如果任何一个子区域的RM小于0.2，则认为contig₁和contig₂之间不存在相邻关系，contig2作为contig1右邻居的得分值为0；如果所有子区域的RM均大于0.2，contig2作为contig1右邻居的得分值为所有子区域的平均RM值。当计算完contig₁和所有其它contig之间的得分值之后，则选择得分最大的contig作为contig₁的右邻居，并在scaffold图中这两个结点之间增加一条边。同理，可以计算contig₁左邻居得分值。

当完成所有contig之间位置顺序的计算后，则需要删除掉一些错误的边。若contig₁的得分值最大的右邻居是contig₂，contig₂得分值最大的左邻居是contig₁，则认为contig₁和contig₂紧邻，两者之间的边保留。否则，删除它们之间的边。

四、FillGap

在scaffold中，两个相连结点之间存在一个空白区域。本步骤利用双端读数对该空白区域进行填充。首先，针对两个contig₁和contig₂，收集能够匹配上它们末端区域的所有双端读数，构成一个局部读数集合。

以该局部读数集合和k值构建一个De Bruijn图。并在该图中找到对应于contig₁右端的结点node₁，和对应于contig₂左端的结点node₂。

以该局部读数集合和k_new值构建一个De Bruijn图，并在该De Bruijn图中找到对应于contig₁右端的结点node₁，和对应于contig₂左端的结点node₂。k_new等于k+5，如果k+5大于r，则k_new取r-1.

从node₁出发，深度优先遍历该De Bruijn图，选择能够连接node₁和node₂的路径。在深度优先遍历中，以contig₁为种子序列，用公式F(ss,ec)＝RM(ss,ec)^MC(ss,ec)对候选扩展序列打分，和公式进行罚分。并以contig₂为种子序列，以公式RM对候选序列打分。如果所有的得分值都大于0.3，则一条路径继续遍历，否则终止该路径。每次都以得分值最高的结点优先遍历，如果候选扩展序列为node₂，则遍历的这条路径作为结果填充该gap区域。整个填充空白区域过程结束。如果没有找到这样的路径，则以k_new＝k_new-2值重复上述步骤，直到k_new小于0或者小于k-5。

五、实验验证

为了验证本方法的有效性，我们在两组模拟数据和四组真实数据上进行了测试，并和目前流行的Velvet和PE-Assembly方法进行比较分析。为了评价拼接结果的连续性和准确性，我们采用N50，N90，contig个数，最长的contig长度，contig平均长度和覆盖度六个指标进行评价。当所有的最终contig按照长度由大到小进行排序时，N50是指中间contig的长度，该contig是所有长度大于等于它的其它contig之和大于等于所有contig长度之和一半的那条contig。N90是大于等于90％的contig长度。覆盖度是contig结果能够覆盖原始基因组(reference)碱基的百分比。

为了验证我们的方法在真实数据上的效果，我们在四组真实数据上进行验证。这四组真实数据是由Illumina技术测序得到，包括两组细菌Staphylococcus aureus(Staph)和Escherichia coli(Ecoli)，两组真菌Schizosaccharomyces pombe(Spombe)和Neurospora crassa(Ncra)。这四组数据均包括两组读数集合，每组的读数长度都比较小，覆盖度在200倍左右。

模拟数据集合在staph和ecoli两个物种上由ART软件产生，其数据特点如下表：

表1模拟数据

表2.真实数据信息

表.3staph拼接结果

Simulated Data(staph-length:2903107)

表.4ecoli拼接结果

Simulated Data(ecoli-length:4638902)

从表2和表3我们可以看出，在模拟数据上，我们方法的六个指标都比其他两个方法有着明显的优势。通过四组真实数据实验结果可以看出(表5-8)，我们的方法PEGA所取得contig个数是最少的，并且覆盖度是最高的，这说明我们的方法PEGA能够更好的克服掉一些复杂重复区和覆盖度不均衡问题。并且我们的方法PEGA在最大contig长度和N50上也领先于另外两种方法，除了在第四组数据上最大contig略短于PE-Assembly方法，这可能是由于第四组的reference并不完整造成的。

表5.Staph拼接结果

True Data(staph-length:2903107)

表6.Ecoli拼接结果

True Data(ecoli-length:4638902)

表.7Spombe拼接结果

True Data(spombe-length:12554318)

表.8Ncra拼接结果

True Data(ncra-length:39225835)

Claims (4)

1.一种基于读数和距离分布的基因组De novo序列拼接方法，其特征在于，包括以下步骤：

1)输入双端读数文库，构建初始De Bruijn图，并对初始De Bruijn图进行优化；对所述初始De Bruijn图进行优化的方法为：首先De Bruijn图中的简单路径合并为一个结点；简单路径是指一个路径中起始结点的出度为1，最终结点的入度为1，并且所有中间结点的出度和入度均为1；删除初始De Bruijn图中长度小于2*k，并且出度或入度为0的结点；

2)以De Bruijn图为基础，选择种子序列，并利用打分函数对候选扩展序列进行打分，选择得分最高的候选扩展序列与种子序列合并，并继续进行扩充，直到结束条件为止，扩展后的每个种子序列即为一条contig，所有的contig构成一个contig集合；contig集合构建的具体步骤为：

2A)将初始De Bruijn图中结点长度大于双端读数文库insertsize的结点为种子序列，以所述种子序列的前驱结点和后继结点为候选扩展序列；其中，insertsize为双端读数文库中每对读数之间的间距；

2B)如果候选扩展序列的长度小于minLen，则以该候选扩展序列为起点，采取深度优先遍历方法，抽取出所有长度大于minLen的子路径，每个子路径对应于一个二级候选扩展序列，对所有的二级候选扩展序列打分，并以最高分作为候选扩展序列的得分值；minLen取值为2*k；

2C)针对一个种子序列ss和一个候选扩展序列ec，采用公式

F(ss,ec)＝RM(ss,ec)^MC(ss,ec)

计算ec的得分值，其中RM是匹配上ss和ec的双端读数的相对值；由匹配上的双端读数得到一个距离集合，MC用于衡量该距离集合和整体insertsize分布的偏移大小；如果在ec中的实际存在读数个数不符合二项式分布B(minLen，p)，则将所述ec的得分值减去一个罚分值 其中a是罚分系数，b为ec中的k- mer平均频次与读数文库中的k-mer平均频次的比值；re是在ec中实际存在的读数个数，p是读数出现在读数文库中的概率；

2D)如果一个种子序列只有一个候选扩展序列，则将该候选扩展序列直接合并到种子序列中；如果一个种子序列存在多个候选扩展序列，并且得分最高的两个候选扩展序列得分值都大于0.7或者相差小于0.2，或者没有候选扩展序列，则扩展终止，否则把得分最高的候选扩展序列合并到种子序列中去；迭代地对种子序列进行左右扩展，直到达到上述的扩展终止条件，最终的种子序列作为一条contig加入到contig集合中；

2E)当扩展完所有的种子序列，则contig集合构建完成；

3)建立scaffold图，每个结点代表一个contig，边代表两条contig在真实基因组序列上的位置紧邻；

4)填充scaffold图中有边相连的两个结点之间的空白区域，通过匹配上的双端读数长生一个局部读数集合，在局部读数集合上构建新的De Bruijn图，并在该De Bruijn图上寻找能够连接两个结点的路径，如果存在这样的路径，则以该路径填充空白区域。

2.根据权利要求1所述的基于读数和距离分布的基因组De novo序列拼接方法，其特征在于，所述步骤1)中，构建初始De Bruijn图的过程如下：依次读取读数文库中的每条读数，每条读数分解成r-k+1个k-mer，在De Bruijn图中查找每个k-mer对应的结点，如果某k-mer对应的结点不存在，则在De Bruijn图中创造一个新的结点并对应于该k-mer；如果一个k-mer最后k-1个碱基和另一个k-mer前k-1个碱基相同，则在这两个k-mer对应的结点之间加一条有向边；当遍历完所有的读数时，初始De Bruijn图建立完成；其中r为每条读数的长度；k为k-mer的长度；且k为奇整数。

3.根据权利要求1所述的基于读数和距离分布的基因组De novo序列拼接方法，其特征在于，所述步骤3)中，确定contig集合中contig位置关系的过程为：

1)建立scaffold图，每个contig作为一个结点；

2)对于contig集合中的任两个contig：contig₁和contig₂，先由与contig₁和contig₂分别匹配的双端读数估计出gap的距离大小，即gapDistance；以contig₁最后insertsize-gapDistance个碱基为测试区域，把该测试区域分成若干个长度为minLen+r-1的子区域，以每一个子区域为候选扩展序列，并以contig₂为种子序列，利用公式对每个子区域打分，其中RN为子区域中存在双端读数的读数个数，SN为双端读数匹配上种子序列的读数个数；如果任何一个子区域的RM小于0.2，则认为contig₁和contig₂之间不存在相邻关系，contig₂作为contig₁右邻居的得分值为0；如果所有子区域的RM均大于0.2，则contig₂作为contig₁右邻居的得分值为所有子区域的平均RM值；当计算完contig₁和所有其它contig之间的得分值之后，选择得分最大的contig作为contig₁的右邻居，并在scaffold图中contig₁和contig₁的右邻居这两个结点之间增加一条边，同理计算contig₁左邻居得分值；

3)当完成所有contig的左右邻居的计算，并得到初始的scaffold图后，则删除错误的边，若contig₁的得分值最大的右邻居是contig₂，contig₂得分值最大的左邻居是contig₁，则认为contig₁和contig₂紧邻，两者之间的边保留；否则，删除它们之间的边。

4.根据权利要求3所述的基于读数和距离分布的基因组De novo序列拼接方法，其特征在于，所述步骤4)中，空白区域填充方法如下：

1)针对两个contig₁和contig₂，收集与所述contig₁和contig₂末端区域分别匹配的所有双端读数，构成一个局部读数集合；

2)以该局部读数集合和k_new构建一个De Bruijn图，并在该De Bruijn图中找到对应于contig₁右端的结点node₁，和对应于contig₂左端的结点node₂；k_new＝k+5，如果k+5大于r，则k_new取r-1；

3)从node₁出发，深度优先遍历该De Bruijn图，选择连接node₁和node₂的路径；在深度优先遍历中，以contig₁为种子序列，用公式F(ss,ec)＝RM(ss,ec)^MC(ss,ec)对候选扩展序列打分，利用公式 进行罚分；并以contig₂为种子序列，以公式RM对候选序列打分，如果所有的得分值都大于0.3，则一条路径继续遍历，否则终止该路径；每次都以得分值最高的结点优先遍历，如果候选扩展序列为node₂，则遍历的这条路径作为结果填充node₁和node₂之间的区域；整个填充空白区域过程结束；如果没有找到这样的路径，则令k_new＝k_new ^*，转到步骤2)，直到k_new小于0或者k_new小于k-5；其中k_new ^*＝k_new-2。