提升基因组组装完整性的方法、装置及其应用

摘要

本发明提供了一种提升基因组组装完整性的方法、装置及其应用。该方法包括：获取目标样本的初步的染色体版本基因组；利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将最优比对短序列根据染色体进行聚类，得到多个类群；在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的组装基因组序列。通过首先利用测序序列，进行常规组装得到初组装基因组序列并进行染色体挂载之后，将三代短序列比对回上述的基因组上，将最优比对的三代短序列根据染色体聚成多个类群，之后在每个类群里做三代局部组装，从而获得完整性更高的基因组序列。

说明书

技术领域

本发明涉及基因组序列组装领域，具体而言，涉及一种提升基因组组装完整性的方法、装置及其应用。

背景技术

自1977年起，DNA测序技术先后经历了三个阶段。第一阶段主要是依据桑格尔(Sanger)和考尔森(Coulson)提出的双脱氧链末端终止法技术，也称为Sanger测序；第二阶段主要是以Roche公司的454测序平台、Illumina公司的Solexa测序系统为代表的第二代测序技术，也称为新一代测序技术("Next-generation" sequencing technology，NGS)；第三个阶段主要是以Pacific Biosciences公司的SMRT（single molecule real time）技术、Oxford Nanopore Technologies公司的纳米孔单分子测序技术，被认为是第三代测序技术。

伴随着测序的产生，科研工作者也开启了对物种基因组的探索。从小至几千碱基的噬菌体基因组，到三十亿碱基的人类基因组，再到最近公布的上百亿碱基的裸子植物（如加州红杉）和脊索动物（如非洲肺鱼）基因组，人们正在一点点的揭开基因组的神秘面纱。在这些研究的背后，基因组从头组装无疑是发挥了基础性的作用。并且随着研究的进一步深入，人们对组装基因组的连续性要求越来越高。在二代测序阶段，基因组组装的一般都比较碎，基因组N50通常在KB量级。后来随着三代长读长测序技术的产生，测序reads读长平均就可以达到几十到几百KB，从而使得组装基因组的连续性得到大幅度提高，重叠群 N50（Contig N50）基本上都在MB量级。并且借助Bionano、Hic等技术，可以将基因组挂载到染色体水平（重叠群之间用N连接，即空缺gap）。

最近，随着Pacbio推出HIFI测序模式（长度大于10kb，准确性大于99%），Nanopore推出Ultra-long模式（短序列 N50在60K以上），人们逐渐开始挑战基因组组装的终极目标，即整个基因组的染色体实现完整的从端粒到端粒组装（记为T2T组装），而没有任何空缺(gap)。然而，由于物种本身复杂的重复特性，以及组装算法的限制和组装策略的问题，目前除了人和水稻等少量物种实现了T2T组装之外，大部分物种基因组的空缺个数还依旧很多。

HiFi reads（High fidelity reads）是PacBio公司推出的基于CCS测序模式的兼具长读长和高准确度的测序序列，又称CCS序列。在这种测序模式下，因酶读长（平均~100Kb以上）远大于插入片段长度（~20 Kb），测序时，聚合酶会绕着DNA模板进行环形测序，使得插入片段被多次测序，产生多条subreads。之后来源于同一条模板链的subreads经过一致性校正，最终得到长度大于10KB，准确性大于99%的HiFi reads，用于基因组组装。

Nanopore Ultralong测序模式是Nanopore公司推出的可以有效提升测序readsN50长度的测序模式。它是新一代基于纳米孔的单分子实时电信号测序技术，其原理是：纳米孔蛋白作为生物传感器，插入聚合物形成的膜中。此外，核酸分子会与马达蛋白（MotorProtein）连接，该马达蛋白一方面对双链进行解链，使核酸单链在电泳的作用下通过特定的纳米孔蛋白，另一方面可控制DNA/RNA分子的移动速度，保证碱基逐一地穿过纳米孔，产生稳定可靠的电信号。由于不同碱基的带电性质不同，通过检测电信号的差异就能检测出通过纳米孔的碱基类别从而实现测序。通常，经过Nanopore Ultralong测序模式产生的reads， Reads N50能达到60-80K以上。

目前广泛报道的组装流程是首先搜索重叠读取以构建重叠群（contig）/支架(scaffold)，然后对重叠群/支架进行排序和定向以进行染色体级组装。长读长测序技术有助于获得更好和更可靠的组装，特别是通过防止同一区域中重复或重复序列引起的错误。比如物种基因组中端粒或者着丝粒区域都包含复杂的重复序列，这些重复序列的长度可能达到几百K甚至几M，属于远端重复，同时有些重复结构非常复杂，当前的三代测序技术是无法测通这些序列的；这些远端重复序列会对组装形成干扰，组装算法会误判远端相似序列与近端相似序列的链接关系，从而导致组装的结果序列断裂或者错误链接。

此外，已经发现通过将测序数据比对回已组装的基因组，会存在有些空缺区域能够被比对跨越，而这些区域富集了重复片段，这也证明了现在的组装方法的局限性。众所周知，在同一染色体上，Hi-C或遗传图谱技术具有比长距离序列更高的接触概率观察近端序列的能力，甚至相距> 200Mb。因此，Hi-C或遗传图谱主要应用于锚定 重叠群或支架到染色体规模的组件上。然而，Hi-C或遗传图谱也具有通过其相邻的非重复序列识别长距离重复序列的潜力。我们可以通过去除低频非邻近重复序列的读取链接噪声来优化装配流程以获得更完整的基因组装配。现在人们开始挑战基因组组装的终级目标，即实现染色体的T2T组装，然而目前大部分的基因组序列中还是包含很多的空缺。

因此，仍需要对现有的染色体组装方案进行改进，以进一步提高组装结果的完整性，达到无空缺组装的目标。

发明内容

本发明的主要目的在于提供一种提升基因组组装完整性的方法、装置及其应用，以解决现有技术中的大部分的基因组序列中仍含有很多空缺的问题。

为了实现上述目的，根据本发明的一个方面，提供了一种提升基因组组装完整性的方法，该方法包括：获取目标样本的初步的染色体版本基因组；利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将将最优比对短序列根据染色体进行聚类，得到多个类群；在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的组装基因组序列。

进一步地，获取目标样本的初步的染色体版本基因组包括：将三代测序短序列组装到重叠群水平，得到重叠群 V1；利用基因组图谱或者HIC数据，将重叠群V1挂载到染色体水平，得到初步染色体版本V1；利用juicerbox软件对初步染色体版本V1进行调整，得到初步的染色体版本基因组。

进一步地，利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将将最优比对短序列根据染色体进行聚类，得到多个类群包括：将三代测序短序列重新比对到初步的染色体版本基因组上，得到初比对文件；过滤掉初比对文件中重复比对的短序列，得到过滤后比对文件；根据过滤后比对文件，挑选出每一个染色体对应的短序列数据集，得到多个类群。

进一步地，在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的组装基因组序列包括：在每个类群中，利用三代测序短序列进行序列组装，得到类群组装结果；在每个类群组装结果中，利用Hi-C或遗传图谱的互作频率信息对重叠群进行初步挂载，得到初挂载结果；在每个类群中，将类群对应的三代测序短序列与初挂载结果进行比对，并根据比对结果对连接点做矫正和空缺填补，得到类群组装序列，类群组装序列即为完整性提升的组装基因组序列。

根据本申请的第二个方面，提供了一种提升基因组组装完整性的装置，该装置包括：获取模块，用于获取目标样本的初步的染色体版本基因组；比对聚类模块，用于利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将最优比对短序列根据染色体进行聚类，得到多个类群；重组装模块，用于在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的基因组序列。

进一步地，获取模块包括：第一组装单元，用于将三代测序短序列组装到重叠群水平，得到重叠群 V1；第一挂载单元，用于利用基因组图谱或者HIC数据，将重叠群V1挂载到染色体水平，得到初步染色体版本V1；第一调整单元，用于利用juicerbox软件对初步染色体版本V1进行调整，得到初步的染色体版本基因组。

进一步地，比对聚类模块包括：第一比对单元，用于将三代测序短序列重新比对到初步的染色体版本基因组上，得到初比对文件；过滤单元，用于过滤掉初比对文件中重复比对的短序列，得到过滤后比对文件，过滤后比对文件包括最优比对短序列；聚类单元，用于根据过滤后比对文件，挑选出每一个染色体对应的短序列数据集，得到多个类群。

进一步地，重组装模块包括：类群组装单元，用于在每个类群中，利用三代测序短序列进行序列组装，得到类群组装结果；第二挂载单元，用于在每个类群组装结果中，利用Hi-C或遗传图谱的互作频率信息对重叠群进行初步挂载，得到初挂载结果；填补单元，用于在每个类群中，将类群对应的三代测序短序列与初挂载结果进行比对，并根据比对结果对连接点做矫正和填补相应空缺处的序列，得到类群组装序列，类群组装序列即为完整性提升的组装基因组序列。

根据本发明的第三个方面，提供了一种计算机可读存储介质，该存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行上述任一种提升基因组组装完整性的方法。

根据本发明的第四个方面，提供了一种处理器，该处理器用于运行程序，其中，程序运行时执行上述任一种提升基因组组装完整性的方法的。

应用本发明的技术方案，通过首先利用测序序列，进行常规组装得到初组装基因之后，将三代短序列 比对（mapping）回染色体版本的初组装基因上（也可以用已发表的近源物种组装版本），将三代短序列根据染色体聚成多个类群，之后在每个类群里做三代局部组装，即可获得完整性更高的基因组序列。

附图说明

构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1示出了根据本发明的实施例1中的提升基因组组装完整性的方法的程示意图；

图2示出了根据本发明的实施例2中的提升基因组组装完整性的方法的详细流程示意图；

图3示出了根据本发明的实施例4中的提升基因组组装完整性的装置的结构示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合实施例来详细说明本发明。

术语解释：

Reads：高通量测序平台产生的短序列，或者PacBio单分子实时测序(singlemolecule real time，SMRT，包括CLR和HIFI数据)，ONT（Oxford Nanopore Technologies）测序产生的长序列，即为短序列。

序列组装（Sequence Assembly）：就是把基因组长的序列打断(shotgunsequencing)，由于不知道基因组整条序列是如何排列（成一条链，最后成为一条染色体）组合（如何区分不同染色体）的，而现有测序技术又无法实现一次把整条长序列完整测序，通过算法，计算机的帮助，把这些短的序列组装起来成为一条完整有序的序列。序列组装包括两种情况，第一种是有参考基因组序列，对重测序的序列进行组装，则可以用MAQ进行组装（即Map to reference genome）；第二种对新物种进行从头测序（de novo sequencing），则可以用相应的拼接软件进行组装（即De novo assembly）。

Contig：拼接软件基于短序列之间的overlap区，拼接获得的序列称为Contig（重叠群），无N。

Contig N50：将所有重叠群按照从长到短进行排序，并顺次相加，当相加的长度达到重叠群总长度的一半时，最后一个加上的重叠群的长度即为重叠群N50。重叠群N50 是基因组拼接结果好坏的一个判定标准。

Scaffold（支架）：基因组 de novo测序，通过短序列拼接获得重叠群后，根据确定的一些重叠群之间的顺序关系而组装成的更长的片段称为Scaffold（含N）。获得重叠群之后，通过利用二代paried-end或者mate-pair序列，三代短序列，10XG，bionano或HIC序列可以确定重叠群的顺序关系和位置关系，最后重叠群按照一定顺序和方向组成Scaffold，其中，形成Scaffold的过程中还需要填补重叠群之间的空缺（gap）。

Scaffold N50：与重叠群 N50的涵义类似。将所有Scaffolds按照从长到短进行排序，并顺次相加，当相加的长度达到Scaffold总长度的一半时，最后一个加上的Scaffold的长度即为Scaffold N50。Scaffold N50 是基因组拼接完整性的一个判定标准。

实施例1

本实施例提供了一种提升基因组组装完整性的方法，如图1所示，该方法包括：

S101，获取目标样本的初步的染色体版本基因组；

S102，利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将将最优比对短序列根据染色体进行聚类，得到多个类群；

S103，在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的组装基因组序列。

本申请的提升基因组组装完整性的方法，通过首先利用测序序列，进行常规组装得到初组装基因之后，将三代短序列比对（mapping）回染色体版本的初组装基因上（也可以用已发表的近源物种组装版本），将三代短序列根据染色体聚成多个类群，之后在每个类群里做三代局部组装，即可获得完整性更高的基因组序列。

上述步骤S101是先获得待提升完整性的目标样本的初步的染色体版本基因组，其可以按照常规的从头组装的方法，将二代或三代测序数据（优选三代测序数据）组装成有一定空缺的基因组。在一种优选的实施例中，S101包括：将三代测序短序列组装到重叠群水平，得到重叠群 V1；利用基因组图谱或者HIC数据，将重叠群V1挂载到染色体水平，得到初步染色体版本V1；利用juicerbox软件对初步染色体版本V1进行调整，得到初步的染色体版本基因组。利用juicerbox主要是为了调整重叠群的顺序和方向，提升组装结果的准确性，方便后续聚类。

为了进一步提升组装得到的基因组序列的连续性和准确性，本申请通过进一步将三代测序数据重新比对回初组装基因组（作为参考基因组进行比对）上，重新对染色体进行挂载的方式来实现更准确和更长的序列组装。

在一种优选的实施例中，S102包括：将三代测序短序列重新比对到初步的染色体版本基因组上，得到初比对文件；过滤掉初比对文件中重复比对的短序列，得到过滤后比对文件，过滤后比对文件包括最优比对短序列；根据过滤后比对文件，挑选出每一个染色体对应的短序列数据集，得到多个类群。

上述优选实施例通过将初组装基因组作为参考基因组序列，并将三代测序数据作为个体重测序数据进行再组装，通过将三代测序数据按照染色体将测序短序列进行分类聚集，得到分属于不同染色体的类群。这样比对有助于使组装序列更准确，将组装结果提升至最终的准染色体级别的基因组上。

更优选地的实施例中，通过对比对参数的设置进行优化，使得同一条短序列不会被分到多个类群中。在一些优选的实施例中，将三代测序数据比对回初步的染色体版本基因组上时，minimap2必须要添加参数--secondary=no（其作用是有效减少一条短序列有多个比对输出的概率）。过滤重复比对的短序列时，采用samtools软件过滤bam文件，添加-F2308参数（F2308参数为samtools软件中的固定参数，在samtools参数介绍中，-F INTonly include reads with none of the FLAGS in INT present [0]，即过滤掉-F后数值所包括的所有比对结果，此处2308=4+256+2048，则会将flag值为4、256和2048的比对结果都过滤掉），即过滤掉flag值为4（该read没比对到参考序列上），256（该read是次优的比对结果）和2048（补充匹配的read）的比对结果，或者使用samtools markdup去除掉重复的比对，这样可以有效避免后续将一条短序列聚类到多个类群中。

在一种优选的实施例中，S103包括：在每个类群里，利用三代测序短序列进行序列组装，得到类群组装结果；在每个类群组装结果中，使用Hi-C或遗传图谱的互作频率信息对重叠群进行初步挂载，得到初挂载结果；在每个类群中，将该类群对应的三代测序短序列与初挂载结果进行比对，并根据比对结果对连接点做矫正和空缺填补，得到类群组装序列，因此，得到了最终的完整性明显提升的组装基因组序列。

上述连接点指的是基于HIC或遗传图谱数据连接起来的重叠群序列的末端；矫正指的是根据比对文件，识别在基因组对应位置上，短序列碱基与基因组碱基的差异，比如基因组序列是A，短序列 序列中一半深度以上都支持是T，则将基因组该位置修改为T。

实施例2

本实施例以基于PacBio 平台的CCS序列的组装为例，结合图2详细描述组装过程。

1）利用CCS 短序列，基于hifiasm/hicanu等软件组装，得到重叠群 V1；

2）将HIC数据比对到重叠群 V1，然后利用ALLHIC软件中的extract，partition，optimize和build模块挂载到染色体水平，得到初步染色体版本pseudochromosome V1；

3）利用juicerbox软件对上述结果调整，得到pseudochromosome V2；

4）利用minimap2软件将三代短序列比对到pseudochromosome V2上，得到alignment bam；

5）利用samtools软件，按照flag值过滤，即samtools view -F2308（2308=4+256+2048），或者使用samtools markdup去除掉重复的比对，从而使每一个短序列只会对应一个比对结果，得到对应的filt.bam；

6）根据比对文件，挑选出每一个染色体对应的短序列数据集，认为是一个类群；

7）在每个类群里，利用hifiasm/hicanu组装；

8）在每个类群中，将HIC 短序列比对到类群的组装结果上，做染色体挂载，得到初挂载结果；

9）在每个类群中，将对应的三代短序列比对到上述8)中初挂载的染色体上，进行矫正和空缺填补。从而得到最终的完整性明显提升的类群组装序列（Final Assembly）。

基于Nanopore Ultralong组装的流程类似。

实施例3

本实施例对一种植物的序列组装进行了测试，对于流程前后对应的重叠群的结果对比发现，在不降低组装质量的前提下，重叠群N50由原来的14M提升到了19M，提升了34%（具体见下表）。

表1：

从以上的描述中，可以看出，本申请上述的实施例具有以下改进点：1）在将三代短序列比对到参考基因组上时，minimap2必须要添加参数--secondary=no，同时对于比对的bam文件，需要利用samtools软件，添加-F2308参数，即过滤掉flag值为4（该read没比对到参考序列上），256（该read是次优的比对结果）和2048（补充匹配的read）的比对结果，或者使用samtools markdup去除掉重复的比对，这样可以有效避免后续将一条短序列聚类到多个类群中；2）对于目前没有已发表参考基因组的物种，无需增加额外的测序数据，只依赖目前自身测序数据，也可以使基因组组装完整性得到有效提升；3）三代长短序列可以解决近端重复，但是对于远端重复（比如重复距离大于200M以上的）无能为力，本方法利用HIC数据，除了进行染色体挂载外，还可以解决部分远端重复序列，有效提升组装完整性。

本申请的实施例相比现有技术，具有以下优点：1）在不增加测序数据的前提下，能够有效提升组装完整性；2）使得原来包含空缺的染色体可以达到T2T组装。

需要说明的是，对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明并不受所描述的动作顺序的限制，因为依据本发明，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明所必须的。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到根据上述实施例的方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质（如ROM/RAM、磁碟、光盘）中，包括若干指令用以使得计算设备执行本发明各个实施例所述的方法，或者是使得处理器来执行本发明各个实施例所述的方法。

实施例4

本实施例提供了一种提升基因组组装完整性的装置，如图3所示，该装置包括：

获取模块10，用于获取目标样本的初步的染色体版本基因组；

比对聚类模块20，用于利用三代测序短序列与初步的染色体版本基因组的序列进行比对，将最优比对短序列根据染色体进行聚类，得到多个类群；

重组装模块30，用于在多个类群里对三代测序短序列的序列进行局部组装，从而获得完整性提升的基因组序列。

可选地，获取模块包括：第一组装单元，用于将三代测序短序列组装到重叠群水平，得到重叠群 V1；第一挂载单元，用于利用基因组图谱或者HIC数据，将重叠群V1挂载到染色体水平，得到初步染色体版本V1；第一调整单元，用于利用juicerbox软件对初步染色体版本V1进行调整，得到初步的染色体版本基因组。

可选地，比对聚类模块包括：第一比对单元，用于将三代测序短序列重新比对到初步的染色体版本基因组上，得到初比对文件；过滤单元，用于过滤掉初比对文件中重复比对的短序列，得到过滤后比对文件，过滤后比对文件包括最优比对短序列；聚类单元，用于根据过滤后比对文件，挑选出每一个染色体对应的短序列数据集，得到多个类群。

可选地，重组装模块包括：类群组装单元，用于在每个类群中，利用三代测序短序列进行序列组装，得到类群组装结果；第二挂载单元，用于在每个类群组装结果中，利用Hi-C或遗传图谱的互作频率信息对重叠群进行初步挂载，得到初挂载结果；填补单元，用于在每个类群中，将类群对应的三代测序短序列与初挂载结果进行比对，并根据比对结果对连接点做矫正和填补相应空缺处的序列，得到类群组装序列，类群组装序列即为完整性提升的组装基因组序列。

实施例5

本实施例提供了一种存储介质，该存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行上述任一种提升基因组组装完整性的方法。

此外，还提供了一种处理器，该处理器用于运行程序，其中，程序运行时执行上述任一种提升基因组组装完整性的方法。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的检测仪器等硬件设备的方式来实现。基于这样的理解，本申请的技术方案中数据处理的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本申请各个实施例或者实施例的某些部分的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

本申请可用于众多通用或专用的计算系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的消费电子设备、网络PC、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。

显然，本领域的技术人员应该明白，上述的本申请的部分模块或步骤可以在通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本申请不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。