CN112133371B - 基于单管长片段测序数据进行骨架组装的方法和装置 - Google Patents

Abstract

一种基于单管长片段测序数据进行骨架组装的方法和装置，该方法包括：获取stLFR测序数据，其包括条形码信息和读长对信息；对stLFR测序数据进行预组装处理，获得种子重叠群；利用条形码信息，确定种子重叠群之间的相对顺序；利用条形码信息，确定近邻种子重叠群之间不同取向的权重；利用读长对信息，将未进行组装的重叠群进行挂载；利用条形码信息和读长对信息特征估计重叠群之间的相对间距大小；输出最终组装得到的骨架序列。本发明引入了stLFR数据所包含的长程信息，即利用来自同一条长片段的读长带有相同条形码信息以及来自同一个短片段的读长对信息来综合判断重叠群的排序和取向及它们之间的间距大小，以连接重叠群得到比较理想的骨架长度。

Description

基于单管长片段测序数据进行骨架组装的方法和装置

技术领域

本发明涉及二代测序技术领域，具体涉及一种基于单管长片段测序数据进行骨架组装的方法和装置。

背景技术

随着高通量以及高精度二代测序技术的成熟与普及，生物信息学技术已经成为系统深入研究和开发各类生物遗传资源的重要手段。其中高质量的基因组序列是保证生物信息下游分析准确性以及完整性的关键。由于二代测序技术读长(read)较短(一般在100bp～300bp之间)，用这些短的读长无法拼接出完整的基因组序列，特别是对于高杂合和高重复率的物种。因此，如何获取基因组的长片段读长(LFR)信息已成为目前基因组测序技术以及生物信息研究领域的热点问题之一。

最近来，深圳华大生命科学研究院自主研发了一种新的条形码(barcode)标记技术，即单管长片段测序技术(Single tube Long Fragment Read，以下简称stLFR)，其参考文献为Wang,O.,et al.Efficient and unique co-barcoding of second-generationsequencing reads from long DNA molecules enabling cost effective and accuratesequencing haplotyping and de novo assembly,Genome Research,2019.。该技术方案，首先将长片段DNA均匀分散并吸附到含有特异条形码的磁珠，然后在转座酶的作用下，磁珠表面的长片段DNA被打断成小片段DNA，同时将磁珠上特异的条形码标记到小片段DNA上，最后使用二代测序仪对标记了条形码的小片段进行测序，产生的数据读长为100碱基对(basepair，bp)，数据具有同二代测序技术一样的高准确率。现有的针对长片段读长数据的组装技术方案主要有两类。第一类是启发式的延伸策略：该策略首先利重叠群(contig)之间共享的条形码信息挑选出局域的候选组装重叠群，然后再利用启发式算法确定候选组装集合中重叠群之间的相对顺序和取向从而实现局域延伸。基于该类技术方案，目前针对10X连接读长(linked read)测序技术，有如下组装工具：Architect(参考文献Kuleshov,V.,Snyder,M.,Batzoglou,S.Genome assembly from synthetic long read clouds,Bioinformatics,2016.)，ARKS(参考文献Coombe,L.,et al.ARKS:chromosome-scalescaffolding of human genome drafts with linked read kmers,BMC Bioinformatics,2018.)，Supernova(参考文献Weisenfeld，N.I.,et al.Direct determination ofdiploid genome sequences,Genome Research,2017.)。第二类是基于图算法的组装策略：该策略首先利用重叠群之间共享的条形码信息构建一个关系图，然后利用最小生成树算法确定关系图中重叠群之间相对顺序和取向从而实现全局的组装。基于该类技术方案，目前针对Illumina的读云测序技术，有以下组装工具：FragScaff(参考文献Adey,A.,et al.Invitro,long-range sequence information for de novo genome assembly viatransposase contiguity,Genome Research,2014.)。以上提到的组装工具中，ARKS仅仅利用了条形码信息进行骨架(scaffold)组装，而Architect、Supernova和FragScaff虽然能够同时利用条形码信息以及读长对(read pair)信息进行组装，但是它们都是基于读长对的插入片段长度满足正态分布的条件。由于stLFR测序技术与其他类似技术所得数据的特征完全不同，且目前没有专门针对stLFR测序数据特征开发的骨架组装工具，因此导致针对特定技术且基于延伸策略的组装软件无法对stLFR数据进行有效组装。而目前基于图理论的技术由于仅仅使用最小生成树(MST)算法并没有考虑全图的结构特征。该类策略过于贪婪导致其算法对输入数据的敏感度太高，因此该类技术的稳定性以及准确度也还有较大的提升空间。

概而言之，现有技术的缺点体现在：1)现有的技术方案，非常依赖于特定的实验数据特征，目前的技术方案仅仅针对读长对的插入片段长度满足正态分布的LFR技术，而stLFR数据中读长对的插入片段长度分布为非正态分布；2)现有的策略将骨架组装中重叠群之间的顺序和取向以及间距大小估计混在一起，模块化程度低，无法快速确定组装错误的来源，不利于组装策略的调整和扩展；3)现有的技术方案，无法屏蔽输入重叠群自身组装错误带来的影响，其会大大降低组装的效率和准确性；4)现有的自下而上的组装策略使骨架组装的速率较低。

发明内容

本申请提供一种基于单管长片段测序数据进行骨架组装的方法和装置，作为一个从头组装的技术方案，引入了stLFR数据所包含的长程信息，即利用来自同一条长片段的读长带有相同条形码信息以及来自同一个短片段的读长对信息来综合判断重叠群的排序和取向及它们之间的间距大小，以连接重叠群得到比较理想的骨架长度。

根据第一方面，一种实施例中提供一种基于单管长片段测序数据进行骨架组装的方法，包括：

获取stLFR测序数据，其包括条形码信息和读长对信息；

对上述stLFR测序数据进行预组装处理，获得种子重叠群；

利用上述条形码信息，确定上述种子重叠群之间的相对顺序；

利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重；

利用上述读长对信息，将未进行组装的重叠群进行挂载；

利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；

输出最终组装得到的骨架序列。

在优选实施例中，上述种子重叠群是其碱基平均覆盖度在设定范围内且序列长度在设定值以上的重叠群。

在优选实施例中，上述对上述stLFR测序数据进行预组装处理，获得种子重叠群，具体包括：

对stLFR对应的二代测序数据进行重叠群组装，并输出重叠群；

将stLFR测序数据的读长比对到输出的重叠群上；

解析上述重叠群所携带的条形码信息和读长对信息，获取可用的映射关系；

对重叠群进行识别和分类，得到种子重叠群和非种子重叠群。

在优选实施例中，上述利用上述条形码信息，确定上述种子重叠群之间的相对顺序，具体包括：

将上述种子重叠群按照设定的长度值进行分割，并计算出重叠群之间共享条形码的属性；

基于上述共享条形码的属性，构建上述种子重叠群之间的关系图；

基于最小生成树算法，获取上述种子重叠群之间的关系图对应的最小生成树图；

基于上述最小生成树图中分叉点的拓扑特征，屏蔽具有组装错误或不满足唯一性的种子重叠群；

对获取最小生成树图和屏蔽种子重叠群的步骤进行循环迭代，确定不包含可疑种子重叠群的关系图；

确定迭代最终得到的重叠群关系图对应的最小生成树图，并确定上述最小生成树图中种子重叠群之间的顺序。

在优选实施例中，上述利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重，具体包括：

获取已排序好的骨架中重叠群两端给定长度区间内的条形码集合；

计算重叠群两端序列区间与其所在骨架中左右两侧n阶近邻重叠群之间的相似度；

比较重叠群两端序列区间与两侧重叠群之间的相似性大小，确定重叠群在骨架中的不同取向对应的权重。

在优选实施例中，上述利用上述读长对信息，将未进行组装的重叠群进行挂载，具体包括：

将设定长度以上的重叠群按照给定长度值进行分割，并计算出重叠群之间共享条形码的属性；

根据上述共享条形码的属性，将骨架中近邻重叠群附近的非骨架中的重叠群进行聚类；

利用重叠群之间的读长对相连关系，对聚类的重叠群集合构建读长对连接关系图；

利用深度优先算法，探索上述读长对连接关系图中两个近邻重叠群之间的最短路径；

基于上述最短路径的信息，将相关的重叠群插入到骨架中，同时验证并修改两个种子重叠群之间的取向。

在优选实施例中，上述利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小，具体包括：

基于长度大于给定阈值的长重叠群，统计不同间距下两个窗口之间的Jaccard相似度；

利用最小二乘法计拟合得到相似度与间距的一一对应关系；

根据相似度与间距的对应关系以及两个重叠群之间的连接属性，估计骨架中相邻两个重叠群的间距。

根据第二方面，一种实施例中提供一种基于单管长片段测序数据进行骨架组装的装置，包括：

数据获取单元，用于获取stLFR测序数据，其包括条形码信息和读长对信息；

预组装处理单元，用于对上述stLFR测序数据进行预组装处理，获得种子重叠群；

顺序确定单元，用于利用上述条形码信息，确定上述种子重叠群之间的相对顺序；

权重确定单元，用于利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重；

重叠群挂载单元，用于利用上述读长对信息，将未进行组装的重叠群进行挂载；

间距估计单元，用于利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；

骨架输出单元，用于输出最终组装得到的骨架序列。

根据第三方面，一种实施例中提供一种基于单管长片段测序数据进行骨架组装的装置，包括：

存储器，用于存储程序；

处理器，用于通过执行上述存储器存储的程序以实现如第一方面的方法。

根据第四方面，一种实施例中提供一种计算机可读存储介质，其特征在于，包括程序，该程序能够被处理器执行以实现如第一方面的方法。

本发明是首个根据stLFR测序数据特征(读长对的插入片段长度呈非正态分布，条形码种类非常多)开发的骨架组装方法，其在使用stLFR数据中的读长对信息时，并不依赖特定读长对的插入片段长度分布。本发明的骨架组装方法可独立使用，前期的重叠群拼接使用的数据和拼接工具不受限制。本发明的骨架组装方法有效地提升了模块化程度以及可扩展性。

在优选实施例中，通过重叠群之间共享条形码的Jaccard相似度构建全局的关系图，然后综合使用最小生成树算法以及树图中节点的拓扑特征，能够确定并屏蔽一些特殊重叠群对于组装的干扰，从而大幅度提高该组装技术的稳定性以及组装结果中重叠群之间排序的准确性。

此外，在优选实施例中，采用根据重叠群之间共享条形码集合之间的相似度随距离单调递减的规律以及常规的读长对所包含的两条读长之间的方向实现对重叠群取向的双层确认，这样有助于提高确认重叠群在骨架中取向的准确性。

附图说明

图1为本发明实施例中基于人类基因组19号染色体的stLFR数据获取的双端测序数据的插入片段长度统计分布图；

图2为本发明实施例中基于单管长片段测序数据进行骨架组装的方法流程图；

图3为本发明实施例中对stLFR测序数据进行预组装处理获得种子重叠群的实现流程图；

图4为本发明实施例中确定种子重叠群之间的相对顺序的实现流程图；

图5为本发明实施例中确定种子重叠群之间相对取向的实现流程图；

图6为本发明实施例中挂载非种子重叠群到骨架中的实现流程图；

图7为本发明实施例中估计重叠群对之间的相对间隔的实现流程图；

图8为本发明实施例中基于单管长片段测序数据进行骨架组装的装置结构框图；

图9为本发明实施例中确定可屏蔽种子重叠群与不可屏蔽种子重叠群在树图中对应分叉点的拓扑特征图；

图10为本发明实施例中相似度随距离变化的统计与线性拟合关系图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。

另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。

通常基于二代数据的重叠群组装软件得到的组装结果其线性长度属性不够理想，无法实现基因组长程信号的分析。因此在获得重叠群之后，一般还需要利用具有大插入片段的双端测序信息或者其它具有长程信息的测序结果做进一步的骨架组装。现阶段利用长片段读长测序技术所携带长程信息做骨架组装的工具，或是需要利用满足正态分布的读长对信息来实现重叠群之间骨架组装，或是无法利用长片段读长测序数据中的读长对信息，因此目前还没有针对插入片段长度分布并不符合常规的正态分布的stLFR数据进行骨架组装的软件。

stLFR数据中读长对的插入片段长度分布如图1所示。本发明解决的技术问题包括以下几点：1)不局限于骨架组装对插入片段长度满足正态分布读长对条件的依赖。2)解决大部分的骨架组装技术中将确定重叠群的顺序、取向以及间距大小的过程混合在一起，没有将骨架过程模块化，不利于骨架组装技术扩展的问题。3)降低输入重叠群集合中具有组装错误的重叠群对于骨架组装效率和准确性的不利影响。4)利用自上而下的组装策略提升组装速率。

本发明的基于单管长片段测序数据进行骨架组装的方法，是一种基于图算法以及统计理论构建的从上至下的骨架组装方法，首先利用重叠群之间共享的条形码信息构建重叠群之间的关系图，然后利用图理论中的最小生成树算法以及于重叠群在树图中的拓扑特征迭代判断以确定重叠群之间的顺序，再利用相邻重叠群之间的首末端共享的条形码信息以及读长对信息确定重叠群之间的相对取向，最后利用长重叠群统计得到的条形码信息与间距的对应列表估计重叠群之间的相对间距，从而实现整个骨架组装过程，结果更为准确以及高效。

如图2所示，本发明的一种实施例中基于单管长片段测序数据进行骨架组装的方法，包括如下步骤：

S110：获取stLFR测序数据，其包括条形码信息和读长对信息；

S120：对上述stLFR测序数据进行预组装处理，获得种子重叠群；

S130：利用上述条形码信息，确定上述种子重叠群之间的相对顺序；

S140：利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重；

S150：利用上述读长对信息，将未进行组装的重叠群进行挂载；

S160：利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；

S170：输出最终组装得到的骨架序列。

如图3所示，本发明的一种实施例中，上述步骤S120：对stLFR测序数据进行预组装处理，获得种子重叠群，具体包括：

S121：对stLFR对应的二代测序数据进行重叠群组装，并输出重叠群。

在一个实施例中，利用成熟的二代数据组装软件如：SOAPdenovo(参考文献Li,R.,De novo assembly of human genomes with massively parallel short readsequencing,Genome Research,2009)，MaSuRCA(参考文献Zimin,A.V.,et al.The MaSuRCAgenome assembler,Bioinformatics,2013)等，将stLFR数据的二代读长信息进行重叠群组装，并生成初始的高质量重叠群集合。

S122：将stLFR测序数据的读长比对到输出的重叠群上。

在一个实施例中，使用成熟的比对工具BWA(参考文献Li H.and Durbin R.Fastand accurate short read alignment with Burrows-WheelerTransform.Bioinformatics(2009)，此工具可以从https://github.com/lh3/bwa免费下载)，将stLFR读长(read)比对到重叠群上。其它实施例中，比对软件不限于BWA，只要可生成或者比对结果可转化为SAM文件格式的比对工具均可。

S123：解析上述重叠群所携带的条形码信息和读长对信息，获取可用的映射关系。

在一个实施例中，基于比对输出的SAM文件，解析重叠群所携带的条形码信息以及读长对信息，具体包括：每个重叠群上所包含的条形码、读长对的种类、个数以及它们在重叠群上的位置分布情况。

S124：对重叠群进行识别和分类，得到种子重叠群和非种子重叠群。

本发明的一种实施例中，种子重叠群是其碱基平均覆盖度在设定范围内且序列长度在设定值以上的重叠群。例如，在一个实施例中，种子重叠群是指，其碱基平均覆盖度在0.5倍总体碱基平均覆盖度与1.5倍总体碱基平均覆盖度之间，且序列长度大于7000bp的重叠群。排除种子重叠群之后剩余的重叠群，为非种子重叠群。

在一个实施例中，通过设定重叠群上碱基平均覆盖度标准以及重叠群的长度标准，将重叠群分成两大类：第一类为种子重叠群，例如其碱基平均覆盖度在0.5倍总体碱基平均覆盖度与1.5倍总体碱基平均覆盖度之间，且序列长度大于7000bp；第二类为非种子重叠群，即排除种子重叠群之后剩余的重叠群。

进行骨架组装的目的是确定重叠群之间的相对顺序，相对取向以及估计相对间隔。考虑到条形码信息本身的对称性以及分辨率，本发明确定了自上而下的组装方案。如图4所示，本发明的一种实施例中，上述步骤S130：利用上述条形码信息，确定上述种子重叠群之间的相对顺序，具体包括：

S131：将上述种子重叠群按照设定的长度值进行分割，并计算出重叠群之间共享条形码的属性。

在一个实施例中，将种子重叠群按照给定长度(例如，7000bp)进行分割，每个等长的片段区间被称之为bin。结合步骤S123的信息，得到每个bin中的条形码集合情况。根据每个bin上的条形码分布集合，利用Jaccard系数计算两个bin之间的Jaccard相似系数，最后选取两个重叠群包含的bin之间的相似值的最高值作为两个重叠群之间的相似值。其中Jaccard相似系数定义为两个集合的交集与两个集合的并集比值：J(A,B)＝|A∩B|/|A∪B|。

S132：基于上述共享条形码的属性，构建上述种子重叠群之间的关系图。

在一个实施例中，根据种子重叠群之间的Jaccard相似系数(即共享条形码的属性)构建种子重叠群之间的关系图。关系图中顶点是种子重叠群，边表示两个种子重叠群之间存在满足超过阈值(例如0.1)的相似性，而相似系数则是对应边的权重。

S133：基于最小生成树算法，获取上述种子重叠群之间的关系图对应的最小生成树图。

在一个实施例中，在上一步构建的相似度关系图中，全部的种子重叠群并不会在同一个连通图中，所以常规的针对单个连通图的最小生成树算法无法直接使用。因此，首先需要利用不相交集合森林的算法(Robert E.Tarjan.Class notes:Disjoin setunion.COS 423,Princeton University,1999.)，将整个关系图划分为不同连通子图。然后使用最小生成树算法对每个连通子图的进行处理，从而得到对应的权重最大的树图，其中树图中边的权重为两个重叠群之间的Jaccard相似度。

S134：基于上述最小生成树图中分叉点的拓扑特征，屏蔽具有组装错误或不满足唯一性的种子重叠群。

在一个实施例中，基于重叠群之间相似性随重叠群之间在基因组上间距增大而减小的特征，最小生成树算法能够从图的全局信息出发获取每个重叠群的最近邻的重叠群。理想状态下每一个种子重叠群应该只有一个最近邻的种子重叠群，但实际上树图中存在大量的分叉点，从而导致组装效果变差。造成树图中存在大量分叉点(即连接度大于3的种子重叠群)，主要有两方面的原因：一个是由于条形码统计信息的涨落导致无法通过相似度区分最近邻和次近邻；另一个原因是该种子重叠群并不满足唯一性或者具有组装错误。基于统计理论分析，两种不同原因导致的分叉点在树图中具有不同的统计特征，如图9所示。第一种原因导致的分叉点，其连接的某一个分叉路径非常短，而其它分叉路径很长，如图9A所示；而第二种原因导致的分叉点，其连接的所有分叉路径都很长，如图9B所示。将图9A以及图9B中的重叠群比对到参考基因组，能得到它们之间的正确排序，分别如图9C和9D所示。同时也发现图9B中的分叉点无法比对回参考基因组，这也就验证了该种子重叠群是具有组装错误的。基于以上的特征和原理，在具体实施例中，首先通过设定长度阈值(例如3)删除树图中较短的枝杈，然后在剩余的树图中识别那些具有多个长枝杈的分叉点，最后标记这些分叉点为可疑种子重叠群，并将其从关系图中屏蔽。

S135：对获取最小生成树图和屏蔽种子重叠群的步骤进行循环迭代，确定不包含可疑种子重叠群的关系图。

在一个实施例中，基于以上屏蔽的种子重叠群关系图，重复S134以及S135步骤，直到所有可疑种子重叠群或者屏蔽的种子重叠群数目超过一定阈值时，重复迭代停止，并且输出最后的树图。

S136：确定迭代最终得到的重叠群关系图对应的最小生成树图，并确定上述最小生成树图中种子重叠群之间的顺序。

在一个实施例中，基于以上输出的树图，将每一个长度超过阈值(例如，枝杈上种子重叠群数目超过4)的枝杈作为已排好序的骨架进行输出，其中骨架的中种子重叠群的顺序定义为它们在枝杈上的拓扑顺序。

如图5所示，本发明的一种实施例中，上述步骤S140：利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重，具体包括：

S141：获取已排序好的骨架中重叠群两端给定长度区间内的条形码集合。

在一个实施例中，截取每个已排序的骨架中包含的种子重叠群两端给定长度(例如，3500bp)区间，统计给定区间上条形码的种类和数目。

S142：计算重叠群两端序列区间与其所在骨架中左右两侧n阶近邻重叠群之间的相似度。

在一个实施例中，针对每一个重叠群，计算出其两端序列与其所在骨架中左右两侧给定阶数(例如，4阶)的近邻种子重叠群之间的相似度，生成四个相似度列表，其包括左侧与左侧重叠群之间的相似度列表，右侧与左侧重叠群之间的相似度列表，左侧与右侧重叠群之间的相似度列表以及右侧与右侧的相似度列表。

S143：比较重叠群两端序列区间与两侧重叠群之间的相似性大小，确定重叠群在骨架中的不同取向对应的权重。

在一个实施例中，根据stLFR数据的统计特征，在理想的骨架中重叠群左侧区间与左侧重叠群的相似度更高，而右侧区间与右侧重叠群之间的相似度更高。基于以上原则，通过比较重叠群左侧与左侧重叠群之间的相似度列表和右侧与左侧重叠群之间的相似度列表，以及左侧与右侧重叠群之间的相似度列表和右侧与右侧的相似度列表，并统计其大小的个数来确定每个重叠群在骨架中的取向权重。

由于以上步骤仅仅对未屏蔽的种子重叠群进行顺序和取向的确认，还存在大量短小且在局域满足唯一性的重叠群并没有参与骨架组装，因此，将通过结合条形码局域筛选以及读长对连通信息实现将未参与骨架组装的重叠群挂载到骨架中。根据步骤S121中得到的读长对在重叠群上分布信息，利用跨重叠群的读长对将部分没有进入到骨架中的重叠群插入；同时利用填充的读长对的信息修正骨架中相邻重叠群之间的相对取向。

如图6所示，步骤S150：利用上述读长对信息，将未进行组装的重叠群进行挂载，具体包括：

S151：将设定长度以上的重叠群按照给定长度值进行分割，并计算出重叠群之间共享条形码的属性。

在一个实施例中，将大于给定长度(例如，1000bp)且未参与骨架的重叠群进行筛选，对筛选出来的重叠群以及骨架中的重叠群按照给定长度(例如，1000bp)进行分割，每个等长的片段区间也被称之为bin。结合步骤S123的信息，得到每个bin中的条形码集合情况。根据每个bin上的条形码分布集合，计算两个bin之间的共享条形码数目，最后选取两个重叠群包含的bin之间的共享条形码的最大值作为两个重叠群之间的共享条形码数目。

S152：根据上述共享条形码的属性，将骨架中近邻重叠群附近的非骨架中的重叠群进行聚类。

在一个实施例中，对于骨架中每一对近邻的重叠群对，基于它们与未参与骨架组装的重叠群之间共享条形码总数，将该近邻重叠群对附近的未参与骨架中的重叠群进行聚类。该聚类集合定义为与重叠群对中任意重叠群具有满足大于给定共享条形码阈值(例如，10)的所有尚未参与骨架组装的重叠群。

S153：利用重叠群之间的读长对相连关系，对聚类的重叠群集合构建读长对连接关系图。

在一个实施例中，读长对连接关系图中顶点是聚类的重叠群集合与骨架中的近邻重叠群对中的任意重叠群，边表示两个重叠群之间存在读长对关系。由于读长对之间存在确定的相对取向，因此该关系图是有向图。

S154：利用深度优先算法，探索上述读长对连接关系图中两个近邻重叠群之间的最短路径。

在一个实施例中，基于重叠群之间通过读长对构建的连通有向图，利用深度优先算法，探索并确认该有向图中两个近邻重叠群之间的最短连通路径。

S155：基于上述最短路径的信息，将相关的重叠群插入到骨架中，同时验证并修改两个种子重叠群之间的取向。

在一个实施例中，实现逻辑如下：如果不存在最短连通路径，那么骨架近邻重叠群对之间不插入任何新的重叠群；如果存在最短连通路径，那么将最短连通路径中重叠群之间排序和取向来替换之前的重叠群对。因此如果最短路径中起始重叠群对之间的取向与骨架中该重叠群对之间的取向一致，即验证了之前结果；如果两者的取向不一致，则是利用读长对信息所提供的取向信息修正之前的取向，最终达到利用读长对信息验证并修正条形码信息给出的两个种子重叠群之间相对取向的目的。

如图7所示，本发明的一种实施例中，上述步骤S160：利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小，具体包括：

S161：基于长度大于给定阈值的长重叠群，统计不同间距下两个窗口之间的Jaccard相似度。

在一个实施例中，利用长度大于给定阈值(例如，3500)的条件，挑选出足够数目的重叠群。对每一个重叠群以给定大小(例如，3500)的窗口进行给定等间距(例如，3500)的滑动，统计不同间距下两个窗口所包含的条形码信息，然后计算并统计出窗口在不同间距的条件下的Jaccard相似系数。

S162：利用最小二乘法计拟合得到相似度与间距的一一对应关系。

在一个实施例中，根据相似度的统计分布的线性属性，利用最小二乘法拟合得到的具体的线性函数形式，其统计分布与线性拟合效果如图10所示。然后根据该函数，在给定相似度值间隔条件下，给出窗口相似度与窗口之间间距的一一对应列表。在具体的实施过程中，删除统计数据中的奇异点，即其与同一相似度下对应的间距值偏离过于明显。

S163：根据相似度与间距的对应关系以及两个重叠群之间的连接属性，估计骨架中相邻两个重叠群的间距。

在一个实施例中，当骨架中两个近邻重叠群之间存在读长对连接信息时，将在重叠群对之间加入固定长度的N(例如，5)；当两个近邻重叠群之间不存在读长对连接信息时，将利用两个重叠群之间的相似度系数对应列表中找到该相似度系数对应的间距值，然后在两个重叠群之间加入该间距值对应个数的N。

对应于本发明实施例的基于单管长片段测序数据进行骨架组装的方法，本发明还提供一种基于单管长片段测序数据进行骨架组装的装置，如图8所示，包括：数据获取单元810，用于获取stLFR测序数据，其包括条形码信息和读长对信息；预组装处理单元820，用于对上述stLFR测序数据进行预组装处理，获得种子重叠群；顺序确定单元830，用于利用上述条形码信息，确定上述种子重叠群之间的相对顺序；权重确定单元840，用于利用上述条形码信息，确定近邻种子重叠群之间不同取向的权重；重叠群挂载单元850，用于利用上述读长对信息，将未进行组装的重叠群进行挂载；间距估计单元860，用于利用上述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；骨架输出单元870，用于输出最终组装得到的骨架序列。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

因此，本发明还提供一种基于单管长片段测序数据进行骨架组装的装置，包括：存储器，用于存储程序；处理器，用于通过执行上述存储器存储的程序以实现如上所述的基于单管长片段测序数据进行骨架组装的方法。

本发明还提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如上所述的基于单管长片段测序数据进行骨架组装的方法。

以下通过实施例详细说明本发明的技术方案，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。

实施例1

本实施例提供一种基于单管长片段测序数据进行骨架组装的方法，其具体步骤如下：

步骤一、stLFR骨架组装数据的预处理和准备

在该实施例中，将通过如下步骤对用于骨架组装的stLFR数据进行预处理和准备：

1.1利用成熟的二代数据组装软件如：SOAPdenovo(参考文献Li,R.,De novoassembly of human genomes with massively parallel short read sequencing,Genome Research,2009)，MaSuRCA(参考文献Zimin,A.V.,et al.The MaSuRCA genomeassembler,Bioinformatics,2013)等，将stLFR数据的二代读长信息进行重叠群组装，并生成初始的高质量重叠群集合。

1.2使用成熟的比对工具BWA(参考文献Li H.and Durbin R.Fast and accurateshort read alignment with Burrows-Wheeler Transform.Bioinformatics(2009)，此工具可以从https://github.com/lh3/bwa免费下载)，将stLFR读长(read)比对到重叠群上。其它实施例中，比对软件不限于BWA，只要可生成或者比对结果可转化为SAM文件格式的比对工具均可。

1.3基于比对输出的SAM文件，解析重叠群所携带的条形码信息以及读长对信息，具体包括：每个重叠群上所包含的条形码、读长对的种类、个数以及它们在重叠群上的位置分布情况。

1.4通过设定重叠群上碱基平均覆盖度标准以及重叠群的长度标准，将重叠群分成两大类：第一类为种子重叠群，例如其碱基平均覆盖度在0.5倍总体碱基平均覆盖度与1.5倍总体碱基平均覆盖度之间，且序列长度大于7000bp；第二类为非种子重叠群，即排除种子重叠群之后剩余的重叠群。

步骤二、利用stLFR数据的条形码信息，确定种子重叠群之间的相对顺序

进行骨架组装的目的是确定重叠群之间的相对顺序，相对取向以及估计相对间隔。考虑到条形码信息本身的对称性以及分辨率，本实施例确定了自上而下的组装方案，该步骤的目的就是确定种子重叠群之间的相对顺序。具体包括如下步骤：

2.1将种子重叠群按照给定长度(例如，7000bp)进行分割，每个等长的片段区间被称之为bin。结合步骤1.3的信息，得到每个bin中的条形码集合情况。根据每个bin上的条形码分布集合，利用Jaccard系数计算两个bin之间的Jaccard相似系数，最后选取两个重叠群包含的bin之间的相似值的最高值作为两个重叠群之间的相似值。其中Jaccard相似系数定义为两个集合的交集与两个集合的并集比值：J(A,B)＝|A∩B|/|A∪B|。

2.2根据种子重叠群之间的Jaccard相似系数构建种子重叠群之间的关系图。关系图中顶点是种子重叠群，边表示两个种子重叠群之间存在满足超过阈值(例如0.1)的相似性，而相似系数则是对应边的权重。

2.3在上一步构建的相似度关系图中，全部的种子重叠群并不会在同一个连通图中，所以常规的针对单个连通图的最小生成树算法无法直接使用。因此，首先需要利用不相交集合森林的算法(Robert E.Tarjan.Class notes:Disjoin set union.COS 423,Princeton University,1999.)，将整个关系图划分为不同连通子图。然后使用最小生成树算法对每个连通子图的进行处理，从而得到对应的权重最大的树图，其中树图中边的权重为两个重叠群之间的Jaccard相似度。

2.4基于树图的拓扑特征，屏蔽关系图中具有组装错误或者不满足唯一性条件的种子重叠群。基于重叠群之间相似性随重叠群之间在基因组上间距增大而减小的特征，最小生成树算法能够从图的全局信息出发获取每个重叠群的最近邻的重叠群。理想状态下每一个种子重叠群应该只有一个最近邻的种子重叠群，但实际上树图中存在大量的分叉点，从而导致组装效果变差。造成树图中存在大量分叉点(即连接度大于3的种子重叠群)，主要有两方面的原因：一个是由于条形码统计信息的涨落导致无法通过相似度区分最近邻和次近邻；另一个原因是该种子重叠群并不满足唯一性或者具有组装错误。基于统计理论分析，两种不同原因导致的分叉点在树图中具有不同的拓扑特征，如图9所示。第一种原因导致的分叉点，其连接的某一个分叉路径非常短，而其它分叉路径很长，如图9A所示；而第二种原因导致的分叉点，其连接的所有分叉路径都很长，如图9B所示。将图9A以及图9B中的重叠群比对到参考基因组，能得到它们之间的正确排序，分别如图9C和9D所示。同时也发现图9B中的分叉点无法比对回参考基因组，这也就验证了该重叠群是具有组装错误的。基于以上的特征和原理，在具体实施例中，首先通过设定长度阈值(例如3)删除树图中较短的枝杈，然后在剩余的树图中识别那些具有多个长枝杈的分叉点，最后标记这些分叉点为可疑种子重叠群，并将其从关系图中屏蔽。

2.5基于以上屏蔽的种子重叠群关系图，重复2.3以及2.4步骤，直到所有可疑种子重叠群或者屏蔽的种子重叠群数目超过一定阈值时，重复迭代停止，并且输出最后的树图。

2.6基于以上输出的树图，将每一个长度超过阈值(例如，枝杈上种子重叠群数目超过4)的枝杈作为已排好序的骨架进行输出，其中骨架中种子重叠群的顺序定义为它们在枝杈上的拓扑顺序。

步骤三、利用条形码信息，确定近邻种子重叠群之间不同取向的权重

基于给定顺序的重叠群集合，对重叠群进行取向的确定将有利于局域信息的进一步使用。

具体包括如下步骤：

3.1截取每个已排序的骨架中包含的种子重叠群两端给定长度(例如，3500bp)区间，统计给定区间上条形码的种类和数目。

3.2针对每一个重叠群，计算出其两端序列与其所在骨架中左右两侧给定阶数(例如，4阶)的近邻种子重叠群之间的相似度，生成四个相似度列表，其包括左侧与左侧重叠群之间的相似度列表，右侧与左侧重叠群之间的相似度列表，左侧与右侧重叠群之间的相似度列表以及右侧与右侧的相似度列表。

3.3根据stLFR数据的统计特征，在理想的骨架中重叠群左侧区间与左侧重叠群的相似度更高，而右侧区间与右侧重叠群之间的相似度更高。基于以上原则，通过比较重叠群左侧与左侧重叠群之间的相似度列表和右侧与左侧重叠群之间的相似度列表，以及左侧与右侧重叠群之间的相似度列表和右侧与右侧的相似度列表，并统计其大小的个数来确定每个重叠群在骨架中的取向权重。

步骤四、利用读长对信息，将未进行骨架组装的重叠群进行挂载

由于以上步骤仅仅对未屏蔽的种子重叠群进行顺序和取向的确认，还存在大量短小且在局域满足唯一性的重叠群并没有参与骨架组装，因此在该实施例中，将通过结合条形码局域筛选以及读长对连通信息实现将未参与骨架组装的重叠群挂载到骨架中。根据步骤1.2中得到的读长对在重叠群上分布信息，利用跨重叠群的读长对将部分没有进入到骨架中的重叠群插入；同时利用填充的读长对的信息修正骨架中相邻重叠群之间的相对取向。具体包括如下步骤：

4.1将大于给定长度(例如，1000bp)且未参与骨架的重叠群进行筛选，对筛选出来的重叠群以及骨架中的重叠群按照给定长度(例如，1000bp)进行分割，每个等长的片段区间也被称之为bin。结合步骤1.3的信息，得到每个bin中的条形码集合情况。根据每个bin上的条形码分布集合，计算两个bin之间的共享条形码数，最后选取两个重叠群包含的bin之间的共享条形码的最大值作为两个重叠群之间的共享条形码数目。

4.2对于骨架中每一对近邻的重叠群对，基于它们与未参与骨架组装的重叠群之间共享条形码总数，将该近邻重叠群对附近的未参与骨架中的重叠群进行聚类。该聚类集合定义为与重叠群对中任意重叠群具有满足大于给定共享条形码阈值(例如，10)的所有尚未参与骨架的重叠群。

4.3利用重叠群之间的读长对相连关系，对聚类的重叠群集合构建读长对关系图。该关系图中顶点是聚类的重叠群集合与骨架中的近邻重叠群对中的任意重叠群，边表示两个重叠群之间存在读长对关系。由于读长对之间存在确定的相对取向，因此该关系图是有向图。

4.4基于重叠群之间通过读长对构建的连通有向图，利用深度优先算法，探索并确认该有向图中两个近邻重叠群之间的最短连通路径。

4.5基于读长对有向图中的最短连通路径的信息，将相关的重叠群挂载到骨架中。具体实现逻辑如下：如果不存在最短连通路径，那么骨架近邻重叠群对之间不插入任何新的重叠群；如果存在最短连通路径，那么将最短连通路径中重叠群之间排序和取向来替换之前的重叠群对。因此如果最短路径中起始重叠群对之间的取向与骨架中该重叠群对之间的取向一致，即验证了之前结果；如果两者的取向不一致，则是利用读长对信息所提供的取向信息修正之前的取向，最终达到利用读长对信息验证并修正条形码信息给出的两个种子重叠群之间相对取向的目的。

步骤五、利用条形码以及读长对信息特征估计重叠群之间的相对间距大小

5.1利用长度大于给定阈值(例如，3500)的条件，挑选出足够数目的重叠群。对每一个重叠群以给定大小(例如，3500)的窗口进行给定等间距(例如，3500)的滑动，统计不同间距下两个窗口所包含的条形码信息，然后计算并统计出窗口在不同间距的条件下的Jaccard相似系数。

5.2根据相似度的统计分布的线性属性，利用最小二乘法计拟合得到的具体的线性函数形式，其统计分布与线性拟合效果如图10所示。然后根据该函数，在给定相似度值间隔条件下，给出窗口相似度与窗口之间间距的一一对应列表。在具体的实施过程中，删除统计数据中的奇异点，即其与同一相似度下对应的间距值偏离过于明显。

5.3根据Jaccard相似系数与窗口间距的对应列表，以及两个重叠群之间的连接属性来估计骨架中相邻两个重叠群的间距。在具体的实施过程中，当骨架中两个近邻重叠群之间存在读长对连接信息时，将在重叠群对之间加入固定长度的N(例如，5)；当两个近邻重叠群之间不存在读长对连接信息时，将利用两个重叠群之间的相似度系数对应列表中找到该相似度系数对应的间距值，然后在两个重叠群之间加入该间距值对应个数的N。

步骤六、输出最终的骨架

该步骤将汇总步骤二、三、四、五的信息，确定重叠群之间最终的顺序、取向和间距大小，然后以标准的fasta文件格式输出最终的骨架。

实施例2

为了确定本发明的基于单管长片段测序数据进行骨架组装的方法，相对于常规无法利用条形码信息进行骨架组装的软件的性能提升，本实施例测试并对比了基于不同二代组装软件与本发明的方法在对人类基因组中19号染色体的stLFR数据组装结果，具体实施过程如下：

1)从中国国家基因库核酸序列归档系统数据库(https://db.cngb.org/cnsa/)下载得到人类基因组对应的stLFR原始测序数据(其对应的归档编号为CNP0000066)，然后基于GRCh38参考基因组，提取能够比对到19号染色体的所有stLFR数据。该数据集包含了78X的双端测序数据，fastq格式，共111G，每一个读长的长度是100bp。

2)利用SOAPdenovo、MaSuRCA对19号染色体的stLFR数据进行组装，得到最终的骨架；上述过程均使用8个线程。

3)利用SOAPdenovo以及MaSuRCA的组装结果准备本发明方法的输入重叠群。由于SOAPdenovo具有独立的重叠群组装模块并输出相应的重叠群，因此能够直接使用其重叠群作为输入，而MaSuRCA并没有独立的重叠群模块，因此需要将其输出的骨架从N处打断，获取相应的重叠群；上述过程均使用8个线程。

4)利用本发明方法对不同的重叠群进行骨架组装，本方案的实验编号为stLFR骨架组装(assembler)。

5)利用评估软件QUAST(参考文献Gurevich,A.,et al.QUAST:qualityassessment tool for genome assemblies,Bioinformatics,2013)，在默认参数的条件下，对不同的组装结果进行评估，具体评估结果见下表1。其中骨架中的重叠群是指由骨架在连续包含十个及以上未知碱基的空位处打断而形成的重叠群。骨架中的无错重叠群是指由骨架在连续包含十个及以上未知碱基的空位处或者组装错误处打断而形成的重叠群。无错骨架是指由初始骨架在组装错误处打断而形成的新骨架。骨架长度/重叠群的NG50表示由一组由骨架/重叠群长度构成的数字集合对应的统计量。将集合中大于NG50的数求和，其和小于给定基因组总长一半；而将大于等于NG50的数求和，其和大于等于给定基因组总长的一半。以上评估的具体实现算法可以参考该软件对应的论文(参考文献Gurevich,A.,etal.QUAST:quality assessment tool for genome assemblies,Bioinformatics,2013)。

表1.stLFR骨架组装对于不同结果的组装提升效果以及资源消耗情况

6)结论：针对19号染色体的数据，相对于SOAPdenovo软件产生的骨架，stLFR骨架组装(assembler)在SOAPdenovo输出的重叠群基础上额外利用了stLFR数据的条形码信息，使骨架中的重叠群的N50提升了2.4倍，其相应的NGA50也提升了5.3倍；同时骨架的N50提升了546倍。相对于MaSuRCA，利用条形码信息进行骨架组装后，使骨架的N50提升了249倍，同时对于相应的NGA50也提升了33倍。这些结果表明本发明的方法能够利用条形码带来的长程信息实现基因组从头组装性能的提升，包括组装的准确性以及组装结果的连续性，从而为下游生物信息分析提供基础保障。

从时间和内存资源消耗方面，从表1可以看到，stLFR骨架组装的时间消耗仅仅是SOAPdenovo完成整个组装消耗时间的130％，是MaSuRCA完成组装所消耗时间的88％，而内存消耗峰值为SOAPdenovo整个组装过程中峰值的8％，是MaSuRCA内存峰值的20％。这些结果表明本发明方法的计算性能非常好，能够适用于各种大小的基因组组装。

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (10)

1.一种基于单管长片段测序数据进行骨架组装的方法，其特征在于，所述方法包括：

获取stLFR测序数据，其包括条形码信息和读长对信息；

对所述stLFR测序数据进行预组装处理，获得种子重叠群；

利用所述条形码信息，确定所述种子重叠群之间的相对顺序；

利用所述条形码信息，确定近邻种子重叠群之间不同取向的权重；

利用所述读长对信息，将未进行组装的重叠群进行挂载；

利用所述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；

输出最终组装得到的骨架序列。

2.根据权利要求1所述的方法，其特征在于，所述种子重叠群是其碱基平均覆盖度在设定范围内且序列长度在设定值以上的重叠群。

3.根据权利要求1所述的方法，其特征在于，所述对所述stLFR测序数据进行预组装处理，获得种子重叠群，具体包括：

对stLFR对应的二代测序数据进行重叠群组装，并输出重叠群；

将stLFR测序数据的读长比对到输出的重叠群上；

解析所述重叠群所携带的条形码信息和读长对信息，获取可用的映射关系；

对重叠群进行识别和分类，得到种子重叠群和非种子重叠群。

4.根据权利要求1所述的方法，其特征在于，所述利用所述条形码信息，确定所述种子重叠群之间的相对顺序，具体包括：

将所述种子重叠群按照设定的长度值进行分割，并计算出重叠群之间共享条形码的属性；

基于所述共享条形码的属性，构建所述种子重叠群之间的关系图；

基于最小生成树算法，获取所述种子重叠群之间的关系图对应的最小生成树图；

基于所述最小生成树图中分叉点的拓扑特征，屏蔽具有组装错误或不满足唯一性的种子重叠群；

对获取最小生成树图和屏蔽种子重叠群的步骤进行循环迭代，确定不包含可疑种子重叠群的关系图；

确定迭代最终得到的重叠群关系图对应的最小生成树图，并确定所述最小生成树图中种子重叠群之间的顺序。

5.根据权利要求1所述的方法，其特征在于，所述利用所述条形码信息，确定近邻种子重叠群之间不同取向的权重，具体包括：

获取已排序好的骨架中重叠群两端给定长度区间内的条形码集合；

计算重叠群两端序列区间与其所在骨架中左右两侧n阶近邻重叠群之间的相似度；

比较重叠群两端序列区间与两侧重叠群之间的相似性大小，确定重叠群在骨架中的不同取向对应的权重。

6.根据权利要求1所述的方法，其特征在于，所述利用所述读长对信息，将未进行组装的重叠群进行挂载，具体包括：

将设定长度以上的重叠群按照给定长度值进行分割，并计算出重叠群之间共享条形码的属性；

根据所述共享条形码的属性，将骨架中近邻重叠群附近的非骨架中的重叠群进行聚类；

利用重叠群之间的读长对相连关系，对聚类的重叠群集合构建读长对连接关系图；

利用深度优先算法，探索所述读长对连接关系图中两个近邻重叠群之间的最短路径；

基于所述最短路径的信息，将相关的重叠群插入到骨架中，同时验证并修改两个种子重叠群之间的取向。

7.根据权利要求1所述的方法，其特征在于，所述利用所述条形码信息和读长对信息特征估计重叠群之间的相对间距大小，具体包括：

基于长度大于给定阈值的长重叠群，统计不同间距下两个窗口之间的Jaccard相似度；

利用最小二乘法计拟合得到相似度与间距的一一对应关系；

根据相似度与间距的对应关系以及两个重叠群之间的连接属性，估计骨架中相邻两个重叠群的间距。

8.一种基于单管长片段测序数据进行骨架组装的装置，其特征在于，所述装置包括：

数据获取单元，用于获取stLFR测序数据，其包括条形码信息和读长对信息；

预组装处理单元，用于对所述stLFR测序数据进行预组装处理，获得种子重叠群；

顺序确定单元，用于利用所述条形码信息，确定所述种子重叠群之间的相对顺序；

权重确定单元，用于利用所述条形码信息，确定近邻种子重叠群之间不同取向的权重；

重叠群挂载单元，用于利用所述读长对信息，将未进行组装的重叠群进行挂载；

间距估计单元，用于利用所述条形码信息和读长对信息特征估计重叠群之间的相对间距大小；

骨架输出单元，用于输出最终组装得到的骨架序列。

9.一种基于单管长片段测序数据进行骨架组装的装置，其特征在于，所述装置包括：

存储器，用于存储程序；

处理器，用于通过执行所述存储器存储的程序以实现如权利要求1-7中任一项所述的方法。

10.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1-7中任一项所述的方法。