CN104850761B - 核酸序列拼接方法及装置 - Google Patents

Abstract

本发明公开了一种核酸序列拼接方法及装置，包括：接收测序序列，测序序列包括读段和测通数据；根据读段构建原始拼接图；将测通数据比对到原始拼接图的边上；从原始拼接图的边集中选择锚点边；构建以锚点边为中心的局部子图；化简局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。依据本发明的方法，通过将测通数据比对到由读段构建得的原始拼接图，从中选择出锚点边并据此构建局部子图，再通过对局部子图的化简和合并得到更长的路径，达到解决锚点边之间的重复序列的路径选择问题，进而可完成测序序列的拼接任务，为提升scaffold构建效果提供可能。

Description

核酸序列拼接方法及装置

技术领域

本发明涉及生物信息处理领域，尤其涉及一种核酸序列拼接方法及装置。

背景技术

在基因组测序技术发展历程中，具有低成本、高测序通量（又称高通量，指一定时间内获得的数据输出量）、更精准特点的第二代测序技术带来了前所未有的应用前景，如基因组组装、结构变异检测等。第二代测序技术得到的read（读段，指随机打乱DNA序列所得的小片段）不仅长度短，数量又极为巨大，覆盖度高，这给序列拼接带来了前所未有的难题，而基因组测序中的一个关键步骤就是序列拼接。序列拼接是将小的reads进行比对合并成大的DNA序列或者RNA序列，这些reads是由碱基A、G、C、T或A、G、C、U组成，reads的长度表示其含有的A、G、C、T或A、G、C、U的个数。含有A、G、C、T或A、G、C、U的reads经过拼接得到更长的reads，该更长的reads再经过拼接或比对最终可得到DNA或RNA序列。

第二代测序技术所面临的基本拼接问题是，由于第二代测序技术read长度偏短的特点，由此给拼接算法带来难题，包括如何处理测序错误，以及如何处理重复序列。在真实的基因组项目中，大基因组尤其是植物基因组的重复序列含量十分高，有的重复序列高达60%以上，在使用现有拼接算法进行序列拼接过程中，由于重复序列的影响，将产生很多短碎的contig（重叠群）序列，使得后续拼接过程的scaffold（支架，是由顺序确定的contigs连接组成）构建、补洞等环节受其影响，效果不佳。

已有研究提出使用长reads来处理拼接过程的重复序列。长reads可通过例如测通的Illumina数据得到。研究表明，长read有利于在contig构建过程中利用新的算法理论与新的拼接策略解决更多的重复序列的问题，得到更长的contig序列，提升contig构建的效果。

发明内容

依据本发明的一方面，提供一种核酸序列拼接方法，包括：接收测序序列，所述测序序列包括读段和测通数据；根据所述读段构建原始拼接图；将所述测通数据比对到所述原始拼接图的边上；从原始拼接图的边集中选择锚点边，所述锚点边的两端没有分叉且跨过该锚点边的读段的路径没有冲突；构建以所述锚点边为中心的局部子图；化简所述局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。

依据本发明的另一方面，提供一种核酸序列拼接装置，包括：接收模块，用于接收测序序列，所述测序序列包括读段和测通数据；原始构建模块，用于根据所述读段构建原始拼接图；比对模块，用于将所述测通数据比对到所述原始拼接图的边上；选择模块，用于从原始拼接图的边集中选择锚点边，所述锚点边的两端没有分叉且跨过该锚点边的读段的路径没有冲突；子图构建模块，用于构建以所述锚点边为中心的局部子图；化简模块，用于化简所述局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；合并模块，用于对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。

依据本发明的又一方面，提供一种包括上述核酸序列拼接装置的基因组测序设备。

依据本发明的再一方面，提供一种计算机可读存储介质，用于存储供计算机执行的程序，本领域普通技术人员可以理解，在执行该程序时，通过指令相关硬件可完成上述核酸序列拼接方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的方法，通过将测通数据比对到由读段构建得的原始拼接图，从原始拼接图中选择出锚点边并据此构建局部子图，再通过对局部子图的化简和合并得到更长的路径，达到解决锚点边之间的重复序列的路径选择问题，进而可完成测序序列的拼接任务，为提升scaffold构建效果提供可能。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解，其中：

图1为依据本发明的一种实施例的核酸序列拼接方法的流程示意图；

图2至图4为由于片段中间的碱基错误形成的环形结构和由测序错误形成的嵌合体示意图；

图5和图6分别为边出现分叉路径和路径存在冲突的示意图；

图7为通过锚点边D获得的所有局部子图及其化简结果示意图；

图8为由锚点边D构成的局部子图与由锚点边F构成的局部子图的合并示意图；

图9为局部子图合并过程中弧关系的处理示意图；

图10为依据本发明的一种实施例的核酸序列拼接方法的流程示意图；

图11为Maize chr1 contig组装结果示意图；

图12为Maize chr1 scaffold构建结果示意图；

图13为Maize chr1 contig评价结果示意图；

图14为Maize基因组补洞后scaftig统计结果示意图；

图15为Maize基因组scaftig评价结果示意图。

具体实施方式

目前已有研究者开发出不同的算法并编写大量的拼接软件。比较成功的拼接软件是velvet、SOAPdenovo等。其中，SOAPdenovo是深圳华大基因开发，是专门针对Illumina高通量测序数据的拼接软件，可从http://soap.genomics.org.cn/soapdenovo.html下载获得。SOAPdenovo基于de Bruijn（德布鲁因）图算法，其拼接流程一般包括如下步骤A-F。

步骤A，构建不同长度的插入读段文库，例如180bp（base pair，碱基对）、500bp等；

步骤B，将所有小读段库（180/500bp）测序得到的reads截断成更小的序列片段，通过它们之间的重叠关系构建de Bruijn图，从而将这些读段连接起来；

步骤C，在步骤B中构建的de Bruijn图非常复杂，为了简化de Bruijn图，需要去除无法继续连接的分支、低覆盖度的分支（一般这两种情况是由于测序错误造成）；并利用reads信息化简重复序列在de Bruijn图的分叉通路；对于少量的杂合位点，采用随机选择策略，合并杂合位点；

步骤D，通过步骤C得到简化后的de Bruijn图，这样的图仍然会有很多分叉位点无法确定真正的连接关系，因此在每个分叉位点将序列截断，得到了最初的contigs；

步骤E，将所有文库测序得到的reads比对回初步得到的contigs，利用reads之间的连接关系和插入读段大小信息，将contigs组装成scaffolds；

步骤F，对contigs之间的gap（空隙）进行补洞，延长contigs。补洞后的contigs长度相比补洞之前一般可以增加6-7倍。

其中，de Bruijn图的构造过程如下：

（1）假定给定的reads集合为F={f1,f2,…,fn}，把这些reads划分为若干连续碱基组成的等长度短序列（称为kmer）。假设每个kmer长度为k，对一个read划分时，先以read的任意一端为起始位置，截取k个碱基，再将起始位置向后移动一个碱基，再截取k个碱基，依此类推，直至截取得到kmer的尾部到达read的另一端，这些kmer组成了de Bruijn图的顶点。

（2）对于任意两个kmer如u和v，如果u的后k-1个碱基序列与v的前k-1个碱基序列相同，则建立一条由u指向v的有向边。

通过以上两步即可构造出一个de Bruijn图。

本发明一种实施例以SOAPdenovo拼接软件为基础，提出一种核酸序列拼接方法及装置，如图1所示，实施例的核酸序列拼接方法包括以下步骤S1-S13。

步骤S1：接收测序序列。

这里，所接收的测序序列包括reads和测通数据。reads通常是指通过高通量测序技术得到的长度较短的reads；测通数据包括长reads和连通数据，连通数据是指测通的Illumina双末端数据（paired-end reads）连接形成的长读段（又称连通read），即插入读长小于两条read长度之和的数据而形成的一条长度更长的read，或者说是利用插入片段长度小于所测得的两条末端读长之和的数据的由两个末端读段连接形成的长read；而测通数据包括的长reads是指长度大于插入片段的reads，例如通过Roche 454测序平台获得的reads。连通数据的获取可采用已有的相关技术实现，例如采用华大基因研发的COPE（可从http://sourceforge.net/projects/coperead/files/src/下载获得）软件实现对连通数据的获取。以下实施例中以测通数据为连通数据为例进行说明，其它实施例中，可以采用长reads与构建的原始拼接图进行比对。

步骤S3：根据读段构建原始拼接图（pregraph）。

这里，原始拼接图可基于de Bruijn图的构造方法得到，可以理解，得到的原始拼接图也可称为kmer图。

步骤S5：将连通数据比对（alignment）到原始拼接图的边上。

在执行本步骤时或者在本步骤之前，还进行去除处理，即对构建得到的原始拼接图，根据原始拼接图的拓扑结构，去除因碱基错误和/或测序错误形成的读段及该读段关联的所有边，例如错误信息出现时拼接得到的分支会非常短，需要去除该分支，从而保证了结果的准确性，同时无冗余存在。

在进行比对时，具体的比对方法可参考常见相关技术实现，这里以SOAPdenovo软件常用的第一步pregraph的比对为基础。一般地，为保证比对的可靠性，要求read和contig之间至少要有两个连续的共同kmer，即连通read和contig至少有（k+1）bp的重叠序列。此外，在本发明的一个实施例中，还可进行信息记录，即，对于一条边，记录跨过该条边的所有的reads的信息，对于一条reads ，记录该条reads所跨过的所有边的信息。

步骤S7：从原始拼接图的边集中选择锚点边。

实施例中本步骤具体包括如下步骤S701-S705。

步骤S701中，根据比对结果和拓扑结构，去除原始拼接图中因测序错误形成的读段及该读段关联的所有边。

具体地，实施例中根据拓扑结构，也就是根据reads的路径层次，进行测序错误去除处理，包括去除kmer图的测序错误以及弱连接分叉，例如图2至图4所示的由于reads中间的碱基错误所形成的泡状结构（bubble）和测序错误所形成的嵌合体，删除由于测序错误造成的read以及该read关联的所有边（即该read跨过的所有边）。

步骤S703中，根据原始拼接图的边集中各条边的边覆盖深度信息和边末端分叉信息，从原始拼接图的边集中确定出各边的类型。

具体地，边类型识别需要用到两个信息，一个是边覆盖深度信息，另一个是边末端分叉信息。首先，根据边覆盖深度信息，将边进行分类，一种举例中，以平均覆盖深度为基准，将边分成错误边（例如，边的覆盖深度小于等于0.1倍平均覆盖深度）、重复边（例如，边的覆盖深度大于1.8倍平均覆盖深度）、正常边（例如，边的覆盖深度在0.8倍平均覆盖深度到1.8倍平均覆盖深度之间）、杂合边（例如，边的覆盖深度在0.1倍平均覆盖深度到0.8倍平均覆盖深度之间）。然后，根据边两端分叉情况，将边进一步进行类别区分和确认。实施例中重点识别杂合边，根据杂合位点在kmer图形成bubble的特点，将组成bubble的两条边进行比对，将比对结果符合标准、且两条边的覆盖深度差不多且不大于平均覆盖深度（例如都是在平均覆盖深度一半左右）的bubble标志为杂合边。这里所说的标准是指，两路径碱基数量相差不超过第一阈值、错配碱基数量不超过第二阈值，所述第一、第二阈值根据经验值产生，比如在一个实施例中的数据，设置最大的读长为500bp，允许的不同碱基的数量在10个以内，并且两条序列之间的匹配碱基数量在90%或者以上，符合这个条件的两条边可以被看成是由于杂合所形成的两条路径。

步骤S705中，选择原始拼接图的边集中的非错误边和非重复边作为锚点边。

在步骤S705中，错误边和重复边直接被排出在候选之列，而锚点边的选择则是除去错误边和重复边的原始拼接图的边集中进行选择，即在正常边和杂合边中选择满足锚点边的条件。锚点边的选择条件首先应满足该边是唯一（unique）的，表现在拓扑图上就是这条边的两端没有任何的分叉路径，然后跨过该条边的reads路径应该是没有冲突的，如图5和图6所示，图5的情况是该条边的右端出现分叉路径，图6中，边C为待选的锚点边，通过待选边C的reads路径存在冲突的情况。

通常选择正常边作为锚点边；当然，也可以根据需要，决定是否将杂合边作为锚点边来进行局部子图构建，杂合边的定义如前述，即根据杂合位点在kmer图形成bubble的特点，将组成bubble的两条边进行比对，将比对结果符合标准、且两条边的覆盖深度都不大于平均覆盖深度（例如都在平均覆盖深度一半左右）的bubble标志为杂合边。

步骤S9：构建以锚点边为中心的局部子图。

这里，通过锚点边上的reads比对信息，获得由该条锚点边往左右两端延伸的局部子图，如图7所示，标号为D的带小黑点的长方块代表锚点边，通过该锚点边上的reads信息获得前后的边界（edge）信息，构成该锚点边的局部子图。

步骤S11：化简局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边。

本步骤中，通过去除支持数较低的分叉以及融合一些由同一锚点边获得的子图的方式化简局部子图。具体地，根据局部子图中各分叉的支持数和/或覆盖度，去除低于预设支持数和/或预设覆盖度的分叉，并根据该分叉对应的读段，去除包含该读段的其余局部子图，此外，对局部子图中的一条锚点边，以该条锚点边为中心，将与该条锚点边关联且不存在路径冲突的所有局部子图融合为一个局部子图。支持数或覆盖度通常是与基因组测序数据的实际情况紧密关联，支持数较低是指支持两条边的弧关系的跨过reads的数量，通常设定3条或3条以下为低支持数，同理，低覆盖度分叉一般也是由于测序深度或测序错误所导致的，通常将分叉长度小于两倍kmer长度的分叉当做低覆盖度分叉。

仍以图7为例，通过锚点边D的一堆reads获得的局部子图包括边B、C、D、E、F（即图示中从上往下的第一条局部子图），然后通过另外一些reads获得另外两个局部子图（即图示中从上往下的第二、三条局部子图），这三个局部子图都是由锚点边D获得的，而且路径上不存在冲突，因此，可以将这三个局部子图合并成图示箭头下方的局部子图，该合并得到的局部子图包括边A、B、C、D、E、F、G、H。

在化简局部子图的过程中还会出现由于测序错误导致的分叉的情况，根据支持该分叉路径的reads数量，例如支持该分叉路径的reads数量过少（例如支持数小于等于3）或者是支持该分叉路径的reads数所占比例小（例如小于0.2），可以将该分叉剔除，同时根据生成该分叉的reads序列，可以将其余包含该reads信息的局部子图的对应信息也删除，更新相应局部子图的信息，从而达到化简子图的目的。本实施例是从全局的角度处理长reads的测序错误造成的影响，即在去除当前局部子图中由于测序错误造成的分叉路径的同时，会将该read路径所关联到的所有局部子图的信息进行更新，凡是由该read路径所确定的局部子图都将被删除。

在化简结果中判断是否存在新的锚点边，如果是，则按该新的锚点边重复执行步骤S9，直至不存在新的锚点边。新的锚点边产生的原因（即触发这种迭代的条件）可以是由于化简后使得边的类型发生改变，例如杂合边变成正常边，或者原来含有多个分叉的边和reads路径存在冲突的边，由于去除测序错误的reads路径后，形成了前后只有一个分叉而且跨过该条边的reads路径没有冲突的符合锚点边条件的边，从而将这些符合条件的正常边或杂合边作为新的锚点边，以此循环进行处理。

步骤S13：对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。

本步骤中，如果某一unique序列代表的局部子图中含有另外一条unique序列，则可以把这两条unique序列代表的局部子图合并，又称路径合并。具体地，对处理后剩余的局部子图，判断其中的第一局部子图是否包含第二局部子图的一个锚点边，如果是，再判断第一局部子图的各边和第二局部子图的各边之间是否存在冲突，如果不存在，则合并第一局部子图和第二局部子图为最终拼接图。这里的“第一”和“第二”仅是为了叙述方便而给出的，并不是将局部子图真正划分为两种类型。如图8所示，在由锚点边D（图示中带小黑点的长方块表示锚点边）获得的局部子图（图示中从上往下第一条）中，其中一条边F刚好是另一个局部子图（图示中从上往下第二条）的锚点边，而且这两个局部子图的边之间没有冲突，因此这两个局部子图可以合并在一起，合并结果为图示中从上往下第三条局部子图。

然而，在进行局部子图的合并过程中，由于局部子图内部的边并不是独立存在的，其会跟其它局部子图外的边发生弧关系，因此在合并局部子图的时候，需要对这些弧关系进行处理，从而在最后由局部子图合并获得的路径是由第一个锚点到最后一个锚点之间的边组成。具体的处理过程如图9所示，包括：选择第一局部子图的第一个锚点边作为最终拼接图的边集的起点，选择第二局部子图的最后一个锚点边作为最终拼接图的边集的终点，将起点与终点之间的第一局部子图和第二局部子图中的边进行复制，并将复制的边加入最终拼接图的边集中，此外，还可保留第一个锚点边之前的边及其弧关系、以及最后一个锚点边之后的边及其弧关系，从而可保证后续分析时数据的完整性。在图9中，带有小黑点标记的长方形表示锚点边，无标记的表示正常边，带有下标r标记的表示重复边，这些边一般都有分叉，需要对其弧关系进行处理。在获得最后的锚点D到锚点L之间的边集的过程中，将其间遇到的重复边复制一份，将复制的内容加入到边集，而原来的边序列留在原处，原来的边弧关系依旧留在原处，重新建立新的边与前后边之间的弧关系。处理后的结果如图9中箭头下方所示，这样得到一条由锚点D开始、L结束的一条长路径，最终连成一条长contig。

综上步骤S1-S13，本实施例的优点在于，首先，在把连通的Illumina数据比对到原始边的过程中，允许一条read产生数条不连续的路径信息，提高数据的利用率；其次，对原始边进行了类型识别后，后续可以针对不同类型的边采取不同的处理措施；然后，根据路径信息去除测序错误，对锚点边进行筛选，提高可用锚点边的数量和准确性，其中锚点边的选择不仅限于单一序列，还包括重复序列里面的特异序列、杂合边等，在选取的方法上，仅通过测序深度排除测序覆盖度远大于平均覆盖度的边（不选为锚点边），而且去除测序错误是通过kmer图的拓扑结构和reads支持路径无冲突等条件来实现；最后，把具有重叠关系的不同锚点边的路径合并，获得更长的路径，达到了解决锚点边之间的重复序列的路径选择问题的目的；从而，在经过后续的其它处理后，可以获得更长contig序列。

图10示出了本实施例的核酸序列拼接方法的详细处理流程，图示中虚线框部分为SOAPdenovo原有的模块，实线框部分为本实施例发明新增的技术特征。可见，本实施例以特异性好的原始边为锚点边，利用连通的Illumina数据的长度优势，把不同的锚点边连接起来。这种策略对于锚点边之间的重复序列形成的拓扑结构没有要求，避免了SOAPdenovo现有策略中对重复序列形成的拓扑结构的要求作出的限制，增加了可以解决重复的重复序列的数量。此外，本实施例还通过迭代过程选择特异性较好的原始边作为锚点边，在每一轮的局部子图处理与更新（即去除测序错误的分叉，去除局部子图中reads测序错误指示的错误路径等），将再进行新一轮的锚点边选取，直至再也无法选取出新的满足条件（即前述锚点边的相关定义）的未被选取过的锚点边为止。

依据本发明的另一方面，还提供与前述核酸序列拼接方法相对应的拼接装置，其包括：

接收模块，用于接收测序序列，所述测序序列包括reads和连通数据；

原始构建模块，用于根据reads构建原始拼接图；

比对模块，用于将连通数据比对到原始拼接图的边上；

选择模块，用于从原始拼接图的边集中选择锚点边，锚点边的两端没有分叉且跨过该锚点边的读段的路径没有冲突；

子图构建模块，用于构建以锚点边为中心的局部子图；

化简模块，用于化简所述局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；

合并模块，用于对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。

上述模块的实现可参考前述方法实施例中相关的描述，不作重述。此外，一种实施例还提供一种采用该核酸序列拼接装置的基因组测序设备。

以下结合具体核酸序列对依据本发明的核酸序列拼接方法或装置的运行结果进行详细的描述。

在具体实例中，选取的是玉米1号染色体模拟数据，玉米1号染色体是高重复基因组序列，重复序列达到60%以上，在模拟数据中，首先模拟了40X插入片段长度在180左右，reads长度为100bp的数据，然后根据SOAPdenovo的特点，模拟了插入片段长度为500和800bp的数据各10X，2K，5K，10K数据各10X。在模拟数据生成以后，利用中国华大基因开发的COPE软件对插入长度为180的40X数据进行连接，用连接后的数据与其它数据一起进行组装。为了说明效果，还进行了与soapdenovo2.0结果的比照，即在相同数据与参数的前提下，分别用两个版本的程序进行组装，比较组装结果。

（1）contig组装结果与结果分析

首先，在相同参数和前期数据的情况下，用SOAPdenovo2.0（这里称为旧版本oldversion）与采用本发明实施例（这里称为新版本New version）分别构建contig，contig组装结果如图11所示的Maize chr1 contig组装结果，图示中参数K表示用于构建原始拼接图所使用的kmer的长度，-d表示去除出现频率小于d参数所设定的值的kmer，-R表示用reads跨过解决原始拼接图中由于重复序列导致的重复结构（该值在旧版本的程序中只能简单处理短reads能跨过的简单的结构，在新版本的程序中是指使用选取锚点构建子图的方式解决原始拼接图中由于重复序列导致的复杂结构），-M表示合并由于杂合形成的bubble结构（根据杂合率的不同，可以设置1、2、3三个档值），readLen表示读段长度（包含使用COPE等工具将测通读段连通后的最长读段长度）。

玉米是高重复基因组，从采用SOAPdenovo2.0版本来看，N50, N90都很短，新版本通过利用连通数据或是测通的长reads解重复处理后，contig N50、N90都有很大提升，N50都从200-300bp提升至2K以上。总体而言，对于玉米这样高重复的基因组，新版本在contig的N50、N90提升是较为明显的。

（2）scaffold组装结果与结果分析

对玉米1号染色体模拟数据，在相同参数和前期数据的情况下，用soapdenovo2.0与新版本分别构建scaffold，scaffold组装结果如图12所示的Maize chr1 scaffold构建结果，图示中Size_includeN表示包含有用N填充洞内序列的总长度，Size_withoutN表示剔除用N表示的序列后的总长度，Scaffold_Num表示支架个数，Mean_Size表示支架的平均长度，Longest_Seq表示最长支架长度，Singleton_Num表示支架结果序列中不含N的序列的个数，这种一般是没有参与构建支架的重叠群序列，Average_length_of_break(N)_in_scaffold表示将支架里面填充的N剔除后的序列平均长度。

从组装结果看，对于不同K值组装的scaffold在N50，N90，平均长度等方面，新版本相对于旧版本都有了很大的提升，K=63的结果中，scaffold N50从616996上升到1753939，N90也是从127bp提升到了119133bp。总的来说，新版本contig结果的改进也带来了scaffold的提升，最终获得更好的组装效果。

（3）组装结果评价与结果分析

首先，对Maize 基因组contig评价结果与分析。将组装后的 contig比对到reference上，获得组装后的contig的比对数量，基因组覆盖度等相关信息，如图13所示的Maize chr1 contig评价结果，图示中contig_num表示contig的个数，contig_size表示所有contig序列的长度总和，contig_coverage表示contig总长与基因组总长度的比值，aligned_num表示比对上参数序列的contig数目，aligned_size表示比对上参数序列的contig的总长度，aligned_coverage表示比对上的contig长度与基因组总长的比值，genome_covered_len表示contig序列覆盖上的基因组序列的长度，genome_coverage表示contig序列比对上基因组的长度与基因组总长的比值，mismatch_base表示比对上的contig的错配碱基数，mismatch_ratio表示错配碱基比例，indel_base表示插入缺失碱基数，indel_ratio表示插入缺失碱基比例。

从比对结果看，旧版本组装结果中，contig 数量更多，即短碎contig居多，最后在统计contig总长时候旧版本的contig总长更长，新版本为了使子图简化，以方便重复路径的处理，在contig前去除部分低覆盖度的分叉，这样操作的结果可能会误删除掉部分正确的路径，导致最后基因覆盖度略有下降，这可以借助后续的补洞把这部分补回来，见补洞后scaftig基因覆盖度（参见图14所示）。此外，新版本在indel与mismatch的数量与比例上有明显下降，这说明新版本相比旧版本准确度有所提升。

其次，对Maize 基因组scaffold补洞后scaftig（将补洞的N去掉后的序列）结果统计与评价分析，如图14所示的Maize基因组补洞后scaftig统计结果，图示中N50表示将序列按照长度从大到小排列后，将长度依照排序累加，当累加和达到大于或者等于总长度一般的那条序列的长度，N50_num表示长度累加到大于或者等于总长度一半时候所累加的序列的数量，N90表示将序列按照长度从大到小排列后，将长度依照排序累加，当累加和达到大于或者等于总长度0.9倍时候的那条序列的长度，N90_num表示长度累加到大于或者等于总长度0.9倍时候所累加的序列的数量，Total number(>=100bp)表示大于或者等于100bp的序列的总数量，Total length表示所有序列长度的总和，Total number(>=2kb)表示大于2kb的序列的数量，Average length表示序列的平均长度，NG50表示以基因组的真实或者估计大小作为总长所计算出来的N50，NG90表示以基因组的真实或者估计大小作为总长所计算出来的N90。

从上述统计结果看，补洞后，新版本的N50,N90相比旧版本提升了3-4倍，在数量上大大减少，说明新版本scaftig普遍较旧版本长，平均长度也较旧版本有较大提升，从一个更客观的数据，NG50及NG90来看，新版本的结果有较明显提升。

图15为Maize基因组scaftig评价结果，图示中含义类似图13。从上述统计结果看出，新版本在基因组覆盖度比旧版本高，从scaftig数量与总长可以看出，新版本整体上更长，完整性也更好。从准确性上看，新版本的mismatch数量与比例都比旧版本略小，indel数量也较旧版本低，因而新版本的准确性较好。

综上可知，本发明一种实施方式的设计思想是：以de Bruijn图为框架，在构建完成的kmer图基础上加入连通的Illumina数据；利用连通的reads长度的优势，将长reads比对到kmer图的边上，获得长reads跨过的边信息集，然后根据一定的条件设定在kmer图上选取符合条件的边作为锚点，以锚点为基准，通过跨过锚点边的reads比对到其它边的信息构建出以锚点边为中心的局部子图；通过对局部子图的去低覆盖度分叉、合并子图等操作，获得进一步化简、优化处理的子图，从而可以输出contig及进行后续其它处理。由此带来的有益效果是：1）更好解决SOAPdenovo组装过程中重复序列的组装难题，获得更长的contig；2）使用本发明产生的contig进行scaffold构建，获得更好的scaffold构建结果；3）改善补洞效果以及提高基因组的总体覆盖率。

本发明各实施例中处理的核酸序列可以是人造基因序列片段或者是通过基因测序仪测序得到的基因序列片段，该基因序列片段可以是DNA片段或RNA片段，本发明各实施例中对DNA片段和对RNA片段的处理方法无任何区别。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的另一方面还提供一种核酸序列拼接装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与上述数据输入单元、数据输出单元及存储单元数据连接，用于执行存储单元中存储的可执行的程序，该程序的执行包括完成上述实施方式中各种方法的全部或部分步骤。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换。

Claims (11)

1.一种核酸序列拼接方法，其特征在于，包括：

接收测序序列，所述测序序列包括读段和测通数据；

根据所述读段构建原始拼接图；

将所述测通数据比对到所述原始拼接图的边上；

从原始拼接图的边集中选择锚点边，所述锚点边的两端没有分叉且跨过该锚点边的读段的路径没有冲突；

构建以所述锚点边为中心的局部子图；

化简所述局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；

对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出；

从原始拼接图的边集中选择锚点边这一步骤包括：

根据比对结果和拓扑结构，去除原始拼接图中含有测序错误的读段及该读段关联的所有边；

根据原始拼接图的边集中各条边的边覆盖深度信息和边末端分叉信息，从原始拼接图的边集中确定出各边的类型；

选择原始拼接图的边集中的非错误边和非重复边作为锚点边。

2.如权利要求1所述的核酸序列拼接方法，其特征在于，将所述测通数据比对到所述原始拼接图的边上这一步骤或这一步骤之前包括：

对构建得到的原始拼接图，根据所述原始拼接图的拓扑结构，去除因碱基错误和/或测序错误形成的读段及该读段关联的所有边。

3.如权利要求1所述的核酸序列拼接方法，其特征在于，所述选择原始拼接图的边集中的非错误边和非重复边作为锚点边这一步骤包括：

选择原始拼接图的边集中的符合标准的杂合边作为锚点边，所述杂合边包括形成为泡状结构且组成泡状结构的两条边的覆盖深度相近且都不大于平均覆盖深度；所述符合标准指，将所述组成泡状结构的两条边进行比对，两条边的碱基数量相差不超过第一阈值且错配碱基数量不超过第二阈值。

4.如权利要求3所述的核酸序列拼接方法，其特征在于，所述组成泡状结构的两条边的覆盖深度为0.1-0.8倍的平均覆盖深度。

5.如权利要求3所述的核酸序列拼接方法，其特征在于，所述第一阈值为10，所述第二阈值为10％。

6.如权利要求1所述的核酸序列拼接方法，其特征在于，构建以所述锚点边为中心的局部子图这一步骤包括：

对每一条所述锚点边，通过跨过所述锚点边的读段的比对信息，得到所述锚点边向两端延伸的边界信息，根据所述边界信息获得所述锚点边的局部子图。

7.如权利要求1所述的核酸序列拼接方法，其特征在于，化简所述局部子图这一步骤包括：

根据局部子图中各分叉的支持数和/或覆盖度，去除低于预设支持数和/或预设覆盖度的分叉，并根据该分叉对应的读段，去除包含该读段的其余局部子图；

和/或，对局部子图中的一条锚点边，以该条锚点边为中心，将与该条锚点边关联且不存在路径冲突的所有局部子图融合为一个局部子图。

8.如权利要求1所述的核酸序列拼接方法，其特征在于，对处理后剩余的局部子图进行合并这一步骤包括：

对处理后剩余的局部子图，判断其中的第一局部子图是否包含第二局部子图的一个锚点边，

如果是，再判断所述第一局部子图的各边和所述第二局部子图的各边之间是否存在冲突，

如果不存在，则合并所述第一局部子图和所述第二局部子图为最终拼接图。

9.如权利要求8所述的核酸序列拼接方法，其特征在于，合并所述第一局部子图和所述第二局部子图为最终拼接图这一步骤包括：

选择所述第一局部子图的第一个锚点边作为最终拼接图的边集的起点，选择所述第二局部子图的最后一个锚点边作为最终拼接图的边集的终点，

将所述起点与所述终点之间的第一局部子图和第二局部子图中的边进行复制，并将复制的边加入最终拼接图的边集中。

10.一种核酸序列拼接装置，其特征在于，包括：

接收模块，用于接收测序序列，所述测序序列包括读段和测通数据；

原始构建模块，用于根据所述读段构建原始拼接图；

比对模块，用于将所述测通数据比对到所述原始拼接图的边上；

选择模块，用于从原始拼接图的边集中选择锚点边，所述锚点边的两端没有分叉且跨过该锚点边的读段的路径没有冲突；

子图构建模块，用于构建以所述锚点边为中心的局部子图；

化简模块，用于化简所述局部子图，在化简结果中重复选择锚点边进行处理直至不存在新的锚点边；

合并模块，用于对处理后剩余的局部子图进行合并，将合并结果作为拼接结果输出。

11.一种包括如权利要求10所述的核酸序列拼接装置的基因组测序设备。