CN107133493B - 基因组序列的组装方法、结构变异探测方法和相应的系统 - Google Patents

Abstract

本发明公开了一种基因组序列的组装方法、和相应的结构变异探测方法，该组装方法包括：通过设计序列映射的唯一性准则，将被测基因组的测序序列向参考基因组进行映射，并对映射结果进行恰当的切割，形成组装叠阵集。然后根据单映序列在组装叠阵集上的坐标和同伴关系估计基因组的构架，并根据组装叠阵集上单映序列的坐标和它们的同伴序列将叠阵向外延拓。延拓对各个叠阵以并行方式计算执行。延拓后的相邻叠阵一致序列如果存在重叠就将它们连接。这样所得到的当前组装基因组作为下一轮的参考基因组，通过调整序列映射的唯一性准则，重复以上拼接步骤，改进基因组的组装结果。所测基因组相对于参考基因组的结构变异，在拼接的过程中同时被探测出来。

Description

基因组序列的组装方法、结构变异探测方法和相应的系统

技术领域

本发明涉及生物信息技术领域，具体来说，涉及一种基因组序列的组装方法、结构变异探测方法和相应的系统。

背景技术

基因组测序是开展分子生物学研究的重要技术。通过对一个物种的基因组进行测序，研究人员可以获得这个物种的基因组碱基序列，它作为这个物种的遗传序列模版，为基因、转录、调控、修饰等层面进行定性或定量的研究，探索生命现象背后的分子机制提供了重要参照。完成测序后，通过将被测物种的基因组与其他物种的基因组进行比较，研究人员可以发现它们在基因组水平上的差异，这为揭示遗传变异、自然或人工选择的机制提供了信息，从而为优质基因的筛选、物种的改良培育提供了指导。此外，基因组测序还可以帮助寻找多倍体物种的杂合位点或杂合区段，是研究杂合性与生命现象的关系的重要基础。

第二代测序技术是目前应用的最广泛的测序技术。和第一代测序技术相比，它具有通量高、成本低的特点。在第二代测序技术发展的最初阶段，所测到的序列长度比较短(碱基对数目通常为30-40bp)而且碱基辨识的质量不够高；随着技术的不断改进，目前能够测量的序列长度大幅度增加(超过100bp)，同时碱基辨识质量也有了很大的改进。第二代测序的一个重要的特点是，它可以从两端对一个很长的片段进行测序，得到这个长片段两端的碱基序列，因此使用第二代测序技术可以获得高通量的双末端测序序列。

将测序序列组装成基因组是计算生物领域的基本问题。因为测序仪所能测量的序列长度远小于基因组长度，所以在测序后需要对所有测序序列进行组装，推断它们的相对位置，还原出被测的基因组。组装基因组面临着以下几个挑战：(1)第二代测序技术的数据具有很高的通量，大量的数据会增加组装的时间和计算设备上被占用的存储空间；(2)基因组上有很多相似度很高、或是重复出现的区段，它们的存在给推测测序序列的相对位置增加了很大的不确定性；(3)对于杂合度高的基因组，需要组装出不同的倍型，同时还要确定不同倍型的位置关系，找到杂合区域。

现有的基因组组装方法在原理上主要分为两类。一类是基于De Bruijn图的方法，该方法的主要操作是：对于每一个测序序列，每隔一个碱基切割出一个特定长度的子序列(通常称为k-mer，k表示子序列的长度)；利用所有被切割出来的子序列构造De Bruijn图；进行一定的纠错操作后，在图上寻找路径，每条路径被推断为被测基因组上的片段。这种方法对杂合度低，重复度低的基因组组装效果会比较好；而对于杂合度高，重复度高的基因组不是很理想。此外，基于De Bruijn图的方法不容易给出测序序列之间的叠落关系，不利于进行统计评估；同时，所切割出的子序列长度明显短于测序序列，会降低特异性，导致在图上寻找路径时出现错误。另一类方法是基于测序序列叠落关系的，该方法对每两条序列进行比对，根据全部比对结果推断序列的叠落关系。这类方法是以测序序列为单位的，而不是k-mer，容易从序列水平进行统计评估。然而，这一类方法需要对每两条测序序列进行比对，时间复杂度高，对测序通量小的第一代测序技术可以应用，但不适用于高通量的第二代测序技术。同时，对于重复度高的基因组，这类方法的效果也不一定理想，会出现拷贝数减少；对于杂合度比较高的区域，当基因组上同一位置的两个倍型差异比较大时，来自两个倍型的测序序列不容易被整合到一个叠阵中，导致倍型的丢失。

针对相关技术中的上述问题，目前尚未提出有效的解决方案。

发明内容

针对相关技术中的上述问题，本发明提出一种基因组序列的组装方法、变异探测方法和组装系统，能够实现了序序列的高效连接，实现测序序列的基因组装。

本发明的技术方案是这样实现的：

根据本发明的一个方面，提供了一种基因组序列的组装方法。

该组装方法包括：

(1)通过预定的映射算法将样品的被测基因组的测序序列映射到参考基因组，得到单映射测序序列叠阵集，其中，样品的测序序列为利用高通量测序技术测得，参考基因组已知并与样品的基因组相近；

(2)基于经过预处理的参考基因组对单映射测序序列叠阵集中的测序序列进行筛选，所得筛选结果根据覆盖度再次筛选，得到筛选后的单映射测序序列叠阵集；

(3)通过单方向测序序列信息对筛选后的单映射测序序列叠阵集进行切割，得到初始预组装叠阵集，将当前预组装叠阵集的初始值设置为初始预组装叠阵集；

(4)确定当前预组装叠阵集中每个叠阵的相对位置，形成组装基因组架构；

(5)对组装基因组架构中的每个叠阵进行延拓，得到每个叠阵的一致序列；

(6)将组装基因组架构中的相邻叠阵的一致序列中符合预定连接规则的一致序列进行连接，得到样品的当前的组装基因组；

(7)根据被测基因组上同源序列的差异调整预定的映射算法的映射参数，通过调整后的该预定的映射算法将样品的被测基因组的测序序列映射到当前的组装基因组，得到当前预组装叠阵集；

对当前预组装叠阵集迭代执行步骤(4)、(5)和(6)，迭代次数为任何非负整数。

其中，步骤(1)中的预定的映射算法中包括预定的映射参数，预定的映射参数包括以下至少之一：

被测基因组与参考基因组的差异预期；

被测基因组的长度、测序序列的长度和质量特征；

其中，映射参数用于提供判别任意一个测序序列和参考基因组上某个位置起始的子序列的相似度是否达到测序序列成功映射到参考基因组位置的准则。

此外，该步骤(1)包括：

在将样品的被测基因组的测序序列映射到参考基因组后，将被测基因组中映射到参考基因组上多个位置的测序序列去除，得到单映射测序序列叠阵集。

另外，该步骤(2)中对参考基因组的预处理包括：

对参考基因组进行自映射，得到参考基因组中的若干唯一性序列区域。

此外，在执行步骤(1)对测序序列进行映射时的映射率低于预定标准的情况下，则进行下述操作：

在执行步骤(2)后，对于筛选后的测序序列叠阵集，在每一个位置，选择最大频数的碱基，用最大频数的碱基更新参考基因组的唯一性序列区域上对应位置的碱基；

调整预定的映射算法的预定的映射参数，基于已经更新过唯一性序列区域的参考基因组，重新执行步骤(1)和步骤(2)。

其中，在执行步骤(1)的映射操作和步骤(2)的筛选操作时，如果被测基因组的双末端测序序列数据集的一对同伴序列的两端都被单映射到当前预组装叠阵集，则一对同伴序列的映射坐标信息在步骤(4)中用于形成组装基因组架构；

如果被测基因组的双末端测序序列数据集的一对同伴序列中的至少一端被单映射到当前预组装叠阵集，则一对同伴序列中的至少一端的映射坐标信息在步骤(5)中用于叠阵延拓；

其中，双末端测序序列数据集包括多个具有不同库长的片段库。

此外，该步骤(3)中的对筛选后的单映射测序序列叠阵集进行切割包括：

对于参考基因组上的每一个碱基，计算覆盖该碱基的所有左向测序序列的尾长的最大值W₁，以及计算覆盖该碱基的所有右向测序序列的尾长的最大值W₂；

如果W₁或者W₂小于一个预定的阈值w，则将该碱基标记为切割位点。预定的阈值w为整数，且0≤w≤L_max，L_max为所有测序序列长度的最大值；

其中，左向测序序列的尾长和右向测序序列的尾长的定义包括：

对于参考基因组上的任一个碱基，该碱基将覆盖该碱基的每一条测序序列分成左右两部分；其中，如果左侧部分的长度大于右侧部分的长度，则称该测序序列为左向测序序列，并且右侧部分的长度为该左向测序序列的尾长；如果右侧部分的长度大于左侧部分的长度，则称该测序序列为右向测序序列，并且左侧部分的长度为该右向测序序列的尾长。

另外，该步骤(4)包括：

利用测序序列的库长信息以及单映的测序序列在当前预组装叠阵集中的坐标，确定当前预组装叠阵集中任意两个叠阵间的距离范围；

对当前预组装叠阵集中的所有叠阵进行排列，使每两个叠阵间的距离与确定的对应该每两个叠阵间的距离范围相匹配。

此外，该步骤(5)包括：

在当前预组装叠阵集中的每个叠阵的每个端点附近设定一个范围，利用单映的测序序列在叠阵中的坐标信息，确定叠阵里范围内的测序序列的同伴序列，同伴序列与叠阵的一致序列共同构成从端点向外延拓叠阵的测序信息库，预定范围与叠阵中的测序序列所属的片段库的库长一致；

对测序信息库中的所有序列按照局部比对的算法进行比对，得到两两对比结果；

利用图论的深度优先算法整合两两比对结果，形成每个端点附近延拓后的叠阵；

基于延拓后的叠阵定义延拓的一致序列。

此外，该步骤(6)包括：

利用局部比对算法判断相邻叠阵一致序列是否存在重叠情况；

在存在重叠的情况下，将该相邻的叠阵的一致序列进行连接，得到样品的当前的组装基因组。

此外，该步骤(7)中的预定的映射算法中包括预定的映射参数，预定的映射参数包括以下至少之一：

被测基因组上同源序列的差异；

被测基因组的长度、测序序列的长度和质量特征；

映射参数用于提供判别任意一个测序序列和当前基因组上某个位置起始的子序列的相似度是否达到测序序列成功映射到参考基因组位置的准则。

该步骤(7)包括：

在将样品的被测基因组的测序序列映射到当前的组装基因组后，将被测基因组中映射到当前的组装基因组上多个位置的测序序列去除，得到当前单映射测序序列叠阵集；

对当前单映射测序序列叠阵集执行步骤(3)，得到当前预组装叠阵集。

此外，基于上述组装方法中的任意一种组装方法还可进行双倍体序列的组装。

根据本发明的另一方面，提供了一种应用上述组装方法中任意一种组装方法的结构变异探测方法。

该结构变异探测方法包括：

根据组装方法对不同样本的基因组之间的结构变异情况进行探测，探测的信息包括步骤(3)中对单映射测序序列叠阵集进行切割时形成的断点信息。

根据本发明的再一方面，提供了一种基因组序列的组装系统。

该组装系统包括：

映射模块，用于通过预定的映射算法将样品的被测基因组的测序序列映射到参考基因组，得到单映射测序序列叠阵集，其中，样品的测序序列为利用高通量测序技术测得，参考基因组已知并与样品的基因组相近；

筛选模块，用于基于经过预处理的参考基因组对单映射测序序列叠阵集中的测序序列进行筛选，所得筛选结果根据覆盖度再次筛选，得到筛选后的单映射测序序列叠阵集；

切割模块，用于通过单方向测序序列信息对筛选后的单映射测序序列叠阵集进行切割，得到初始预组装叠阵集；

架构模块，用于确定初始预组装叠阵集中每个叠阵的相对位置，形成组装基因组架构；

延拓模块，用于对组装基因组架构中的每个叠阵进行延拓，得到每个叠阵的一致序列；

连接模块，用于将组装基因组架构中的相邻叠阵的一致序列中符合预定连接规则的一致序列进行连接，得到样品的当前的组装基因组；

调整映射模块，用于根据被测基因组上同源序列的差异调整预定的映射算法的预定的映射参数，通过调整的该预定的映射算法将样品的被测基因组的测序序列映射到当前的组装基因组，得到当前预组装叠阵集；

延拓模块和连接模块进一步用于对调整映射模块中的当前预组装叠阵集进行操作。

本发明通过将被测基因组的测序序列和参考基因组进行映射，并对映射结果进行切割，以及将切割后的叠阵进行组装和延拓，从而实现了测序序列的高效连接，实现测序序列的基因组装。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是根据本发明实施例的基因组序列的组装方法的总体流程图；

图2是根据本发明实施例的基因组组装方法的详细流程图；

图3是根据本发明实施例的确定参考基因组上自映射唯一性区域的示意图；

图4是根据本发明实施例的根据唯一性条件进行筛选的示意图；

图5是根据本发明实施例的利用单方向测序序列信息确定叠阵切割位点方法的示意图；

图6是根据本发明实施例的利用双末端信息连接叠阵估计相邻叠阵距离的示意图；

图7是根据本发明实施例的收集用于延拓叠阵一端所需要的测序序列的方法示意图；

图8是根据本发明实施例的计算两个测序序列的最优位移的示意图；

图9是根据本发明实施例的对叠阵进行延拓的方法的示意图；

图10(a)～图10(d)是根据本发明实施例的连接相邻叠阵的一致序列示意图；

图11是根据本发明实施例的通过再映射和局部组装获取双倍体序列方法示意图。

图12是根据本发明实施例的基因组序列的组装系统的框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。

根据本发明的实施例，提供了一种基因组序列的组装方法。

为了便于理解本发明实施例的基因组序列的组装方法，下面先对本发明实施例的相关技术名词作出如下定义和解释：

基因组：包含一个生物体所有遗传信息的遗传序列，由四种碱基排列而成，通常表示为含有A，C，G，T四种字符的序列，每一种字符代表一种碱基。

双倍体：大多生物个体体细胞含有两套基因组，一套来自于父本，另一套来自于母本；这两套基因组基本一致，但也存在差异；存在差异的位置为杂合位点或者杂合区域；在杂合位点或者杂合区域，两套基因组序列不一致，具有两个倍型。

被测基因组：需要确定碱基序列的基因组。

测序序列：由测序平台输出的某一生物体基因组子序列的测量结果，由若干个代表四种碱基的字符(A，C，G，T)组成，代表基因组的一个子序列。测序序列的一端为5’端，另一端为3’端；通常在测序平台给出的文件中，测序序列的左端为5’端，右端为3’端。

bp：英语basepair的缩写。100bp指一个碱基序列有100个碱基。

第一代测序技术：最早出现的测序技术，利用若干生物化学技术，将被测基因组序列切割成多个子序列，并确定每个子序列碱基组成；测序序列长度可以到500bp左右，但产生的测序序列数据量较小。

第二代测序技术：也称作新一代测序技术，英文名称为Next GenerationSequencing，英文缩写为NGS；其特点是能并行地产生大量的双末端测序序列数据，每一条测序序列的长度大多为100bp-200bp。

双末端测序序列：第二代测序技术在进行操作时，将被测的序列切割成若干片段，片段的长度可小可大；之后在片段的两端进行测序，在一端得到一个测序序列，在另一端得到一个测序序列；如果片段过长，中间部分一般无法测到；所得到的两个测序序列对应于同一个片段，称为双末端测序序列。

片段库、库长：在测序平台给出的同一个双末端测序序列文件中，所有双末端测序序列所属的片段长度相近，测序平台会估计出一个平均值；在下文中，称测序平台给出的每一个双末端测序序列文件为片段库，称其估计出片段长度平均值为库长。

碱基误读：由于测序技术的误差使得测序序列上某些碱基不同于真实碱基。

质量值：在测序过程中，测序序列的每一个碱基都有被误读的可能性，质量值是这一可能性的反应。测序序列的每一个碱基都对应一个这样的质量值，质量值越高，被误读的可能性越低。

参考基因组：一个已完成测序或者组装，明确了碱基序列的基因组，在研究过程中被用作模板或者参照。

插入：指相对于参考基因组，测序序列的某两个相邻的碱基之间额外插入了一段碱基序列。

删失：指相对于参考基因组，测序序列丢失一段或若干段碱基序列。

INDEL：一个插入或者删失。

测序序列映射：对于一个测序序列，在参考基因组上寻找一个与其长度大体一致的子序列，该子序列与测序序列完全一样，或者存在差异，但差异在预先设定的标准之内。通常以碱基替换个数，INDEL长度来衡量二者的差异。

碱基替换：测序序列被映射到参考基因组上后某些碱基和参考基因组上对应的碱基不相同。

映射结果：对于一个测序序列，如果参考基因组上存在子序列，且子序列与测序序列的差异在设定的标准之内，则该测序序列映射成功。映射结果包括：(1)参考基因组上的子序列最左端碱基的坐标，该值作为测序序列在参考基因组上的映射坐标；(2)测序序列的映射方向，可以按5’至3’的方向映射，也可以按3’至5’的方向映射；(3)测序序列与参考基因组上子序列的比对信息，包括两者之间存在的碱基替换，INDEL情况。当参考基因组上存在多个与测序序列差异在设定标准之内的子序列时，该测序序列具有多个映射结果，每一个映射结果对应于参考基因组上的一个子序列。

基因组组装：利用测序序列还原被测基因组，组装后得到的基因组称为组装基因组。

测序序列叠阵(可以简称为叠阵)：英文术语为contig或layout；表示由若干测序序列形成的集合，在该集合当中，每两条测序序列之间的叠落关系或者相对距离是确定的。测序序列叠阵可以用含有A，C，G，T，Ф五种字符的矩阵表示；矩阵的某个元素为Ф表示矩阵的该位置上不存在碱基；矩阵的每一列对应于叠阵的每一个位点且每一列至少含有一个非Ф字符；矩阵的列数为叠阵的长度；矩阵的每一行代表一个测序序列，每行第一个非Ф字符所在列的列标为该行所代表的测序序列在叠阵中的坐标。

叠阵集：由叠阵组成的集合。

架构(或组装基因组架构)：英文术语为scaffold；它由一系列按照一定顺序排列后的叠阵组成，其中所有叠阵的相对位置关系已知，作为组装基因组的一个框架。

一致序列：英文术语为consensus；对于叠阵的每一列，根据其包含的所有碱基信息推断出一个碱基，并将所推断出的碱基连接后得到的序列；推断碱基的一个可行的方法为，对叠阵的每一列，在其包含的所有碱基中取频率最大者。

在清楚了相关技术名词的定义和解释后，在本发明的一个实施例中，为使得基因组组装能够高效地进行，同时提高具有高杂合度，或者高重复度的基因组的组装准确性，提出了一种基于高通量测序技术的基因组组装方法。

具体的，本发明的实施例中提出的基因组组装方法以双末端测序序列和一个参考基因组作为输入。输入的参考基因组可以是来自相近物种的已被测序的基因组；也可以是使用其他组装方法对测序序列进行组装，但效果不够理想、需要进行修正的组装基因组。其中，在一个实施例中，双末端测序序列来自被测基因组，双末端测序序列可以包含多个具有不同库长的片段库。鉴于被测基因组与参考基因组来自相同或者相近物种，二者具有一定的相似度，本发明借助参考基因组，利用测序序列映射工具，实现被测基因组的组装。如图1所示，本发明包括下述步骤：

步骤1、将测序序列映射至参考基因组，并经过后续操作得到初始预组装叠阵集；所述初始预组装叠阵集为一个叠阵集，它是步骤2进行组装的基础；

其中，在一个实施例中，可根据参考基因组、测序序列的特征，测序序列所属基因组与参考基因组差异大小的先验信息，以及对映射的灵敏度、特异度的预期，设计映射所用的参数；

其中，测序序列的特征指测序序列的长度分布情况和质量值分布情况；二者所属物种基因组差异大小指两个基因组之间出现SNP、INDEL的频率以及INDEL的长度分布情况；而对于库长不同，或者测序序列长度不同的片段库，可以设计不同的映射参数；

在一个实施例中，在完成映射后，将映射到参考基因组上多个位置的测序序列去除，获得单映射测序序列叠阵集；

另外，在一个实施例中，对参考基因组进行自映射，获得参考基因组上具有自映射唯一性的区域，并利用唯一性条件对所述单映射测序序列叠阵进行筛选，得到筛选后的叠阵集；如果映射率不够理想，可以根据映射率的具体情况调整映射参数重新进行映射，也可以基于筛选后得到的叠阵集对参考基因组进行更新，并重新执行步骤1；获得筛选的叠阵集之后，评估其中每个叠阵的连续性，在必要的位置进行切割，得到初始预组装叠阵集；

步骤2、基于步骤1所得初始预组装叠阵集进行基因组的组装，得到组装基因组；

利用双末端序列的库长信息估计初始预组装叠阵集中每两个叠阵的距离，重新排列叠阵使每两个叠阵的距离与估计值吻合，得到组装基因组架构；对于每个叠阵，将位于其两端的测序序列的同伴序列，与位于叠阵两端的测序序列进行比对，使叠阵向两端延拓，并推断一致序列；比对相邻的叠阵的一致序列，根据比对结果精确连接一致序列，得到组装基因组。

步骤3、将测序序列映射至组装基因组，获得当前预组装叠阵集，基于当前预组装叠阵集重新执行步骤2，实现迭代组装。步骤3可以执行，也可以不执行。

完成上述步骤1，步骤2和步骤3后，将组装基因组输出至标准格式的文件。

以下详细阐述本发明的方法和原理。图2详细地示出了本发明实施例的相关处理流程。

本发明实施例的步骤1根据参考基因组，双末端测序序列的具体情况，以及对二者差异的先验知识设计映射参数，将测序序列映射至参考基因组，获得相对于参考基因组的单映的测序序列叠阵集，并利用唯一性条件对叠阵进行筛选，得到筛选后的叠阵集，再对其中的叠阵进行必要的切割，得到初始预组装叠阵集。其具体包括以下步骤：

步骤11、设计映射参数。映射所用的主要为：

A:测序序列与参考基因组之间的错误匹配个数上界M。该参数由参考基因组和测序序列所属物种的差异率γ，测序序列的长度l，以及测序序列碱基误读的平均比率。一般地，可以将M设为测序序列长度的5％或6％。如果预测参考基因组和测序序列所属物种的差异较大，可以将M调整为测序序列长度的10％至15％。

B：可以被探测的INDEL的长度的最大值MAX_INDEL，该参数反映参考基因组和测序序列所属物种的差异；同时，该参数也会影响映射的速度。一般地，MAX_INDEL可以选择为5。

C：参数S，表示一个测序序列最多可以被映射到参考基因组上S个位置。如果一个测序序列被映射到了参考基因组上超过S个位置，即认为该测序序列映射失败。优选地，可以将S设为大于10的整数。

如果采用单个子序列-延拓的映射算法，还需设计以下参数：

D：完全匹配子序列的长度下限值k。该参数表示在映射时需要在测序序列上找到一个长度不小于k的子序列，同时在参考基因组上能够找到一个与之完全相同的子序列，以该子序列在参考基因组上的位置作为测序序列的初步定位。该参数的选取依赖于参考基因组和测序序列所属物种的差异率，测序序列的长度，测序序列碱基辨识的平均错误率，以及参考基因组的总长度。

E：每个测序序列的子序列搜寻数目上限值U。该参数设得越大，所寻找的子序列数目越多，但也会增加计算时间。优选地，可以将其设为20。

在设计参数时，可以对同一个片段库的所有测序序列设计共同的参数。此时，可以以所有测序序列长度的众数L代替l。对于测序序列长度分布，或者碱基辨识质量值分布具有明显差异的不同片段库，可以设置不同的参数。

步骤12、对于每一个片段库，利用步骤11所设计的参数，将该片段库的双末端测序序列映射至参考基因组，将映射到多个位置的测序序列从映射结果中去除，得到单映射测序序列叠阵集。如果映射率低于预定的标准，则进入步骤13，否则结束步骤1的操作。

步骤13、利用唯一性条件对单映测序序列叠阵进行筛选，得到筛选后的叠阵集。其具体包括以下步骤：

步骤131、对参考基因组进行自映射，以获得参考基因组上具有自映射唯一性的区段。对于长度分布具有明显差异的不同序列片段库，可以分别进行操作，得到不同自映射唯一性区域。步骤131具体包括以下步骤：

步骤1311、从参考基因组上每隔一个碱基截取一个长度为L的序列，将这些序列映射到参考基因组上。优选地，可以利用步骤11中设定的参数完成此步的映射。

步骤1312、对参考基因组上的每一个碱基，设定一个深度值，将所有碱基的深度值设定为0。遍历步骤1311中所有序列的映射结果，进行以下操作：如果一条序列映射到了参考基因组上唯一的位置，则将参考基因组上该序列所覆盖住的所有碱基的深度值加1；如果该序列被映射到了参考基因组上多于一个位置，则对于每一个成功映射的位置，该序列所覆盖的所有碱基的深度值加1。

步骤1313、记录参考基因组上所有具有自映射唯一性的区域。所述自映射唯一性的区域为参考基因组上的一个区间，满足以下条件：该区间中所有碱基的深度值均等于L；任意包含该区间的其它区间均含有深度值不等于L的碱基。

图3为确定参考基因组上自映射唯一性区域的方法示意图，如图3所示，完成自映射后，所有连续的且深度为L的碱基形成一个自映射唯一性区域。

步骤132、遍历步骤12中得到的所有单映的测序序列，按照和步骤1312相同的方法重新计算参考基因组上每一个碱基的深度值；取D_α为深度值分布的上α分位数，其中α为大于0小于0.5的数，优选地，以将其取为0.05；

步骤133、对于步骤12中得到的每一个单映的测序序列，检验其是否同时满足以下两个唯一性条件。如果其不能同时满足，则将其从所处的单映射测序序列叠阵中去除，最后得到一系列新的叠阵，称这些叠阵组成的叠阵集为“经过唯一性条件筛选后的单映射测序序列叠阵集”；

其中，第一唯一性条件：在参考基因组上，所述测序序列完全包含于具有自映射唯一性的区域；第二唯一性条件：步骤132中计算出的参考基因组上被所述测序序列覆盖的所有碱基的深度值小于D_α。

图4为根据唯一性条件进行筛选的示意图，如图4所示，不满足唯一性条件的测序序列从单映测序序列叠阵中被去除。

步骤14、如果步骤12的映射率、或者单映射的测序序列所占比率低于预定的标准，可以进行以下操作；两个操作可以全部进行，也可以择一进行：

操作一：调整映射参数；优选地，可以将测序序列与参考基因组之间的错误匹配个数上界M调大；返回步骤12；

操作二：设定频率值下界θ和深度值上界d，对于参考基因组上的每个碱基，基于步骤133得到的叠阵，计算该碱基的深度值以及覆盖住该碱基位点的A，C，G，T四种核苷酸的频率；如果该碱基的深度值大于d，且最大的频率值q_max大于θ，则将参考基因组上该碱基替换为q_max所对应的核苷酸，返回步骤11。

步骤15、利用单方向测序序列信息评估步骤133所得叠阵集中的叠阵的连续性，在必要位置对叠阵进行切割，得到初始预组装叠阵集。

在具体阐述步骤15的操作之前，引入下述定义：

任取参考基因组上的一个碱基，考查覆盖该碱基的每一条测序序列，碱基将该测序序列分成左右两部分。如果左侧部分的长度大于右侧，则称该测序序列为“左向测序序列”，同时称右侧部分的长度为该左向测序序列的尾长；如果右侧部分的长度大于左侧，则称该测序序列为“右向测序序列”，同时称左侧部分的长度为该右向测序序列的尾长。

步骤15具体包括以下步骤：

步骤151、对于参考基因组上的每一个碱基，计算覆盖该碱基的所有左向测序序列的尾长的最大值，记之为W₁；同时计算覆盖该碱基的所有右向测序序列的尾长的最大值，记之为W₂.如果W₁或者W₂小于一个预定的阈值w，则将所述碱基标记为切割位点。所述w为一整数，0≤w≤L_max，L_max为所有测序序列长度的最大值。

图5以参考基因组上的碱基b为例示出了W₁和W₂的计算方式。图中的星号代表每条测序序列的中点，如果该中点位于b的左侧，则该测序序列为一左向测序序列，反之则为一右向测序序列。

步骤152、在步骤151所得的所有切割位点处，对步骤133得到的叠阵集中的叠阵进行切割，得到切割后的叠阵集。具体的操作为：在参考基因组上从左至右扫描所有公共切割位点，对每一个切割位点，从步骤133得到的叠阵集中找到覆盖该切割位点的叠阵；将该叠阵分割为两个叠阵，其中一个包含该叠阵中所有映射至所述切割位点左侧的测序序列，另一个包含该叠阵中所有映射至所述切割位点右侧的测序序列；分割之后得到的每个叠阵中，测序序列的相对位置关系与其在进行分割之前的叠阵中的相对位置关系一致。

下文中，至步骤3之前，若无特殊说明，出现的所有“叠阵”二字均指步骤15所得叠阵集，即切割之后得到初始预组装叠阵集中的叠阵。

本发明的步骤2基于步骤1得到的初始预组装叠阵集实现基因组组装，得到组装基因组。其具体包括以下步骤：

步骤21、利用单映测序序列在叠阵上的坐标信息以及双末端测序序列的库长信息估计任意两个叠阵的距离，并将所有叠阵进行排列，使每两个叠阵的距离与估计值相吻合。优选地，两个叠阵X₁和X₂的距离可以采用步骤211和步骤212进行估计：

步骤211、对两个叠阵中的测序序列进行扫描，如果X₁中的一个测序序列R₁与X₂中的一个测序序列R₂为双末端测序序列，则计算R₁的左端到X₁的右端的距离，以及R₂的右端到X₂的左端的距离；用R₁和R₂所在的片段库的库长减去两个距离之和，得到一个差值，该差值作为X₁与X₂的距离的一个观测值。

步骤212、统计步骤211得到的X₁与X₂的距离的观测值个数，如果个数大于设定的下界，则取这些观测值的中位数作为X₁与X₂的距离的估计值。

完成每两个叠阵的距离估计之后，将所有叠阵进行排列，使得排列后的每两个叠阵的距离与估计值相吻合。

图6为步骤21所述方法的示意图。如图6所示，I，II和III为三个叠阵，每条虚线连接的两个箭头表示一对双末端测序序列。根据双末端信息，将叠阵I，II和III排列为I→III→II，并估计相邻叠阵之间的距离，得到组装基因组架构。

步骤22、采用基于测序序列两两比对的方法在左右两端延拓每个叠阵。首先阐述将一个叠阵向右延拓的具体步骤，向左延拓可以用类似的方法实现。将一个叠阵X向右延拓的具体操作为：

步骤221、基于单映射测序序列在叠阵上的坐标信息，收集将向右延拓叠阵X时所需的测序序列，计算每个测序序列对于X的先验坐标，具体包括以下步骤：

步骤2211、建立集合SET并将其初始化为空集。记LEN_X为叠阵X的长度。以X的最左端为坐标原点，计算X中所有测序序列的起始坐标。

步骤2212、遍历X中全体测序序列，对每一个测序序列进行以下操作：

记R为所述测序序列，记pos_R为其在X上的起始坐标，记insert_size_R为R所属片段库的库长。检验R是否满足以下两个条件：(1)LEN_X-t·insert_size_R≤pos_R，t为预先设定的不小于1的数；(2)R在X上从左至右为5’端至3’端。如果R同时满足上述两个条件，则将R的同伴序列R′的反向互补序列添加到SET中，以pos_R+insert_size_R-L_R′作为R′对于X的先验坐标，其中L_R′为R′的长度。。

步骤2213、将位于X右端的一部分序列添加到SET中，以这些序列在X上的坐标作为先验坐标。

所得集合SET保存所有向右延拓叠阵X所需的测序序列。

图7为步骤221的一个简单示意图。图中所有实线箭头表示靠近叠阵X右端的测序序列，其中方向向右的测序序列的同伴序列表示为虚线箭头，这些同伴序列即为所要收集用于延拓X的测序序列。

步骤222、设定两个测序序列重叠部分的碱基匹配数目下界match_bound(优选地将其设定为大于20的整数)；设定所述匹配数目与重叠部分长度的比值下界ratio_bound(优选地将其设定为大于0.9且小于1的数)；设定先验坐标之差的阈值pos_bound(该阈值为不大于测序序列长度的非负数)。对于SET中任意两个测序序列R₁和R₂，如果二者先验坐标之差的绝对值小于pos_bound，则计算R₁相对于R₂的最优位移，以及R₂相对于R₁的相对位移。

记R₁的长度为

记R₂的长度为

所述R₁相对于R₂的最优位移按照以下步骤计算：

步骤2221、设定三个变量s，max_ratio和optimal_shift(R₁，R₂)。将s初始化为

将max_ratio初始化为ratio_bound，将optimal_shift初始化为正无穷。

步骤2222、计算两个数值overlap(R₁，R₂，s)和match(R₁，R₂，s)，具体的计算方法如下：

当s＜0时，

当s≥0时，

和的定义方式如下：如果R₁[i-s]＝R₂[i]，则否则如果R₁[i]＝R₂[i+s],则

否则

步骤2223、如果match(R₁,R₂,s)≥match_bound,则计算比值

ratio＝match(R₁，R₂，s)/overlap(R₁，R₂，s)。

如果ratio≥max_ratio，则将max_ratio替换为ratio，将optimal_shift替换为s。

步骤2224、如果

则结束操作，否则用s+1代替s，返回步骤2222。

完成上述操作后得到的optimal_shift(R₁，R₂)即为测序序列R₁相对于R₂的最优位移，其可以为正无穷。

图8为计算最优位移的一个实施例。如图8所示，若将R₂相对于R₁向右平移5bp，二者的不匹配数目为15；若将R₂相对于R₁向右平移10bp，二者的不匹配数目为9；若将R₂相对于R₁向右平移8bp，二者的不匹配数目为0。因此R₂相对于R₁的最优位移为8，R₁相对于R₂的最优位移为-8。

步骤223、构建有向图G，G的每一个结点为集合SET中每一个测序序列；对于R中任意两条测序序列R₁和R₂，如果optimal_shift(R₁，R₂)为负，则在G中添加一条由R₁指向R₂的边；如果optimal_shift(R₁，R₂)为正且不等于正无穷，则在G中添加一条由R₂指向R₁的边。

步骤224、在步骤223构建的有向图G中，找出入度为0且先验坐标最小的结点。如果G中不存在入度为0的结点，则找出先验坐标最小的测序序列所对应的结点。以所述被找出的结点作为初始结点，进行以下操作：

步骤2241、对G中所有结点设定访问状态，将初始结点的访问状态设定为已访问，将其余结点的访问状态设定为未访问；

步骤2242、从初始结点起，按照深度优先原则对G进行遍历；在遍历过程中，对于每一个被访问到的结点，将其访问状态更改为已访问，并在该结点所指向的所有未被访问的结点中，选取相对于该结点最优位移最小的一个进行下一步的访问；如果该结点的出度为0，或者该结点所指向的所有结点均已被访问，则根据遍历过程记录从初始结点到该结点的路径；

步骤2243、在步骤2242记录的所有路径中，选出含结点数最多的一条。

步骤225、利用步骤2243选出的路径中所包含的测序序列构造一个叠阵，记这个叠阵为Y，具体的构造方式如下：

记所述路径所包含的结点依次为R₁，…，R_n，其中n为结点总数。将R₁在Y中的坐标设为1；对于任意的正整数i(2≤i≤n)，将R_i-1在Y中的坐标与optimal_shift(R_i，R_i-1)相加，作为R_i在Y中的坐标。

步骤226、整合X和Y，将X向右延拓，推断一致序列。其具体操作包括以下步骤：

步骤2261、如果Y中含有来自X的测序序列，则寻找一个同时属于X和Y的测序序列，记它在X和Y中的坐标分别为COOR_X和COOR_Y；对于Y中的每一个测序序列，将其在Y中的坐标与(COOR_X-COOR_Y)相加，作为其在延拓后的叠阵X中的坐标；对于同时包含于集合SET以及X中的测序序列，如果其不在Y中，则将其从X中去掉；如果Y中不包含来自X的测序序列，则以X自身作为延拓后的叠阵。

步骤2262、推断X的一致序列，记作C(X)。

上述步骤221至步骤226为向右延拓叠阵X，推断得到一条或两条一致序列的方法。向左延拓X时，可以将X整体取反向互补，之后同样进行步骤221至步骤226的操作，最后将所得到的一条或者两条一致序列取反向互补。

在具体实施过程中，可以并行地对每一个叠阵执行步骤22的操作。

步骤23、连接相邻两个叠阵的一致序列，得到组装基因组。

首先阐述两条序列seq1和seq2，且seq1在左，seq2在右的连接方法。一个优选的方案如下：

截取seq1不同长度的后缀序列，以及seq2不同长度的前缀序列。如果seq1的某个后缀序列suffix_seq1，和seq2的某个前缀序列prefix_seq2能够完全匹配，或者仅存在很少的碱基替换或者拆入/删失，则认为seq1和seq2能够成功连接，否则认为二者不能成功连接。如果seq1和seq2能够成功连接，则将seq2中位于prefix_seq2右侧的部分接到seq1右端，形成连接后的序列，记prefix_seq2的长度为二者公共部分的长度；如果seq1和seq2不能成功连接，则在seq1右端接上若干字符N，将seq2接在这些N的右端，得到将seq1和seq2连接后的序列。seq1和seq2连接后的序列记作seq1οseq2。

从左至右顺次扫描组装基因组架构，对每两个相邻、且已被延拓的叠阵的一致序列按照上述方案进行连接，得到组装基因组。

步骤3、如果步骤2所得组装基因组不够理想，可以以其作为参考基因组，将测序序列进行映射，获得当前预组装叠阵集，并基于当前预组装叠阵集执行步骤2，实现迭代组装，从而提高参考基因组的评估指标。其具体包括以下步骤：

步骤31、根据同源序列的差异和测序错误率设计映射所用的参数，具体的方法和步骤11相同；特别地，可以将测序序列与参考基因组之间的错误匹配个数上界M调小；

步骤32、按照步骤31设计的参数，将测序序列映射至组装基因组，将所有映射到多个位置的测序序列去除，得到单映测序序列形成的叠阵集；

步骤33、通过下述两种方法之一获得当前预组装叠阵集：

方法一：按照和步骤15相同的方法切割相对于当前组装基因组的单映测序序列叠阵；

方法二：找到组装基因组上相邻两个连续为N的片段，截取位于它们中间的不包含N的碱基片段，以该碱基片段和单映射至该碱基片段的测序序列作为需要被延拓的叠阵，这些叠阵构成当前预组装叠阵集；

步骤34、基于当前预组装叠阵集，利用单映射测序序列在当前预组装叠阵集上的坐标信息，执行步骤2实现迭代组装。

在经过一定次数的迭代组装之后，输出所得组装基因组，作为方法的最终输出。

以上为本发明的基本步骤。如按下述方案进行修改，本发明在组装基因组的同时，还可以构建出组装基因组上高杂合区域的双倍体序列，并且同时输出组装基因组、双倍体序列、双倍体序列与组装基因组的位置关系信息。所述修改方案为：以下述步骤a代替步骤226；以下述步骤b代替步骤23；在步骤3之后增加步骤c。

步骤a、整合X和Y，将X向右延拓，推断一致序列。完成此步骤的操作后可以得到一条一致序列，也可以得到两条一致序列。如果得到两条一致序列，每一条对应于一个倍型。图9示出了步骤a的步骤以及主要符号的含义。步骤a具体操作包括以下步骤：

如果Y中含有来自X的测序序列，则进入步骤a1，否则进入步骤a6。

步骤a1、寻找一个同时属于X和Y的测序序列，记它在X和Y中的坐标为COOR_X和COOR_Y；对于Y中的每一个测序序列，将其在Y中的坐标与(COOR_X-COOR_Y)相加，作为其在延拓后的叠阵X中的坐标；对于同时包含于集合SET以及X中的测序序列，如果其不在Y中，则将其从X中去掉；

步骤a2、推断X的一致序列，记作C(X)；

步骤a3、将Y中的测序序列从G中去除，如果G还有剩余测序序列，且其中至少有一条与X中的某个测序序列存在有向边，则进入步骤a4；否则认为X经延拓后得到一条一致序列，结束步骤a的操作；

步骤a4、利用G中剩余的测序序列构造一个叠阵H_Xr；记这些测序序列为r₁，...，r_m，一个优选的构造方式如下：

步骤a41、设定整数K₁，...，K_m，将其初始化为0；将H_Xr设定为空集；

步骤a42、遍历测序序列r₁，...，r_m，对于测序序列r_i(1≤i≤m)，若r_i与X中的一条测序序列v之间存在一条有向边，则将K_i更新为optimal_shift(r_i，v)与v在X中的坐标之和，并将r_i加入H_Xr中，将其在H_Xr中的坐标初始化为1；否则若r_i与已存在于H_Xr中的一条测序序列r_j之间存在一条有向边，则将K_i更新为K_j+optimal_shift(r_i，r_j)，并将r_i加入H_Xr中，将其在H_Xr中的坐标初始化为1；

步骤a43、重复执行步骤a42直到H_Xr中不再加入新的测序序列；

步骤a44、将K₁，...，K_m中未被更新的值去掉，在余下的整数中找到最小值，记作K_min；对于测序序列r_i(1≤i≤m)，如果r_i位于K_Xr中，则将r_i在H_Xr中的坐标更新为K_i-K_min+1；

步骤a5、推断H_Xr的一致序列，记作C(H_Xr)；记录数值K_min，作为H_Xr左端相对于X的平移量，也作为C(X)和C(H_Xr)的位置关系信息；将C(X)的第K_min个碱基标记为分叉点，记之以符号P_Xr；结束步骤a的操作。

如果Y中无来自X的测序序列，则执行步骤a6：

步骤a6、分别推断X和Y的一致序列；如果Y最左端的测序序列相对于X的先验坐标值大于X的长度，则将Y的一致序列接在X的一致序列右边，中间以若干个字符N相隔，形成一条一致序列，将其记为C(X)；否则保留两条一致序列，将X的一致序列记为C(X)，将Y的一致序列记为C(H_Xr)，同时记录位于Y最左端的测序序列的先验坐标值prior_{Y_left}，作为两条一致序列的位置关系信息，将C(X)的第prior_{Y_left}个碱基标记为分叉点，记之以符号P_Xr。

步骤b、按照步骤24中描述的连接两条序列的方案对相邻的每两个叠阵的一致序列进行连接，得到组装基因组、连接区域的双倍体序列，以及二者的位置关系信息，并将其输出。

以下结合图10(a)至图10(d)，以两个相邻的叠阵A和B为例，阐述一致序列的连接，以及三项结果的输出方法。以下分情况进行阐述：

图10(a)所示，A右端和B左端均存在双倍型：判断[C(A)]_r和[C(B)]_l、[C(A)]_r和C(H_Bl)、C(H_Ar)和[C(B)]_l、C(H_Ar)和C(H_Bl)能否成功连接；如果存在可以成功连接的序列对，则将连接部分长度最大的序列对进行连接，并将余下的两个序列进行连接；如果不存在可以成功连接的序列对，则任取一个序列对用N连接，将余下的两个序列用N连接；如图10(a)所示，被连接得到的两条序列分别为[C(A)]_rο[C(B)]_l和C(H_Ar)οC(H_Bl)；；连续输出[C(A)]_c、[C(A)]_rο[C(B)]_l和[C(B)]_c的每一个碱基作为组装基因组；输出C(H_Ar)οC(H_Bl)作为双倍体序列；计算P_Ar和P_Bl在组装基因组上的坐标并输出，作为双倍体序列与组装基因组的位置关系信息。

图10(b)所示，A右端存在双倍型，B左端不存在双倍型：连接[C(A)]_r和[C(B)]_c得到[C(A)]_rο[C(B)]_c(或者连接C(H_Ar)和[C(B)]_c得到C(H_Ar)ο[C(B)]_c)，将[C(A)]_rο[C(B)]_c(或者C(H_Ar)ο[C(B)]_c)输出作为组装基因组；输出C(H_Ar)(或者[C(A)]_r)作为双倍体序列；计算P_Ar在组装基因组上的坐标并输出，作为双倍体序列与组装基因组的位置关系信息。

图10(c)所示，A右端不存在双倍型，B左端存在双倍型：连接[C(A)]_c和[C(B)]_l得到[C(A)]_cο[C(B)]_l(或者连接[C(A)]_c和C(H_Bl)得到[C(A)]_cοC(H_Bl))，将[C(A)]_cο[C(B)]_l(或者[C(A)]_cοC(H_Bl))输出作为组装基因组；输出C(H_Bl)(或者[C(B)]_l)作为双倍体序列；计算P_Bl在组装基因组上的坐标并输出，作为双倍体序列与组装基因组的位置关系信息。

图10(d)所示，A右端和B左端均不存在双倍型：连接[C(A)]_c和[C(B)]_c得到[C(A)]_cο[C(B)]_c，输出[C(A)]_cο[C(B)]_c作为组装基因组。

从左至右顺次扫描组装基因组架构，对每两个相邻、且已被延拓的叠阵的一致序列进行上述操作，即可得到组装基因组，双倍体序列，以及二者的位置关系信息。

步骤c、利用映射和基于序列两两比对的组装方法构建非连接处的双倍体序列，作为步骤b所构建的双倍体序列的重要补充。图11为步骤c的操作示意图，主要操作包括测序序列再映射、局部化组装、推断局部化组装所得叠阵的一致序列得到双倍体序列。步骤c具体包括以下步骤：

步骤c1、将全体测序序列映射至步骤2得到的组装基因组，和所有已构建出的双倍体序列。在映射所用的参数中，错误匹配数目上界不应设得过大，可以设为测序序列长度的5％至6％。

步骤c2、对于步骤c1映射失败的序列，再映射至组装基因组；进行再映射时，将错误匹配数目的上界适度调大。如果进行再映射后，映射失败的序列比例仍然很大，可以将错误匹配数目的上界继续调大，对映射失败的序列重复再映射。

步骤c3、对于任意一条在步骤c2中映射成功的测序序列r，进行以下操作，以决定保留或者舍弃其映射结果：

步骤c31、如果r的同伴序列r′在步骤c1或者步骤c2中被成功映射至组装基因组(而非步骤b构建的双倍体序列)，则进入步骤c31；否则舍弃r的映射结果，结束步骤c3的操作。

步骤c32、记组装基因组上r被成功映射到的位置为pos_i(1≤i≤S)，s表示r被成功映射到的位置数目；记组装基因组上r′被成功映射到的位置数目为pos′_j(1≤j≤S′)，S′表示r′被成功映射到的位置数目。如果存在唯一的地址对满足：

(1)在和

这两个位置上，r与r′方向相反；

(2)

和

二者距离与r所属片段库的库长相近；

则保留r在

的映射结果，舍弃余下S-1个映射结果；否则将r的全部S个映射结果舍弃。

步骤c4、对经过步骤c3的筛选后被保留的测序序列按照在组装基因组上的坐标进行分类。一个可行的分类方法如下：将所有被保留的测序序列按照其在组装基因组上的坐标从小到大进行排序；将坐标最小的测序序列分入第一类；按坐标从小到大的顺序依次扫描排好序的所有测序序列，如果一个测序序列的坐标与上一个被扫描的测序序列坐标之差小于设定的阈值，则将其与上一个被扫描的测序序列分入同一类，否则将其分入新的类中。

步骤c5、对于步骤c4得到的每一类，将其中所有测序序列按步骤232、步骤233、步骤234、步骤235所述方法构建一个叠阵。构建叠阵时所用到的先验序关系可以根据测序序列在组装基因组上的坐标获得。推断所构建的叠阵的一致序列，作为一条双倍体序列；记录位于所构建叠阵最左端与最右端的测序序列在组装基因组上的坐标，作为该双倍体序列与组装基因组的位置关系信息。

执行完步骤c后，输出以下内容，作为方法的最终输出包括：组装基因组；步骤b和步骤c构建的所有双倍体序列；步骤b和步骤c构建的每一条双倍体序列与组装基因组的位置关系信息。

本发明提出的方法适用基因组组装。就基因组而言，本方法可以用于高杂合度、高重复度的基因组组装；就测序序列而言，本方法适用于高通量、双末端测序序列的组装。如果存在与被测基因组相近的参考基因组，本方法可以通过映射、局部化组装等步骤直接实现基因组组装，以及双倍体序列的构建；在操作时，可以选择一个参考基因组进行上述操作，也可以选择多个参考基因组，利用每一个参考基因组进行上述操作，最后将结果整合，得到组装基因组。如果不能找到合适的参考基因组，而且使用其他组装方法得到的组装基因组效果不理想时，可以将其他方法给出的组装基因组作为参考基因组，利用本方法对被测基因组进行重新组装，进行矫正。

本发明提出的方法还适用于不同基因组之间的比较。通过对叠阵进行重新排列以及后续的局部化组装，可以获得被测基因组与参考基因组在结构层面的差异信息，找到大片段的结构变异，以及两个物种基因组相似、或者变异密集的区段。

目前人类基因组已经有若干版本的参考基因组，作为一个特别的应用，本发明提出的方法可以基于人类参考基因组实现一个人的基因组组装，并探测其与参考基因组之间存在的结构层面的差异。这在个体化医疗中，对于实现基因组层面的疾病预测、诊断，可以起到重要的作用。

本发明提出的上述方法，具有以下优点：

1、在设计思路上，本方法并不等同地看待所有测序序列，一次性完成组装，而是先组装出唯一性强、可靠性高的叠阵，通过这些叠阵获得组装基因组架构，将重复区域和未组装的部分局部化，之后进行局部化的组装，并构建双倍体序列；这样有助于降低重复区域或者高杂合度给组装带来的不确定性；

2、通过将测序序列映射至参考基因组，根据映射结果得到一部分测序序列的叠落关系，用于实现被测基因组上与参考基因组相似度较高的地方的组装，对组装操作实现了部分程度的简化；

3、映射的参数设计有对应的定量评估方法来指导，可以根据对被测基因组与参考基因组的差异预期、测序序列的长度和质量特征，定量地设计不同的参数，以实现灵敏度与特异度的平衡；

4、在映射后，仅保留映射到参考基因组上一个位置的测序序列，减少测序序列叠落关系中的不确定性；

5、测序序列映射、映射结果的过滤，对于不同测序序列都是独立的，可以并行地实现；

6、在估计叠阵之间的距离时，同时利用不同片段库的测序序列，增加估计时所用的样本量；

7、利用参考基因组上具有自映射唯一性的区域过滤映射成功的测序序列，进一步提高特异度，减少叠阵间的错误连接；

8、排列叠阵后，被测基因组上的重复区域以及未被组装的区域被局部化，表现为相邻叠阵之间的空缺部分；已存在的叠阵可靠性较高，利用这些叠阵中的测序序列和双末端信息，寻找属于空缺部分的测序序列，使测序序列和空缺部分能够更加准确地对应；

9、经过局部化后，在延拓每一个叠阵时，所需组装的测序序列数目大幅度减少，使得进行延拓可以采用基于测序序列叠落关系的组装方法进行，从而保留了测序序列整体信息，无需将其切为k-mer；同时可以并行地对每一个叠阵进行延拓，增加方案执行的效率；

10、延拓叠阵、连接一致序列后，可以得到双倍体序列，以及双倍体和组装基因组之间的位置关系信息；

11、构建非连接处的双倍体序列时，基于再映射和双末端信息寻找属于另一个倍型的测序序列；在映射时允许的错误匹配数目较多，双末端信息的应用有助于纠正由此导致的错误映射；另外，找到的属于另一个倍型的测序序列具有在组装基因组上的再映射坐标，双倍体序列与组装基因组的位置关系可以较为准确地获得；

12、延拓叠阵，连接一致序列可以并行地实现；寻找非连接处的双倍体序列，也可以将测序序列根据再映射坐标进行分类，并行地实现。

根据本发明的实施例，还提供了一种应用上述组装方法中任意一种组装方法的变异探测方法。

根据本发明实施例的变异探测方法包括：

根据组装方法对不同样本的基因组之间的结构变异情况进行探测，探测的信息包括对单映射测序序列叠阵集进行切割时形成的断点信息。

根据本发明的实施例，还提供了一种基因组序列的组装系统。

如图12所示，根据本发明实施例的组装系统包括：

映射模块121，用于通过预定的映射算法将样品的被测基因组的测序序列映射到参考基因组，得到单映射测序序列叠阵集，其中，样品的测序序列为利用高通量测序技术测得，参考基因组已知并与样品的基因组相近；

筛选模块122，用于基于经过预处理的参考基因组对单映射测序序列叠阵集中的测序序列进行筛选，所得筛选结果根据覆盖度再次筛选，得到筛选后的单映射测序序列叠阵集；

切割模块123，用于通过单方向测序序列信息对筛选后的单映射测序序列叠阵集进行切割，得到初始预组装叠阵集；

架构模块124，用于确定初始预组装叠阵集中每个叠阵的相对位置，形成组装基因组架构；

延拓模块125，用于对组装基因组架构中的每个叠阵进行延拓，得到每个叠阵的一致序列；

连接模块126，用于将组装基因组架构中的相邻叠阵的一致序列中符合预定连接规则的一致序列进行连接，得到样品的当前的组装基因组；

调整映射模块127，用于根据被测基因组上同源序列的差异调整预定的映射算法的预定的映射参数，通过调整的该预定的映射算法将样品的被测基因组的测序序列映射到当前的组装基因组，得到当前预组装叠阵集；

架构模块124，延拓模块125和连接模块126进一步用于对调整映射模块127中的当前预组装叠阵集进行操作。

综上所述，借助于本发明的上述技术方案，通过将被测基因组的测序序列和参考基因组进行映射，并对映射结果进行切割，以及将切割后的叠阵进行组装和延拓，从而实现了测序序列的高效连接，实现测序序列的基因组装。

总之，本发明公开了一种基因组序列的组装方法、和相应的结构变异探测方法，以及基因组组装系统，该组装方法包括：通过设计序列映射的唯一性准则，将被测基因组的测序序列和参考基因组进行映射，并对映射结果进行恰当的切割，形成预组装叠阵集。然后根据单映射序列在组装叠阵集上的坐标和同伴关系估计基因组的架构，并根据组装叠阵集上单映序列的坐标和它们的同伴序列将叠阵向外延拓。延拓采用以下三部曲算法：1.序列两两比对；2.用图论方法整合两两比对结果形成延拓后的叠阵；3.基于所述延拓后的叠阵定义延拓的一致序列。上述延拓对各个叠阵以并行方式计算执行。延拓后的相邻叠阵一致序列经过比对判别，如果存在重叠，就将它们连接，从而完成一轮的基因组拼接。所得到的当前组装基因组作为下一轮的参考基因组，通过调整序列映射的唯一性准则，重复以上拼接步骤，改进基因组的组装结果。所测基因组相对于参考基因组的结构变异，在拼接的过程中同时被探测出来。

这项知识产权的研发得到了中国科学院B类先导专项“动物复杂性状的进化解析与调控"课题XDB13040600的资助，得到了国家自然科学基金委员会重大研究计划培育项目91530105、91130008的资助，以及中国科学院国家数学与交叉科学中心的各种支持。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (15)

1.一种基因组序列的组装方法，其特征在于，包括：

(1)通过预定的映射算法将样品的被测基因组的测序序列映射到参考基因组，得到单映射测序序列叠阵集，其中，所述样品的测序序列为利用高通量测序技术测得，所述参考基因组已知并与所述样品的基因组相近；

(2)基于经过预处理的参考基因组对所述单映射测序序列叠阵集中的测序序列进行筛选，所得筛选结果根据覆盖度再次筛选，得到筛选后的单映射测序序列叠阵集；

(3)通过单方向测序序列信息对所述筛选后的单映射测序序列叠阵集进行切割，得到初始预组装叠阵集，将当前预组装叠阵集的初始值设置为所述初始预组装叠阵集；

(4)确定所述当前预组装叠阵集中每个叠阵的相对位置，形成组装基因组架构；

(5)对所述组装基因组架构中的每个叠阵进行延拓，得到每个叠阵的一致序列；

(6)将所述组装基因组架构中的相邻叠阵的一致序列中符合预定连接规则的一致序列进行连接，得到所述样品的当前的组装基因组；

(7)根据所述被测基因组上同源序列的差异调整所述预定的映射算法的映射参数，通过调整后的该预定的映射算法将所述样品的被测基因组的测序序列映射到所述当前的组装基因组，得到当前预组装叠阵集；

对所述当前预组装叠阵集迭代执行所述步骤(4)、(5)和(6)，迭代次数为任何非负整数。

2.根据权利要求1所述的组装方法，其特征在于，所述步骤(1)中的所述预定的映射算法中包括预定的映射参数，所述预定的映射参数包括以下至少之一：

所述被测基因组与所述参考基因组的差异预期；

所述被测基因组的长度、测序序列的长度和质量特征；

其中，所述映射参数用于提供判别任意一个测序序列是否能够成功映射到参考基因组上某个位置起始的子序列的准则。

3.根据权利要求1所述的组装方法，其特征在于，所述步骤(1)包括：

在将样品的被测基因组的测序序列映射到参考基因组后，将所述被测基因组中映射到所述参考基因组上多个位置的测序序列去除，得到所述单映射测序序列叠阵集。

4.根据权利要求1所述的组装方法，其特征在于，所述步骤(2)中对所述参考基因组的预处理包括：

对所述参考基因组进行自映射，得到所述参考基因组中的若干唯一性序列区域。

5.根据权利要求4所述的组装方法，其特征在于，在执行权利要求1步骤(1)对所述测序序列进行映射时的映射率低于预定标准的情况下，则进行下述操作：

在执行权利要求1步骤(2)后，对于筛选后的所述测序序列叠阵集，在每一个位置，选择最大频数的碱基，用所述最大频数的碱基更新参考基因组的唯一性序列区域上对应位置的碱基；

调整所述预定的映射算法的预定的映射参数，基于已经更新过唯一性序列区域的参考基因组，重新执行所述权利要求1步骤(1)和所述步骤(2)。

6.根据权利要求1所述的组装方法，其特征在于，在执行所述步骤(1)的所述映射操作和所述步骤(2)的筛选操作时，如果被测基因组的双末端测序序列数据集的一对同伴序列的两端都被单映射到所述当前预组装叠阵集，则所述一对同伴序列的映射坐标信息在所述步骤(4)中用于形成所述组装基因组架构；

如果所述被测基因组的双末端测序序列数据集的所述一对同伴序列中的至少一端被单映射到当前预组装叠阵集，则所述一对同伴序列中的所述至少一端的映射坐标信息在所述步骤(5)中用于叠阵延拓；

其中，所述双末端测序序列数据集包括多个具有不同库长的片段库。

7.根据权利要求1所述的组装方法，其特征在于，所述步骤(3)中的对所述筛选后的单映射测序序列叠阵集进行切割包括：

对于所述参考基因组上的每一个碱基，计算覆盖该碱基的所有左向测序序列的尾长的最大值W₁，以及计算覆盖该碱基的所有右向测序序列的尾长的最大值W₂；

如果W₁或者W₂小于一个预定的阈值w，则将该碱基标记为切割位点，所述预定的阈值w为整数，且0≤w≤L_max，L_max为所有测序序列长度的最大值；

其中，所述左向测序序列的尾长和所述右向测序序列的尾长的定义包括：

对于所述参考基因组上的任一个碱基，该碱基将覆盖该碱基的每一条测序序列分成左右两部分；其中，如果左侧部分的长度大于右侧部分的长度，则称该测序序列为左向测序序列，并且所述右侧部分的长度为该左向测序序列的尾长；如果右侧部分的长度大于左侧部分的长度，则称该测序序列为右向测序序列，并且所述左侧部分的长度为该右向测序序列的尾长。

8.根据权利要求1所述的组装方法，其特征在于，所述步骤(4)包括：

利用所述测序序列的库长信息以及单映的测序序列在当前预组装叠阵集中的坐标，确定所述当前预组装叠阵集中任意两个叠阵间的距离范围；

对所述当前预组装叠阵集中的所有叠阵进行排列，使每两个叠阵间的距离与所述确定的对应两个叠阵间的距离范围相匹配。

9.根据权利要求1所述的组装方法，其特征在于，所述步骤(5)包括：

在当前预组装叠阵集中的每个叠阵的每个端点附近设定一个预定范围，利用单映的测序序列在所述叠阵中的坐标信息，确定所述叠阵里所述范围内的测序序列的同伴序列，所述同伴序列与所述叠阵的一致序列共同构成从所述端点向外延拓所述叠阵的测序信息库，所述预定范围与所述叠阵中的测序序列所属的片段库的库长一致；

对所述测序信息库中的所有序列按照局部比对的算法进行比对，得到两两对比结果；

利用图论的深度优先算法整合所述两两比对结果，形成所述每个端点附近延拓后的叠阵；

基于所述延拓后的叠阵定义延拓的一致序列。

10.根据权利要求1所述的组装方法，其特征在于，所述步骤(6)包括：

利用局部比对算法判断相邻叠阵一致序列是否存在重叠情况；

在存在重叠的情况下，将该相邻的叠阵的一致序列进行连接，得到所述样品的当前的组装基因组。

11.根据权利要求1所述的组装方法，其特征在于，所述步骤(7)中的所述预定的映射算法中包括预定的映射参数，所述预定的映射参数包括以下至少之一：

所述被测基因组上同源序列的差异；

所述被测基因组的长度、测序序列的长度和质量特征；

所述映射参数用于提供判别任意一个测序序列是否能够成功映射到当前基因组上某个位置起始的子序列的准则。

12.根据权利要求1所述的组装方法，其特征在于，所述步骤(7)包括：

在将所述样品的被测基因组的测序序列映射到所述当前的组装基因组后，将所述被测基因组中映射到所述当前的组装基因组上多个位置的测序序列去除，得到当前单映射测序序列叠阵集；

对所述当前单映射测序序列叠阵集执行所述步骤(3)，得到所述当前预组装叠阵集。

13.基于权利要求1～12任意一项所述的组装方法的双倍体序列的组装方法。

14.一种应用如权利要求1～12任意一项所述的组装方法的变异探测方法，其特征在于，包括：

根据所述组装方法对不同样本的基因组之间的结构变异情况进行探测，探测的信息包括所述权利要求1步骤(3)中对所述单映射测序序列叠阵集进行切割时形成的断点信息。

15.一种基因组序列的组装系统，其特征在于，包括：

映射模块，用于通过预定的映射算法将样品的被测基因组的测序序列映射到参考基因组，得到单映射测序序列叠阵集，其中，所述样品的测序序列为利用高通量测序技术测得，所述参考基因组已知并与所述样品的基因组相近；

筛选模块，用于基于经过预处理的参考基因组对所述单映射测序序列叠阵集中的测序序列进行筛选，所得筛选结果根据覆盖度再次筛选，得到筛选后的单映射测序序列叠阵集；

切割模块，用于通过单方向测序序列信息对所述筛选后的单映射测序序列叠阵集进行切割，得到初始预组装叠阵集；

架构模块，用于确定所述初始预组装叠阵集中每个叠阵的相对位置，形成组装基因组架构；

延拓模块，用于对所述组装基因组架构中的每个叠阵进行延拓，得到每个叠阵的一致序列；

连接模块，用于将所述组装基因组架构中的相邻叠阵的一致序列中符合预定连接规则的一致序列进行连接，得到所述样品的当前的组装基因组；

调整映射模块，用于根据所述被测基因组上同源序列的差异调整所述预定的映射算法的预定的映射参数，通过调整后的该预定的映射算法将所述样品的被测基因组的测序序列映射到所述当前的组装基因组，得到当前预组装叠阵集；

所述架构模块、所述延拓模块和所述连接模块进一步用于对所述调整映射模块中的当前预组装叠阵集进行操作。

Images