CN115910199A - 一种基于比对框架的三代测序数据结构变异检测方法 - Google Patents

Abstract

一种基于比对框架的三代测序数据结构变异检测方法，涉及基因生物学技术领域。本发明是为了解决现有针对三代测序数据的变异检测方法还存在检测时间长、计算量大的问题。本发明包括：利用测序序列获取精确匹配片段：采用模糊比对方法对精确匹配片段进行延伸，获得模糊匹配片段；采用稀疏动态规划算法利用模糊匹配片段构建比对框架；对比对框架上的模糊匹配片段进行扩展获取最大匹配块，对最大匹配块的一致性关系解析，获得结构变异；对结构变异进行聚类，将同一类型的结构变异进行合并获得每个类型的结构变异簇；根据每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列。本发明用于检测结构变异。

Description

一种基于比对框架的三代测序数据结构变异检测方法

技术领域

本发明涉及基因生物学技术领域，特别涉及一种基于比对框架的三代测序数据结构变异检测方法。

背景技术

现如今新冠病毒横行，未来也有可能出现新的未知疾病，科学家需要在极短时间内研究出应对之法。所以说如何在短时间内通过人类样本数据检测出未知变异变得尤为重要。

目前，针对三代测序数据的变异检测技术有很多，其中主要以Sniffles、cuteSV、SVIM和PBSM为代表。这类检测工具都是以序列比对算法为前提的变异检测工具，这就导致了整个变异检测需要花费相当一部分时间用于前期的数据处理。同时碱基序列比对结果的后处理、结构变异的识别也都是需要大量计算，这就会导致整个流程时间消耗大大增加。

发明内容

本发明目的是为了解决现有针对三代测序数据的变异检测方法还存在检测时间长、计算量大的问题，而提出了一种基于比对框架的三代测序数据结构变异检测方法。

一种基于比对框架的三代测序数据结构变异检测方法具体过程为：

步骤一、获取测序序列，利用测序序列获取精确匹配片段：

步骤一一、获取测序序列，按照预设步长K划分每条测序序列，利用每个长度为k的k-mer片段与预先建立的GRM图索引结构进行对比，获得每一个碱基都能够完全比对成功的k-mer片段作为候选种子片段；

所述GRM图索引结构通过测序序列建立；

步骤一二、对候选种子片段的上下游进行扩展直到上下游都遇到无法匹配的碱基为止，获得精确匹配片段；

步骤二、采用Landau-Vishkin模糊比对方法对精确匹配片段进行延伸，直至匹配错误的碱基数量大于预设错误匹配阈值时结束延伸，获得模糊匹配片段；

采用Landau-Vishkin模糊比对方法对每一个精确匹配片段进行延伸是在GRM图索引结构中进行的；

所述精确匹配片段通过五元组形式记录，具体为：染色体编号，在参考基因组上的起始位置，在参考基因组上的终止位置，在测序片段上的起始位置，在测序片段上的终止位置；

步骤三、采用稀疏动态规划算法利用步骤二获得的模糊匹配片段构建比对框架；

步骤四、对步骤三获得的比对框架上的模糊匹配片段进行扩展获取最大匹配块，对最大匹配块的一致性关系解析，获得结构变异；

所述一致性关系包括：一致性关系、非一致性关系；

步骤五、对步骤四获得的结构变异进行聚类，将同一类型的结构变异进行合并获得每个类型的结构变异簇；

步骤六、根据步骤五获取每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列。

进一步地，所述步骤三中的采用稀疏动态规划算法利用步骤二获得的模糊匹配片段构建比对框架，包括以下步骤：

步骤三一、将步骤二获得的模糊匹配片段进行排序：

首先，按照染色体编号将模糊匹配片段从大到小进行排序；

然后，将染色体编号相同的模糊匹配片段按照基因组起始位置的顺序进行排序；

步骤三二、利用排序后的模糊匹配片段构建无回路的有向图；

步骤三三、对步骤三二获得的无回路的有向图上的每个顶点进行打分；

步骤三四、利用启发式动态规划算法按照分数从高到低依次选择顶点组成最优路径，一条或多条最优路径组成比对框架。

进一步地，所述步骤三二中的利用排序后的模糊匹配片段构建无回路的有向图，具体为：

将模糊匹配片段作为顶点，如果两个顶点位于同一条染色体上，且在参考基因组上的位置偏差和在测序片段上的位置偏差均小于预设最大跨度阈值，则两个顶点之间的连线构成一条有向边；

所述有向边的方向从在参考基因组上的起始位置值小的顶点指向起始位置值大的顶点。

进一步地，所述步骤三三中的对步骤三二获得的无回路的有向图上的每个顶点进行打分，如下式：

S(Aj)＝max{S(Ai)+W-P}

其中，Aj和Ai表示连接后的边上的两个锚点j和i，W是锚点Aj的权重值，P表示Ai至Aj的罚分。

进一步地，所述步骤四中的对步骤三获得的比对框架上的模糊匹配片段进行扩展获取最大匹配块，对最大匹配块的一致性关系解析，获得结构变异，包括以下步骤：

步骤四一、对步骤三构建的比对框架中顶点的上下游进行扩展，获取最大匹配块，并获取最大匹配块之间的一致性关系，具体为：

步骤四一一、在比对框架中按照顶点长度从大到小依次选取顶点Ci作为延伸目标；

步骤四一二、预设窗口w，获取预设窗口w内顶点Ci的扩展比对片段A在参考基因组和测序序列上的坐标，并根据获取的坐标填补顶点Ci与顶点Cj之间的孔隙，获得最大匹配块及最大匹配块之间的一致性关系；

其中，Cj是Ci下游相邻的顶点；

所述最大匹配块之间的一致性关系包括：一致性关系、非一致性关系；

步骤四二、从满足不一致性关系的最大匹配块中解析出结构变异。

进一步地，所述步骤四一二中的预设窗口w，获取预设窗口w内顶点Ci的扩展比对片段A在参考基因组和测序序列上的坐标，然后利用A的坐标填补顶点Ci与Cj之间的孔隙，获得最大匹配块及最大匹配块之间的一致性关系，包括以下步骤：

S1、将顶点Ci和顶点Cj之间的间隙按照预设窗口长度w分割为多个窗口区域；

S2、利用哈希结构将窗口w内的测序序列片段映射到参考基因组上，获得能够与参考基因组比对成功的测序序列片段即顶点Ci在窗口w内的扩展比对片段A，获得A在参考基因组上的坐标；

S3、获取A在测序序列上的坐标，根据A在参考基因组上的坐标和A在测序序列片段上的坐标获得A在参考基因组和测序序列片段的位置偏差，如果A在参考基因组和测序序列片段的位置偏差小于预设距离阈值，则将A填补到Ci与Cj之间的孔隙中，获得最大匹配块：

S031、首先，获取A在参考基因组和测序序列片段的位置偏差：

将A在参考基因组上位置坐标与Ci在参考基因组上的位置坐标相减获得位置偏差a；

将A在测序序列上位置坐标与Ci在测序序列上的位置坐标相减获得位置偏差b；

获取b-a的值，若b-a<距离阈值15bp，则将扩展比对片段A填补到Ci与Cj之间的孔隙中，执行S302；

S302，将扩展比对片段A加入顶点Ci并更新Ci在参考基因组上的坐标信息和Ci在测序序列上的坐标信息，如果是上游扩展，取A在参考基因组上的起始位置作为Ci在参考基因组上的起始位置，A在测序序列上的起始位置作为Ci在测序序列上的终止位置，如果是下游扩展，则取A在参考基因组上的终止坐标作为Ci在参考基因组上的起始位置，A在测序序列上的终止坐标作为Ci在测序序列上的终止坐标；

所述上游扩展为：A在参考基因组上的坐标小于更新后的Ci在参考基因组上的坐标；

A在测序序列上的坐标小于更新后的Ci在测序序列上的坐标；

所述下游扩展为：A在参考基因组上的坐标大于更新后的Ci在参考基因组上的坐标；A在测序序列上的坐标大于更新后Ci在测序序列上的坐标；

S4，按照长度w选取下一个窗口，重复执行S2-S3直至选取的窗口与顶点Cj相邻，获得当前加入A后的顶点即为最大匹配块；若不同顶点形成的最大匹配块之间有交叠，则最大匹配块满足一致性关系；若不同顶点形成的最大匹配块之间没有重叠，则最大匹配块满足非一致性关系。

进一步地，所述步骤四二中的从满足不一致性关系的最大匹配块中解析出结构变异，具体如下：

插入结构变异：变异区间内有一个非一致性单元，同时非一致性单元在测序序列上的位置相较于非一致性单元在参考基因组上的位置偏差长度大于+50bp；

所述非一致性单元为满足两个非一致性关系的最大匹配块构成的单元；

删除结构变异：变异区间内有一个非一致性单元，非一致性单元在测序序列上的位置相较于非一致性单元在参考基因组上的位置偏差长度大于-50bp；

重复结构变异：构成的非一致性单元在参考序列上存在大于20bp的重叠区域；

倒位结构变异：在框架上最大匹配块与参考基因组比对方向相反，但是相对于参考基因组是一致性单元或者是非一致性单元；

易位结构变异：发生在两条染色体上且相对于参考基因组个体片段位置发生交换。

进一步地，所述步骤五中的对步骤四获得的结构变异进行聚类，将同一类型的结构变异进行合并获得每个类型的结构变异簇；

步骤五一、对步骤四获得的结构变异进行聚类，将同一类型的结构变异放到一个集合中；

步骤五二、将每个集合中的结构变异进行排序，并根据排序后的结构变异获得结构变异簇：

在每个集合中将染色体一致的结构变异放在一起，然后按照染色体起始位置从小到大进行排序；

所述结构变异簇为结构变异排序后产生的位置重叠的结构变异；

步骤五三、在每个集合中对步骤五二获得结构变异簇进行合并，获得每个类型的结构变异簇：

删除结构变异、重复结构变异、倒位结构变异三种类型变异簇满足以下条件进行合并：

对于两个同一类型的结构变异簇，两个结构变异簇来自同一条染色体，且结构变异簇的位置中心坐标之间的距离不超过预先设定的中心坐标阈值，重叠片段的比值不小于预先设定的重叠片段比值阈值；

删除结构变异、重复结构变异、倒位结构变异三种类型变异簇通过三元组的形式记录，三元组内包括：染色体信息以及起始、终止位置信息；

插入结构变异簇满足以下条件进行合并：

两个结构变异簇来自同一条染色体，并且起始位置差值小于预先设定的第一起始坐标阈值，插入片段长度差值小于预先设定的长度阈值；

易位结构变异簇满足以下条件进行合并：

染色体A与染色B值相同，同时两组染色体的起始位置差值小于预先设定的第二起始坐标阈值2000bp；

所述染色体A和染色体B为两个变异结构簇所属的染色体；

所述易位结构变异簇以四元组的形式记录，四元组中的四个元素分别为：染色体A，染色体A变异起始位置，染色体B，染色体B变异起始位置。

进一步地，所述步骤六中的根据步骤五获取每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列，包括以下步骤：

步骤六一、利用每个类型的结构变异簇获得变异位点：

对于同一类型的结构变异簇，首先计算每个结构变异簇中的所有结构变异坐标的平均值，选择其中变异坐标与均值偏差小于预设第一偏差阈值的结构变异，即为变异位点；

步骤六二、获取结构变异的变异长度；

步骤六三、获取结构变异的碱基序列：

将步骤六一得到的结构变异簇的位点比对到参考基因组上获得结构变异的碱基序列。

进一步地，所述步骤六二中的获取结构变异的变异长度，具体为：

对于删除结构变异、插入结构变异以及重复结构变异：用终止坐标减去起始坐标获得变异长度；

对于倒位结构变异，采用以下公式获得变异长度：

其中，N表示每个类型的结构变异集合中结构变异的数量，Li表示倒位结构变异的长度。

本发明的有益效果为：

本发明提出的TTSV框架(本发明的方案即为TTSV框架)有针对性的对人类序列进行三代结构变异检测，本发明提出的TTSV省略了传统检测方法中序列比对的文件格式转换、片段信息排序、数据索引建立过程，减少了结构变异检测的计算量，进而提升了三代测序数据的速度。同时本发明在不损失准确性和敏感性的前提下，减少了变异检测的时间，提高了整个检测流程的效率。

附图说明

图1为本发明示意图。

具体实施方式

具体实施方式一：如图1所示，本实施方式一种基于比对框架的三代测序数据结构变异检测方法具体过程为：

步骤一、获取测序序列，利用测序序列获取精确匹配片段：

步骤一一、获取测序序列，按照预设步长K＝20bp划分每一条测序序列，划分为多个长度20bp的k-mer，利用每个长度为20bp的k-mer片段与预先通过测序序列构建的GRM图索引结构进行比对，获取k-mer片段中每一个碱基都能够完全比对成功的候选种子片段；

所述建立GRM图索引结构的过程即为测序序列片段与标准NCBI上能够下载的hg19人类基因组参考序列比对的过程；

步骤一二、对候选种子片段的上下游进行扩展，能够比对成功就继续延伸，直到上下游都遇到无法匹配的碱基为止，获得精确匹配片段；

步骤二、采用Landau-Vishkin模糊比对方法对每一个精确匹配片段继续进行延伸，直到比对不上碱基的数量超过我们事先设定的错误匹配阈值，结束延伸，称此时延伸得到的片段为模糊匹配片段。

Landau-Vishkin对于精确匹配片段的延伸发生在GRM图索引结构中。我们预先设定不匹配碱基数量阈值为3，对于每个精确匹配片段，我们记录一个五元组：染色体编号，参考基因组起始位置，参考基因组终止位置，测序片段的起始位置，测序片段的终止位置。

步骤三、采用稀疏动态规划算法利用步骤二获得的模糊匹配片段构建比对框架，包括以下步骤：

步骤三一、根据步骤二获得的模糊匹配片段的位置信息进行排序，先按照染色体编号从小到大排序。同一条染色体的模糊匹配片段按照参考基因组起始位置进行排序。

步骤三二、利用排序后模糊匹配片段构建无回路的有向图。我们将模糊匹配片段作为图的顶点，如果两个顶点位于同一条染色体上，且参考基因组上位置偏差和测序片段上的位置偏差小于预设最大跨度阈值，则构成一条有向边。方向是从参考基因组起始位置小的顶点指向起始位置大的顶点。

步骤三三、对于步骤三二构建的无回路有向图上每个顶点进行打分，打分基于以下公式：

S(Aj)＝max{S(Ai)+W+P}

其中，Aj和Ai表示有向边上的两个顶点j和i，其中i是j的前驱节点，W是顶点的权重(连接边的数量)，公式中表示Aj锚点顶点的权重值。Aj可能存在多个前驱节点，我们通过上述递推方程遍历图中每个顶点，给每个点给出一个得分。其中起始位置最小的点由于没有前驱节点，得分设置为0。

得分是一个递推关系，第一个P通过两个顶点参考序列以和测序序列之间偏差值的最小值计乘上测序错误率得出。

步骤三四、利用启发式动态规划算法按照分数从高到低依次选择顶点组成最优路径。最优路径可能存在多条，一条或多条最优路径组成比对框架；

其中，构造过程中如果存在多条最优比对骨架，TTSV(本发明)最多输出其中5条。

步骤四、对步骤三获得的比对框架上的模糊匹配片段进行扩展获取最大匹配块，通过最大匹配块一致性或者是非一致性关系解析其中蕴含的结构变异信息，包括以下步骤：

步骤四一、对步骤三构建的比对框架中顶点(模糊匹配片段)的上下游进行扩展获取最大匹配块，并获取最大匹配块之间的一致性关系，具体为：

我们使用固定大小的滑动窗口在两个顶点之间进行扩展；

步骤四一一、在比对框架中按照顶点长度从大到小依次选取顶点Ci作为延伸目标；

步骤四一二、预设窗口长度w，获取预设窗口w内顶点Ci的扩展比对片段在参考基因组、测序序列上的坐标，并根据坐标填补顶点Ci与顶点Cj之间的孔隙，获得最大匹配块及最大匹配块之间的一致性关系：

S1，将顶点Ci和与其下游相邻顶点Cj之间的间隙按照预设窗口长度w分割为多个窗口区域；

S2，利用局部哈希结构将窗口区域内的测序序列片段映射到参考基因组上，此时来自测序序列上能够比对到参考基因组上的片段称其为Ci在窗口内的扩展比对片段A，记录A在参考基因组上的坐标；

S3，根据A在参考基因组上的坐标和A在测序序列片段上的坐标获得A在参考基因组和测序序列片段的位置偏差，如果A在参考基因组和测序序列片段的位置偏差小于距离阈值15bp，则将扩展比对片段A填补到Ci与Cj之间的孔隙中：

S301，首先，获取A在参考基因组和测序序列片段的位置偏差：

将A在参考基因组上位置坐标与Ci在参考基因组上的位置坐标相减获得位置偏差a；

将A在测序序列上位置坐标与Ci在测序序列上的位置坐标相减获得位置偏差b；

获取b-a的值，若b-a<距离阈值15bp，则将扩展比对片段A填补到Ci与Cj之间的孔隙中，执行S302；

S302，将扩展比对片段加入顶点Ci并更新Ci在参考基因组上的坐标信息和Ci在测序序列上的坐标信息，如果是上游扩展，取A在参考基因组上的起始位置作为Ci在参考基因组上的起始位置，A在测序序列上的起始位置作为Ci在测序序列上的终止位置，如果是下游扩展，则取A在参考基因组上的终止坐标作为Ci在参考基因组上的起始位置，A在测序序列上的终止坐标作为Ci在测序序列上的终止坐标；

所述上游扩展为：A在参考基因组上的坐标小于更新后的Ci在参考基因组上的坐标；

A在测序序列上的坐标小于更新后的Ci在测序序列上的坐标；

所述下游扩展为：A在参考基因组上的坐标大于更新后的Ci在参考基因组上的坐标，A在测序序列上的坐标大于更新后Ci在测序序列上的坐标；

S4、按照固定长度w选取下一个窗口，重复窗口选取直到窗口w与Cj顶点相邻，获得当前加入扩展比对片段后的顶点(模糊匹配片段)称之为最大匹配块。通过上述步骤尽可能多的填补相连顶点之间的空隙，如果不同顶点形成的最大匹配块有交叠则称其满足一致性关系，并最大匹配块。如果两个最大匹配块无重叠，则称其满足非一致性关系。

步骤四二、从满足非一致性关系的最大匹配块中解析出结构变异：

插入结构变异(Insertion变异)是参考基因组上冗余的序列片段。对于此类序列，变异区间内会形成一个非一致性单元。同时非一致性单元在测序序列上的位置相较于非一致性单元在参考基因组序列上的位置偏差长度大于+50bp；

所述非一致性单元为满足两个非一致性关系的最大匹配块构成的单元；

其中，+50bp的含义为测序序列减去参考基因组序列等于50bp；

删除结构变异(Delection变异)是指参考基因组序列片段缺失的现象。对于此类序列，变异区间内会形成一个非一致性单元，长度偏差大于-50bp；

其中，-50bp含义为测序序列减去参考基因组序列等于-50bp，即参考基因序列减去测序序列等于50bp；

重复结构变异(Duplication变异)是一种特殊的插入变异，是指参考基因组上的冗余序列，这些序列又有较高的重复性。构成的非一致性单元在参考序列上存在大于20bp的重叠区域；

倒位结构变异(Inversion变异)是指相对参考基因组序列存在序列的反转。倒位在框架上最大匹配块比对方向恰好是相反的，但是相对于参考基因组可以是一致性单元或者是非一致性单元；

易位结构变异(Translocation变异)是指相对于参考基因组序列，个体片段交换位置的现象。易位类似于插入结构变异，但是相较于插入变异，易位发生在两条染色体上；

步骤五、对步骤四解析出的结构变异聚类，然后对同一类型结构变异进行合并得到每个类型的结构变异簇：

步骤五一、对步骤四解析出的结构变异进行聚类，同一类型的结构变异放到一个集合中；

步骤五二、对于分类的结构变异，我们需要进行排序操作。我们按照染色体和染色体起始位置作为排序的主要依据。染色体一样的结构变异放到一起，然后通过结构变异染色体起始位置从小到大进行排序操作，获得结构变异簇；

所述结构变异簇为排序后产生位置重叠的结构变异；

步骤五三、在每个集合中对步骤五二获得的结构变异簇进行合并，获得每个类型的结构变异簇：

删除结构变异、重复结构变异、倒位结构变异三种类型结构变异簇，都需要维护一个三元组信息。三元组内包含染色体信息以及起始、终止位置信息，记录为Sig＝(chr,s,e)。对于两个同一类型的结构变异簇，必须满足两个结构变异簇来自同一条染色体，位置中心坐标((起始位置+终止位置)/2)不超过预先设定的中心坐标阈值200，重叠片段的比值不小于预先设定的重叠片段比值阈值0.7；

插入结构变异簇对于同一组的两个结构变异簇需要保证来自同一条染色体，并且起始位置差值需要小于预先设定的第一起始坐标阈值50，插入片段长度差值需要小于预先设定的长度阈值20；

易位结构变异簇需要维护一个四元组信息，四个元素分别是染色体A，染色体A变异起始位置，染色体B，染色体B变异起始位置，记录为Sig＝(chrA,SA,chrB,SB)。我们需要保证统一组的两个结构变异簇，染色体A与染色B值相同。同时两组染色体的起始位置差值小于预先设定的第二起始坐标阈值2000bp。

步骤六、根据步骤五获取每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列，包括以下步骤：

步骤六一、利用每个类型的结构变异簇获得变异位点：

对于同一类型的结构变异簇，我们首先计算每个结构变异簇中的所有结构变异坐标的平均值，选择其中变异坐标与均值偏差小于预设第一偏差阈值的变异信号，即为变异位点；

步骤六二、获取结构变异的变异长度：

对于删除结构变异，插入结构变异以及重复结构变异我们只需要用终止坐标减去起始坐标就可以得到长度信息；

对于倒位结构变异，处理起来相对复杂，我们需要利用如下公式进行计算长度：

其中，N表示每个类型的结构变异簇中结构变异的数量，Li表示倒位结构变异的长度；

步骤六三、获取结构变异的碱基序列：

我们依据步骤六一得到的变异簇的位点位置信息，比对到参考基因组上得到最终的变异序列信息。

Claims (10)

1.一种基于比对框架的三代测序数据结构变异检测方法，其特征在于所述方法具体过程为：

步骤一、获取测序序列，利用测序序列获取精确匹配片段：

步骤一一、获取测序序列，按照预设步长K划分每条测序序列，利用每个长度为k的k-mer片段与预先建立的GRM图索引结构进行对比，获得每一个碱基都能够完全比对成功的k-mer片段作为候选种子片段；

所述GRM图索引结构通过测序序列建立；

步骤一二、对候选种子片段的上下游进行扩展直到上下游都遇到无法匹配的碱基为止，获得精确匹配片段；

步骤二、采用Landau-Vishkin模糊比对方法对精确匹配片段进行延伸，直至匹配错误的碱基数量大于预设错误匹配阈值时结束延伸，获得模糊匹配片段；

采用Landau-Vishkin模糊比对方法对每一个精确匹配片段进行延伸是在GRM图索引结构中进行的；

所述精确匹配片段通过五元组形式记录，具体为：染色体编号，在参考基因组上的起始位置，在参考基因组上的终止位置，在测序片段上的起始位置，在测序片段上的终止位置；

步骤三、采用稀疏动态规划算法利用步骤二获得的模糊匹配片段构建比对框架；

步骤四、对步骤三获得的比对框架上的模糊匹配片段进行扩展获取最大匹配块，对最大匹配块的一致性关系解析，获得结构变异；

所述一致性关系包括：一致性关系、非一致性关系；

步骤五、对步骤四获得的结构变异进行聚类，将同一类型的结构变异进行合并获得每个类型的结构变异簇；

步骤六、根据步骤五获取每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列。

2.根据权利要求1所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤三中的采用稀疏动态规划算法利用步骤二获得的模糊匹配片段构建比对框架，包括以下步骤：

步骤三一、将步骤二获得的模糊匹配片段进行排序：

首先，按照染色体编号将模糊匹配片段从大到小进行排序；

然后，将染色体编号相同的模糊匹配片段按照基因组起始位置的顺序进行排序；

步骤三二、利用排序后的模糊匹配片段构建无回路的有向图；

步骤三三、对步骤三二获得的无回路的有向图上的每个顶点进行打分；

步骤三四、利用启发式动态规划算法按照分数从高到低依次选择顶点组成最优路径，一条或多条最优路径组成比对框架。

3.根据权利要求2所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤三二中的利用排序后的模糊匹配片段构建无回路的有向图，具体为：

将模糊匹配片段作为顶点，如果两个顶点位于同一条染色体上，且在参考基因组上的位置偏差和在测序片段上的位置偏差均小于预设最大跨度阈值，则两个顶点之间的连线构成一条有向边；

所述有向边的方向从在参考基因组上的起始位置值小的顶点指向起始位置值大的顶点。

4.根据权利要求3所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤三三中的对步骤三二获得的无回路的有向图上的每个顶点进行打分，如下式：

S(Aj)＝max{S(Ai)+W-P}

其中，Aj和Ai表示连接后的边上的两个锚点j和i，W是锚点Aj的权重值，P表示Ai至Aj的罚分。

5.根据权利要求4所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤四中的对步骤三获得的比对框架上的模糊匹配片段进行扩展获取最大匹配块，对最大匹配块的一致性关系解析，获得结构变异，包括以下步骤：

步骤四一、对步骤三构建的比对框架中顶点的上下游进行扩展，获取最大匹配块，并获取最大匹配块之间的一致性关系，具体为：

步骤四一一、在比对框架中按照顶点长度从大到小依次选取顶点Ci作为延伸目标；

步骤四一二、预设窗口w，获取预设窗口w内顶点Ci的扩展比对片段A在参考基因组和测序序列上的坐标，并根据获取的坐标填补顶点Ci与顶点Cj之间的孔隙，获得最大匹配块及最大匹配块之间的一致性关系；

其中，Cj是Ci下游相邻的顶点；

所述最大匹配块之间的一致性关系包括：一致性关系、非一致性关系；

步骤四二、从满足不一致性关系的最大匹配块中解析出结构变异。

6.根据权利要求5所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤四一二中的预设窗口w，获取预设窗口w内顶点Ci的扩展比对片段A在参考基因组和测序序列上的坐标，然后利用A的坐标填补顶点Ci与Cj之间的孔隙，获得最大匹配块及最大匹配块之间的一致性关系，包括以下步骤：

S1、将顶点Ci和顶点Cj之间的间隙按照预设窗口长度w分割为多个窗口区域；

S2、利用哈希结构将窗口w内的测序序列片段映射到参考基因组上，获得能够与参考基因组比对成功的测序序列片段即顶点Ci在窗口w内的扩展比对片段A，获得A在参考基因组上的坐标；

S3、获取A在测序序列上的坐标，根据A在参考基因组上的坐标和A在测序序列片段上的坐标获得A在参考基因组和测序序列片段的位置偏差，如果A在参考基因组和测序序列片段的位置偏差小于预设距离阈值，则将A填补到Ci与Cj之间的孔隙中，获得最大匹配块：

S031、首先，获取A在参考基因组和测序序列片段的位置偏差：

将A在参考基因组上位置坐标与Ci在参考基因组上的位置坐标相减获得位置偏差a；

将A在测序序列上位置坐标与Ci在测序序列上的位置坐标相减获得位置偏差b；

获取b-a的值，若b-a<距离阈值15bp，则将扩展比对片段A填补到Ci与Cj之间的孔隙中，执行S302；

S302，将扩展比对片段A加入顶点Ci并更新Ci在参考基因组上的坐标信息和Ci在测序序列上的坐标信息，如果是上游扩展，取A在参考基因组上的起始位置作为Ci在参考基因组上的起始位置，A在测序序列上的起始位置作为Ci在测序序列上的终止位置，如果是下游扩展，则取A在参考基因组上的终止坐标作为Ci在参考基因组上的起始位置，A在测序序列上的终止坐标作为Ci在测序序列上的终止坐标；

所述上游扩展为：A在参考基因组上的坐标小于更新后的Ci在参考基因组上的坐标；

A在测序序列上的坐标小于更新后的Ci在测序序列上的坐标；

所述下游扩展为：A在参考基因组上的坐标大于更新后的Ci在参考基因组上的坐标；A在测序序列上的坐标大于更新后Ci在测序序列上的坐标；

S4，按照长度w选取下一个窗口，重复执行S2-S3直至选取的窗口与顶点Cj相邻，获得当前加入A后的顶点即为最大匹配块；若不同顶点形成的最大匹配块之间有交叠，则最大匹配块满足一致性关系；若不同顶点形成的最大匹配块之间没有重叠，则最大匹配块满足非一致性关系。

7.根据权利要求6所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤四二中的从满足不一致性关系的最大匹配块中解析出结构变异，具体如下：

插入结构变异：变异区间内有一个非一致性单元，同时非一致性单元在测序序列上的位置相较于非一致性单元在参考基因组上的位置偏差长度大于+50bp；

所述非一致性单元为满足两个非一致性关系的最大匹配块构成的单元；

删除结构变异：变异区间内有一个非一致性单元，非一致性单元在测序序列上的位置相较于非一致性单元在参考基因组上的位置偏差长度大于-50bp；

重复结构变异：构成的非一致性单元在参考序列上存在大于20bp的重叠区域；

倒位结构变异：在框架上最大匹配块与参考基因组比对方向相反，但是相对于参考基因组是一致性单元或者是非一致性单元；

易位结构变异：发生在两条染色体上且相对于参考基因组个体片段位置发生交换。

8.根据权利要求7所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤五中的对步骤四获得的结构变异进行聚类，将同一类型的结构变异进行合并获得每个类型的结构变异簇；

步骤五一、对步骤四获得的结构变异进行聚类，将同一类型的结构变异放到一个集合中；

步骤五二、将每个集合中的结构变异进行排序，并根据排序后的结构变异获得结构变异簇：

在每个集合中将染色体一致的结构变异放在一起，然后按照染色体起始位置从小到大进行排序；

所述结构变异簇为结构变异排序后产生的位置重叠的结构变异；

步骤五三、在每个集合中对步骤五二获得结构变异簇进行合并，获得每个类型的结构变异簇：

删除结构变异、重复结构变异、倒位结构变异三种类型变异簇满足以下条件进行合并：

对于两个同一类型的结构变异簇，两个结构变异簇来自同一条染色体，且结构变异簇的位置中心坐标之间的距离不超过预先设定的中心坐标阈值，重叠片段的比值不小于预先设定的重叠片段比值阈值；

删除结构变异、重复结构变异、倒位结构变异三种类型变异簇通过三元组的形式记录，三元组内包括：染色体信息以及起始、终止位置信息；

插入结构变异簇满足以下条件进行合并：

两个结构变异簇来自同一条染色体，并且起始位置差值小于预先设定的第一起始坐标阈值，插入片段长度差值小于预先设定的长度阈值；

易位结构变异簇满足以下条件进行合并：

染色体A与染色B值相同，同时两组染色体的起始位置差值小于预先设定的第二起始坐标阈值2000bp；

所述染色体A和染色体B为两个变异结构簇所属的染色体；

所述易位结构变异簇以四元组的形式记录，四元组中的四个元素分别为：染色体A，染色体A变异起始位置，染色体B，染色体B变异起始位置。

9.根据权利要求8所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤六中的根据步骤五获取每个类型的结构变异簇进行结构变异检测，获得变异位点、变异长度、变异序列，包括以下步骤：

步骤六一、利用每个类型的结构变异簇获得变异位点：

对于同一类型的结构变异簇，首先计算每个结构变异簇中的所有结构变异坐标的平均值，选择其中变异坐标与均值偏差小于预设第一偏差阈值的结构变异，即为变异位点；

步骤六二、获取结构变异的变异长度；

步骤六三、获取结构变异的碱基序列：

将步骤六一得到的结构变异簇的位点比对到参考基因组上获得结构变异的碱基序列。

10.根据权利要求9所述的一种基于比对框架的三代测序数据结构变异检测方法，其特征在于：所述步骤六二中的获取结构变异的变异长度，具体为：

对于删除结构变异、插入结构变异以及重复结构变异：用终止坐标减去起始坐标获得变异长度；

对于倒位结构变异，采用以下公式获得变异长度：

其中，N表示每个类型的结构变异集合中结构变异的数量，Li表示倒位结构变异的长度。

Images