CN115762633B - 一种基于三代测序的基因组结构变异基因型校正方法 - Google Patents

Abstract

一种基于三代测序的基因组结构变异基因型校正方法，涉及生物信息处理领域。本发明是为了解决现有结构变异检测方法还存在由于测序错误噪声、测序深度导致的基因变异检测结果存在误差的问题。本发明包括：获取个体测序数据比对片段，根据比对片段信息提取异常比对片段，并将其划分为比对内部信号和比对间信号，从异常比对片段中提取不同类型的结构变异信号，并获取所有比对片段的起始位置、终止位置、片段名称；获取待分型结构变异，利用待分型结构变异对结构变异信号进行筛选获得支持待分型结构变异的结构变异信号集；利用所有比对片段信息和支持待分型结构变异的结构变异信号集获取待分型结构变异的基因型。本发明用于对结构变异基因型进行校正。

Description

一种基于三代测序的基因组结构变异基因型校正方法

技术领域

本发明涉及生物信息处理领域，特别涉及一种基于三代测序的基因组结构变异基因型校正方法。

背景技术

基因虽然十分稳定，能在细胞分裂时精确地复制自己，但这种稳定性是相对的。在一定的条件下基因也可以从原来的存在形式突然改变成另一种新的存在形式，就是在一个位点上，突然出现了一个新基因，代替了原有基因，这就是基因变异。从分子水平上看，基因变异是指基因在结构上发生碱基对组成或排列顺序的改变。根据人体的基因组特性，人类基因组是二倍体，是由两个单倍体组成的，某个变异在二倍体上可能有三种出现方式：在两个单倍体上都包含此变异、只有一个单倍体上包含此变异、两条单倍体上都不包含此变异，分别对应基因型为“1/1”、“0/1”、“0/0”。对变异基因型的精准检测能够精确绘制变异在人群中的图谱，同时根据孟德尔遗传学定律可以进行家系上遗传、进化的深入分析。

目前研究中的基因型检测多与个体变异检测结合在一起，在检测个体SV(Structural Variations，基因组结构变异)的同时判断检测结构变异的基因型，但由于存在测序错误噪声、测序深度不足等情况，导致基因变异的基因型检测结果存在一定误差。同时，当将检测的对象扩大为群体时，会对个体检测中得到的结构变异信息进行一定的修正，使得在个体结构变异检测中得到的基因型信息不能够准确地应用到群体上。目前的结构变异基因型的校正方法中一般将基因型的校正与测序片段比对、个体结构变异检测结合到了一起，导致校正基因型结果的准确性与精确度仍有缺陷，同时，增加了不必要的运算量，大量增加了所花费的时间，导致了校正速度缓慢的问题。

发明内容

本发明目的是为了解决现有的基因型校正方法还存在准确率低、精度低、校正速度缓慢的问题，而提出了一种基于三代测序的基因组结构变异基因型校正方法。

一种基于三代测序的基因组结构变异基因型校正方法具体过程为：

步骤一、获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并从异常比对片段中提取不同类型的结构变异信号，获取提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称；

步骤二、获取待分型结构变异，并基于空间相似度利用待分型结构变异对步骤一获得的结构变异信号进行筛选获得支持待分型结构变异的结构变异信号集；

步骤三、利用比对片段的信息和结构变异信号集获取待分型结构变异的基因型；

所述比对片段信息是步骤一获得的提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称；

所述结构变异信号集为步骤二获得的支持待分型结构变异的结构变异信号集。

本发明的有益效果为：

本发明提出的基于局部筛选优化的结构变异基因型校正方法，面向异常比对片段进行结构变异信号的提取，通过度量结构变异信号之间的空间相似性，筛选支持待分型结构变异的信号，根据支持待分型结构变异信号的数量与测序时覆盖待分型结构变异的所有信号数量进行似然估计，根据似然估计值为待分型结构变异进行基因型校正。本发明通过对结构变异信号进行全面提取与精准筛选来进行基于最大似然估计的基因型计算校正，避免了噪声信号影响真正的结构变异信号，提高了计算精度和所得基因型的准确性；同时，本发明提出的结构变异基因型校正方法直接使用个体的比对结果作为输入，减少了重复比对所花费的大量时间，并利用分块并行策略，利用计算资源的多核特性，大幅提升基因型校正速度。

附图说明

图1为本发明的流程图。

具体实施方式

具体实施方式一：本实施方式一种基于三代测序的基因组结构变异基因型校正方法具体过程如图1所示，具体包括以下步骤：

步骤一、获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，从异常比对片段中提取不同类型的结构变异信号，并获取提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称；

步骤二、获取待分型结构变异，并基于空间相似度利用步骤一获得的结构变异信号对待分型结构变异进行筛选获得支持待分型结构变异的结构变异信号集；

步骤三、利用步骤一获得的提取不同类型的结构变异信号时比对片段的信息和步骤二获得的支持待分型结构变异的结构变异信号集获得待分型结构变异的基因型。

具体实施方式二：所述步骤一的获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并从异常比对片段中提取不同类型的结构变异信号，获取提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称，包括以下步骤：

步骤一一、获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并将异常比对片段按照比对形式划分为比对内部信号(Intra-alignment signatures)和比对间信号(Inter-alignment signatures)；

基因组测序片段一般可以表示为被测个体基因组局部区域的一段连续DNA序列。当被测个体基因组发生了一段DNA片段序列的改变，根据其DNA变化形式，相应的测序片段比对也会随之发生中断、冗余、片段丢失等现象，使得其在比对到参考基因组上时发生异常，这类片段比对结果称之为异常片段比对；

所述比对内部信号为个体测序数据的基因组测序片段面向参考基因组比对时，获得的不匹配信息以碱基级比对结果(CIGAR域)的形式记录在比对片段内部的信息；

所述比对间信号为个体测序数据的基因组测序片段面向参考基因组比对时，将个体测序数据的基因组测序片段拆分成多个小片段比对到参考基因组上的信息，其中一个小片段即为一个拆分比对；

所述比对内部信号和比对间信号划分方法如下：

在对测序片段进行面向参考基因组的比对时，会识别与测序片段最为相似的基因组区域，并将其比对上去，而由于基因组变异、测序错误等的存在，导致测序片段无法与参考序列进行完全一致的匹配，例如出现测序片段相比参考基因组增加或减少了一些碱基、测序片段在参考基因组上的比对出现了中断等等。为了更好的完成比对，比对软件会对比对情况计算一个惩罚分数，当遇到不匹配信息时，增加比对的惩罚分数，越多的不匹配信息会带来越大的惩罚分数，分数越高代表比对的情况越差。比对中的不匹配信息以碱基级比对结果(CIGAR域)的形式记录在比对片段内部，形成了比对内部的信号；而当比对中的不匹配信息累计增多，惩罚得分不断增大时，会导致比对情况变差，因此，当惩罚得分超过比对软件指定阈值时，将测序片段进行一定的拆分，拆分后的每段小片段会比对到基因组不同的区域，使得拆分后的片段有着相对更低的惩罚得分，形成了拆分比对。

步骤一二、根据步骤一一获得的异常比对片段的划分结果确定异常比对片段中是否存在比对间信号，若存在比对间信号则执行步骤一三与步骤一四，若不存在比对间信号则执行步骤一四；

步骤一三、从拆分比对中提取插入结构变异信号、删除结构变异信号、倒位结构变异信号、重复结构变异信号、易位结构变异信号：

步骤一三一、从拆分比对中提取插入结构变异信号、删除结构变异信号：

获取相邻两段拆分比对在参考基因组上的位置分别为A(s1,e1)、B(s2,e2)和同一相邻两段拆分比对在比对片段上的位置(s1’,e1’)以及(s2’,e2’)，当两段拆分比对在参考基因组上位置接近(|e1-s2|<a)，而在比对片段中间隔了较长的一段碱基(s2’–e1’>a)，可以认为这一段碱基(s2’–e1’)是一个插入变异，将这一段碱基的起始位置与长度记录为插入结构变异信号。当两段拆分比对在比对片段上位置接近(|e1’-s2’|<a)，而在参考基因组上间隔了较长的一段碱基(s2–e1>a)，可以认为参考基因组上这一段碱基(s2–e1)被删除了，将这一段碱基的起始位置与长度记录为删除结构变异信号；

其中，a是预设的长度阈值，为常数。

步骤一三二、从拆分比对中提取倒位结构变异信号：

倒位变异是指比对片段相对于参考基因组有一段碱基序列发生了反转，根据倒位变异的结构，由于其碱基序列发生了反转，在拆分比对中一定存在一个比对片段比对方向与其他比对片段比对方向相反，将比对方向相反的比对片段的起始与终止坐标提取出来记录为倒位结构变异信号。

步骤一三三、从拆分比对中提取重复结构变异信号：

重复变异是指比对片段相对于参考基因组多了一段碱基序列的重复，它是一种插入变异，但插入的序列与参考基因组上已存在的序列相同，由于它涉及到具体的序列信息，不能单独靠比对内信号提取，要使用比对间信号，当拆分比对中，位于被测个体测序数据的基因组测序片段上游的片段在参考基因组上比对到了其下游片段的后方，则参考基因组上以上游片段和下游片段为起始点和终点的这段序列发生了重复，将这段序列的起始位点与长度提取出来记录为重复结构变异信号。

步骤一三四、从拆分比对中提取易位结构变异信号：

易位变异是指比对片段相对于参考基因组有两段碱基序列之间交换了位置，由于它发生在不同的染色体之间，因此只能通过拆分比对来提取对应信号，当相邻两段拆分比对被比对到不同的染色体上时，代表这里发生了染色体间的易位，将两段拆分比对在参考基因组上的染色体与坐标位点提取出来记录为易位结构变异信号。

步骤一四、提取比对片段的CIGAR域，并从CIGAR域中提取插入结构变异信号、删除结构变异信号：

CIGAR域是一组以数字+字母组成的字符串，其中不同的字母代表不同的比对方式，例如“8＝2I4X5D”代表这条比对片段中，有8个碱基比对到了参考基因组上，接下来的2个碱基是参考基因组上没有的而比对片段上出现的碱基，即插入片段，接下来4个碱基与参考基因组匹配失败，接下来5个碱基是参考基因组上存在而在比对片段上不存在的碱基，即删除片段。在提取结构变异信号时，我们关注CIGAR域中的I与D信息，将插入片段与删除片段的信息提取出来，即插入或删除片段的长度，同时从比对文件中提取该比对片段比对到的染色体、在参考基因组上出现的位置，将这三者信息以三元组的形式记录为插入结构变异信号或删除结构变异信号；

所述比对文件为个体测序数据的测序片段与参考基因组比对结果的集合。

步骤一五、获取步骤一三到步骤一四提取所有不同类型的结构变异信号时比对片段的起始位置、终止位置、片段名称，并将信息写入磁盘临时文件。

具体实施方式三：所述步骤二中的获取待分型结构变异，并基于空间相似度利用步骤一获得的结构变异信号对待分型结构变异进行筛选，获得支持待分型结构变异的结构变异信号集：

步骤二一、获取待分型结构变异，并按照待分型结构变异所在染色体对待分型结构变异进行划分，每条染色体上的待分型结构变异组成一个待分型结构变异集合；

步骤二二、从步骤一获取的结构变异信号中提取与每个待分型结构变异集合染色体相关的结构变异信号作为每个待分型结构变异集合对应的待筛选信号集合；

步骤二三、对每一个待分型结构变异，获得待筛选信号集合中首个起始位点大于当前待分型结构变异位置的结构变异信号：

在筛选过程中，由于结构变异信号是按照基因组坐标有序的，因此采用二分查找的方法，在较低的时间复杂度内搜索得到结构变异信号中首个起始位点恰好大于待分型结构变异坐标的信号，然后在其上下游扩展，筛选符合支持要求的结构变异信号。具体来说，在二分查找中，搜索区间初始化为待筛选信号集合中的结构变异信号，每次取区间中位数位置的结构变异信号，比较其基因组坐标与待分型结构变异的基因组坐标，当其基因组坐标小于待分型结构变异的基因组坐标时，说明待分型结构变异基因组坐标大于此中位数结构变异信号基因组坐标，因此将搜索区间减半，缩小为后半段区间；当其基因组坐标大于或等于待分型结构变异的基因组坐标时，说明待分型结构变异基因组坐标不大于此中位数结构变异基因组坐标，因此将搜索区间缩小为前半段区间，直至搜索区间大小为1，找到结构变异信号中首个起始位点恰好大于待分型结构变异基因组坐标的信号位置。

步骤二四、根据步骤二三中获得的待筛选信号集合中首个起始位点大于待分型结构变异的结构变异信号，在其上下游进行扩展，在待筛选信号集合中筛选出支持待分型结构变异的所有结构变异信号，从而获得支持待分型结构变异的结构变异信号集，支持待分型结构变异的所有结构变异信号，满足如下公式：

其中SV_pos、SV_len代表待分型结构变异的基因组坐标位点与长度，sig_pos、sig_len代表待筛选信号集合中结构变异信号的基因组坐标位点与长度；

在任一待分型结构变异对应的待筛选结构变异信号集合的上下游进行扩展，满足如下式公式停止扩展：

其中，sig′_pos是当前结构变异信号上/下游临近的结构变异信号的基因组坐标。

具体实施方式四：所述步骤三中的利用步骤一获得的所有比对片段的信息和步骤二获得的支持待分型结构变异的结构变异信号集计算待分型结构变异的基因型：

步骤三一、利用步骤一获得的提取不同类型的结构变异信号时比对片段的比对起始位置、终止位置、片段名称信息获取覆盖待分型结构变异的所有比对片段，并统计比对片段数量：

比对片段完全覆盖结构变异满足以下公式：

其中read_s、read_e分别代表比对片段在参考基因组上的起始、终止位置，SV_pos、SV_end分别代表结构变异在参考基因组上的起始、终止位置；

当比对片段完全覆盖结构变异时将其记录下来并统计数量。

步骤三二、将步骤三一获得的覆盖待分型结构变异的所有比对片段与步骤二获得的支持待分型结构变异的结构变异信号所来源的比对片段做差集获得不支持待分型结构变异的结构信号集合，利用支持待分型结构变异的结构变异信号集合和不支持待分型结构变异的结构变异信号集合计算待分型结构变异不同基因型的似然概率：

其中，DR代表不支持待分型结构变异的比对片段数量，DV代表支持待分型结构变异的比对片段数量，ε代表比对片段被错误地比对到单倍体的概率。

步骤三三、获取步骤三二获得的待分型结构变异三种基因型的似然概率中最大的似然概率对应的基因型即待分型结构变异的基因型。

Claims (7)

1.一种基于三代测序的基因组结构变异基因型校正方法，其特征在于所述方法具体过程为：

步骤一、获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并从异常比对片段中提取不同类型的结构变异信号，获取提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称；

步骤二、获取待分型结构变异，并基于空间相似度利用待分型结构变异对步骤一获得的结构变异信号进行筛选获得支持待分型结构变异的结构变异信号集；

步骤三、利用比对片段的信息和结构变异信号集获取待分型结构变异的基因型；

所述比对片段信息是步骤一获得的提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称；

所述结构变异信号集为步骤二获得的支持待分型结构变异的结构变异信号集；

步骤三包括以下步骤：

步骤三一、利用步骤一获得的提取不同类型的结构变异信号时比对片段的比对起始位置、终止位置、片段名称信息获取覆盖待分型结构变异的所有比对片段，并统计比对片段数量；

获取覆盖待分型结构变异的所有比对片段，需满足以下公式：

其中read_s、read_e分别代表异常比对片段在参考基因组上的起始、结束位置，SV_pos、SV_end分别代表结构变异在参考基因组上的起始、结束位置；

步骤三二、将步骤三一获得的覆盖待分型结构变异的所有比对片段与步骤二获得的支持待分型结构变异的结构变异信号所来源的比对片段做差集获得不支持待分型结构变异的结构信号集合，利用支持待分型结构变异的结构变异信号集合和不支持待分型结构变异的结构变异信号集合计算待分型结构变异不同基因型的似然概率，如下式：

其中，DR代表不支持待分型结构变异的比对片段数量，DV代表支持待分型结构变异的比对片段数量，ε代表比对片段被错误地比对到单倍体的概率；

步骤三三、获取步骤三二获得的待分型结构变异不同基因型的似然概率中最大的似然概率对应的基因型即待分型结构变异的基因型。

2.根据权利要求1所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤一中的获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并从异常比对片段中提取不同类型的结构变异信号，获取提取不同类型的结构变异信号时比对片段的起始位置、终止位置以及片段名称，包括以下步骤：

步骤一一、获取个体测序数据的比对片段，根据比对片段的信息提取包含结构变异的异常比对片段，并将异常比对片段划分为比对内部信号和比对间信号；

所述比对内部信号为个体测序数据的基因组测序片段面向参考基因组比对时，获得的不匹配信息，不匹配信息以CIGAR域的形式记录在比对片段内部；

所述比对间信号为个体测序数据的基因组测序片段面向参考基因组比对时，将个体测序数据的基因组测序片段拆分成多个小片段比对到参考基因组上的信息，其中一个小片段即为一个拆分比对；

步骤一二、根据步骤一一获得异常比对片段的划分结果确定异常比对片段中是否存在比对间信号，若存在比对间信号则执行步骤一三、步骤一四，若不存在比对间信号则执行步骤一四；

步骤一三、从拆分比对中提取插入结构变异信号、删除结构变异信号、倒位结构变异信号、重复结构变异信号、易位结构变异信号；

步骤一四、提取比对片段的CIGAR域，并从CIGAR域中提取插入结构变异信号、删除结构变异信号；

步骤一五、获取步骤一三到步骤一四提取所有不同类型的结构变异信号时比对片段的起始位置、终止位置、片段名称。

3.根据权利要求2所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤一三中的从拆分比对中提取插入结构变异信号、删除结构变异信号、倒位结构变异信号、重复结构变异信号、易位结构变异信号，具体包括以下步骤：

步骤一三一、从拆分比对中提取插入结构变异信号和删除结构变异信号：

首先，获取相邻两段拆分比对m和n在参考基因组上的位置A(s1,e1)、B(s2,e2)，以及m和n在比对片段上的位置A’(s1’,e1’)以及B’(s2’,e2’)；

然后，利用A、B、A’、B’获取插入结构变异信号和删除结构变异信号：

若|e1-s2|<a且s2’-e1’>a，则s2’-e1’段碱基是一个插入变异，将这一段碱基的起始位置与长度记录为插入结构变异信号；

若|e1’-s2’|<a且s2-e1>a，则参考基因组上s2-e1段碱基被删除了，将这一段碱基的起始位置与长度记录为删除结构变异信号；

其中，a是长度阈值；

步骤一三二、从拆分比对中提取倒位结构变异信号：

在拆分比对中存在一个比对片段比对方向与其他比对片段比对方向相反，将比对方向相反的比对片段的起始与终止坐标提取出来记录为倒位结构变异信号；

步骤一三三、从拆分比对中提取重复结构变异信号：

若拆分比对中位于被测个体测序数据的基因组测序片段上游的片段在参考基因组上比对到了其下游片段的后方，则参考基因组上以上游片段和下游片段为始终的这段序列为重复序列，将这段重复序列的起始位点与长度提取出来记录为重复结构变异信号；

步骤一三四、从拆分比对中提取易位结构变异信号：

当相邻两段拆分比对被比对到不同的染色体上时，代表这里发生了染色体间的易位，将两段拆分比对在参考基因组上的染色体与坐标位点提取出来记录为易位结构变异信号。

4.根据权利要求3所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤一四中的提取比对片段的CIGAR域，并从CIGAR域中提取插入结构变异信号、删除结构变异信号，具体如下：

首先，根据CIGAR域中I和D的信息获得比对片段中插入和删除片段的长度：

其中，CIGAR域中I前面的数字即为插入片段的长度，D前面的数字即为删除片段的长度；

然后，从比对文件中提取比对片段比对到的染色体、以及比对片段在参考基因组上出现的位置；

所述比对文件为个体测序数据的测序片段与参考基因组比对结果的集合；

最后，将删除或插入片段的长度、比对片段比对到的染色体、比对片段在参考基因组上出现的位置以三元组的形式记录为插入结构变异信号或删除结构变异信号。

5.根据权利要求4所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤二中的获取待分型结构变异，并基于空间相似度利用待分型结构变异对步骤一获得的结构变异信号进行筛选获得支持待分型结构变异的结构变异信号集，包括以下步骤：

步骤二一、获取待分型结构变异，并按照待分型结构变异所在染色体对待分型结构变异进行划分，每条染色体上的待分型结构变异组成一个待分型结构变异集合；

步骤二二、从步骤一获取的结构变异信号中提取与每个待分型结构变异集合染色体相关的结构变异信号作为每个待分型结构变异集合对应的待筛选信号集合；

步骤二三、对每一个待分型结构变异，获得待筛选信号集合中首个起始位点大于当前待分型结构变异基因组坐标的结构变异信号；

步骤二四、在待筛选信号集合中首个起始位点大于当前待分型结构变异基因组坐标的结构变异信号的上下游进行扩展，在待筛选信号集合中筛选出支持待分型结构变异的所有结构变异信号，从而获得支持待分型结构变异的结构变异信号集。

6.根据权利要求5所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤二三中的对每一个待分型结构变异，获得待筛选信号集合中首个起始位点大于当前待分型结构变异基因组坐标的结构变异信号，具体如下：

将搜索区间初始化为待筛选信号集合，每次取搜索区间中位数位置的结构变异信号，比较中位数位置的结构变异信号基因组坐标与待分型结构变异的基因组坐标，当结构变异信号基因组坐标小于待分型结构变异基因组坐标时，说明待分型结构变异基因组坐标大于此中位数结构变异信号的基因组坐标，因此将搜索区间减半，缩小为后半段区间，再次查找；当结构变异信号基因组坐标大于或等于待分型结构变异基因组坐标时，说明待分型结构变异基因组坐标不大于此中位数结构变异信号基因组坐标，因此将搜索区间缩小为前半段区间再次查找，直至搜索区间大小为1，找到待筛选信号集中首个起始位点大于当前待分型结构变异基因组坐标的结构变异信号。

7.根据权利要求6所述的一种基于三代测序的基因组结构变异基因型校正方法，其特征在于：所述步骤二四中的在待筛选信号集合中筛选出支持待分型结构变异的所有结构变异信号，满足以下公式：

其中SV_pos、SV_len代表待分型结构变异的基因组坐标位点与长度，sig_pos、sig_len代表待筛选信号集合中结构变异信号的基因组坐标位点与长度；

所述对每一个待分型结构变异，在待筛选信号集合中首个起始位点大于当前待分型结构变异基因组坐标的结构变异信号的上下游进行扩展，满足以下公式停止扩展：

其中，sig′_pos是当前结构变异信号上/下游临近的结构变异信号的基因组坐标。