CN111445950B - 一种基于滤波策略的高容错基因组复杂结构变异检测方法 - Google Patents

Abstract

本发明公开了一种基于滤波策略的高容错基因组复杂结构变异检测方法，对SAM格式的输入文件进行预处理，遍历最优质量比对读段中的CIGAR字段；根据比对后的CIGAR字段和变异分数计算准则，计算出各个位点在当前读段对应的变异分数，并将其预先保存在每个位点的变异分数集合中；统计每个位点的变异分数集合中的平均数当作该位点最终的变异分数并得到此样本的变异分数函数；对变异分数函数进行卡尔曼或高斯滤波，得到滤波降噪后的变异分数函数；依照滤波后的变异分数函数，设定阈值并分离出结构变异区域，提取特征；训练支持向量机(SVM)模型，再用训练好的SVM模型对结构变异区域分类并得到复杂indel结果集。本发明解决测序错误对结构变异的确定产生的干扰。

Description

一种基于滤波策略的高容错基因组复杂结构变异检测方法

技术领域

本发明属于第三代核酸序列测序(Single Molecule Real Time，SMRT)技术领域，具体涉及一种基于滤波策略的高容错基因组复杂结构变异检测方法。

背景技术

复杂indel(Complex insertion-deletion)是一种在人群中相对罕见但在肿瘤基因组中较多存在的基因组结构变异。复杂indel表现为在某一基因上DNA片段发生了缺失变异，由于DNA分子的自我修复机制，随后在同一位点上插入了其他的DNA片段并且插入片段有可能发生倒置的一种复合变异。目前已发现的复杂indel的表现形式就有数十种。作为一种重要结构变异，复杂indel的检测是下游分析肿瘤易感性与表型相关性等研究的基础。扩大复杂indel检测范围不仅有助于确定复杂indel的基因型及其表型效应，而且也能够推动研究肿瘤复杂indel之间的关系，加快个性化医疗的步伐。

检测复杂indel主要通过基因组测序数据。目前已有的复杂indel检测算法主要有三种，分别是INDELseek、Pindel-C和SV-Bay。INDELseek对比对质量字段进行重译选择出其中变异位点在读段内进行聚类，将聚类出的变异区域按照阈值进行过滤选择出复杂indel。Pindel-C找出一端完美比对另一端不能完美比对的读对，使用模式增长算法和分裂split-read的思想，筛选出不能完美比对的读段，并将该读段作为一个复杂indel发生的可疑区域。SV-Bay主要利用贝叶斯方法识别结构变异。

但是，上述算法都只适用于第二代测序数据，第二代测序数据的单一读段长度仅有100bps。读段长度短也意味着读段携带的信息较少，会影响复杂indel的检测效果，第三代测序数据具有读段长度长的显著优势。然而，使用第三代测序数据来检测复杂indel，主要的阻碍来自于第三代测序数据的读段高错误率。高错误率导致传统算法在处理读段时，被迫丢弃大量的高错误率的读段，或对读段做出剪裁等，以期获得错误率低的读段数据。此时，传统算法主要有以下两方面的不足：

(1)在检测结构变异时，测序数据中的测序错误会对结构变异的识别产生干扰，导致算法误把本该忽略的测序错误当作结构变异，或者将本该识别出来的结构变异当作测序错误而忽略；

(2)在区分不同类型的结构变异时，如区分复杂indel和其他结构变异时，测序错误会干扰正确的分类。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于滤波策略的高容错基因组复杂结构变异检测方法，解决测序错误对结构变异的确定产生的干扰，以及解决测序错误对复杂indel和其他结构变异的区分产生的干扰。

本发明采用以下技术方案：

一种基于滤波策略的高容错基因组复杂结构变异检测方法，包括以下步骤：

S1、对SAM格式的输入文件进行预处理，遍历最优质量比对读段中的CIGAR字段；

S2、根据比对后的CIGAR字段和变异分数计算准则，计算出各个位点在当前读段对应的变异分数，并将其预先保存在每个位点的变异分数集合中；

S3、统计每个位点的变异分数集合中的平均数当作该位点最终的变异分数并得到此样本的变异分数函数；

S4、对变异分数函数进行卡尔曼或高斯滤波，得到滤波降噪后的变异分数函数；

S5、依照滤波后的变异分数函数，设定阈值并分离出结构变异区域，提取特征；

S6、训练支持向量机SVM模型，再用训练好的SVM模型对结构变异区域分类并得到复杂indel结果集。

具体的，步骤S2中，位点i的变异分数S_i为：

其中，多次计算位点i的变异分数并分别记作{S_i1,S_i2,...,S_ij}，j表示位点i的测序覆盖度，k为窗口半径，位点i变异分数依赖于(i-k,i+k)范围的变异程度；C_j(i)为逻辑函数，表示在第j个读段中位点i在某个具体的比对结果中的变异情况，如果位点i发生了变异，则C_j(i)的值为1，反之为0。

具体的，步骤S3中，变异分数函数即计算出的所有位点的变异分数为离散函数，能够描述染色体位点对于参考序列的变异程度，值域为[0,1]。

具体的，步骤S4中，变异分数函数卡尔曼滤波过程如下：

S4011、通过位点i-1的最优预估值S_i'_-1确定位点i的变异分数预测值

S4012、根据上一次计算的误差方差P_i-1和预测过程噪声ω，预测当前位点i的误差方差

S4013、结合测量系统参数H和测量噪声ε，计算当前位点i的卡尔曼增益K_i；

S4014、结合当前位点i的测量值y_i计算当前位点i的最优预估值S_i'；

S4015、更新误差方差，并重复迭代以上步骤。

进一步的，变异分数预测值

如下：

误差方差

卡尔曼增益K_i：

最优预估值S′_i：

其中，A，B是系统参数，U是系统控制量，P_i为误差方差，H为测量系统参数，I为单位矩阵，ε为测量噪声。

具体的，步骤S4中，变异分数函数高斯滤波过程如下：

S4021、取变异位点i，变异位点i坐标代表其本身变异与否，若变异取1，否则取0；取距离i点最近的k个坐标点，即与位点i相邻的k个位点，根据数据量和计算能力确定k的取值，对于普通计算机建议取值为k＝8；

S4022、设定μ和σ，计算出9个位点对应的高斯模板取值G_σj，j＝1,2,...,9，即位点i的邻域点的变异与否对位点i是否为伪变异的影响权重；

其中，G_σ是标准差为σ的高斯核；

S4023、对高斯模板归一化，即

得到最终的高斯模板；

S4024、计算9个位点的高斯滤波I_σj，并相加得到位点i的高斯滤波的值：

I_σ＝I*G_σ

其中，*是卷积操作。

进一步的，如果某个位点的变异为伪变异，伪变异及其周围的分数将通过滤波平滑掉；如果某个位点是结构变异中的位点，结构变异与伪变异区分开。

具体的，步骤S5中，设定阈值并分离结构变异区域具体如下：

将滤波后变异分数的均值和方差分别记作μ'和σ'，以μ'+3σ'为阈值，将变异分数函数和滤波后的变异分数函数进行上下分割并取出大于μ'+3σ'的区间，分别记作：

SV＝{(a₁,b₁),(a₂,b₂),...,(a_n,b_n)}

SV'＝{(A₁,B₁),(A₂,B₂),...,(A_n,B_n)}

将左端点A当作对SV集合的筛选的标准，对于(a,b)∈SV，如果

且A∈(a,b)，保留，反之则丢弃，筛选后SV集合中的每个元素代表一个结构变异区域。

进一步的，提取特征具体如下：

选取结构变异区域的插入比率、删除比率、替换比率、变异分数最大值以及四等分位点对应的变异分数作为特征值，以上七个特征值分别记作insRatio、delRatio、subRatio、score₀、score₁、score₂、score₃，分别给出区域(a,b)内的计算公式：

其中，T_ij、P_ij、Q_ij分别是第j个读段内目标序列，比对信息和读段序列的第i个字符的值，c为insert-size正常区域半径，I为指示函数，S_i是位点i的变异分数。

具体的，步骤S6中，训练SVM模型如下：

S601、将结构变异集合中每个变异看作高维空间中的一个点，结构变异的每个特征代表一个维度；

S602、用高斯核函数将结构变异集合中的元素映射到结构变异特征个数组成的高维空间中，然后在高维空间中训练得到复杂indel和其他结构变异的最优分离超平面；

S603、把得到的模型用于结构变异的复杂indel分类中。

与现有技术相比，本发明至少具有以下有益效果：

本发明是一种基于滤波策略的高容错基因组复杂结构变异的检测方法，遍历SAM格式的输入文件中所有比对结果，因为SAM文件已经给出了读段的所有比对质量并以此排序，因此每条读段只处理比对质量最优的结果，计算各个位点在当前读段对应的变异分数，并得到当前样本的变异分数函数。由于比对器会将测序错误当作变异的异常点进行比对并输出结果，进而产生“伪变异”来干扰结构变异的确定，本发明针对此问题拟合出一个函数用来表达每个位点与该位点周围的变异程度之间的关系，即变异分数。由于结构变异是由大量连续的密集单点变异组成并且数量较少，而“伪变异”是广泛随机地存在于所有读段中并且长度非常短，针对此特点，变异分数可以大致区分出结构变异与“伪变异”。由变异分数进而得到该样本的变异分数函数，解决由于测序错误对于结构变异以及复杂indel的干扰而影响其检测的问题。

进一步的，计算出各个位点在当前读段对应的变异分数，并将其预先保存在每个位点的变异分数集合中，实现了拟合出一个函数来量化每个位点与该位点周围变异程度之间的关系，从而反映“伪变异”与结构变异的区别，进而将两者区分开来。

进一步的，统计每个位点的变异分数集合中的平均数当作该位点最终的变异分数并得到此样本的变异分数函数，由此得到的变异分数函数是一个离散函数，目的是将其看作是原始序列的变异分数与“伪变异”产生的变异分数的一种叠加，使得该位点的变异程度不只由该点与参考序列的异同所决定，而是综合考虑了周围各点的变异，更适合本专利检测复杂indel的目的。

进一步的，对变异分数函数进行卡尔曼滤波降噪处理，当测序错误率加大时，只依靠变异分数函数不再能准确描述真正的变异情况，本发明使用滤波算法对变异分数函数进行滤波处理以降低“伪变异”对结构变异的影响。

进一步的，基于数据特的和使用者具备的不同的计算能力，卡尔曼滤波可以由高斯滤波代替；从数学原理层面，一方面变异位点的变异分数特征与高斯滤波原理相符，另一方面根据高斯噪声的定义——概率密度函数服从高斯分布(即正态分布)的一类噪声，也符合“伪变异”的分布特征；因此，当使用者能够确定测序错误的概率分布，且拥有较好的计算资源时，可以根据测序错误的概率分布设置参数，采用高斯滤波。

进一步的，分离出结构变异区域，提取特征，通过设定阈值分离出结构变异区域，并且提取结构变异区域的插入比率、删除比率、替换比率、变异分数最大值以及四等分位点对应的变异分数作为特征值。

进一步的，训练支持向量机(Support Vector Machine，SVM)模型，对结构变异区域分类并最终得到复杂indel的结果集。针对测序错误对复杂indel和其他结构变异的区分产生的干扰，本发明考虑到在结构变异中区分出复杂indel和其他结构变异本质上是一个二分类问题，在有测序错误的情况下，两者的特征界限变得模糊，此时只是单纯通过阈值进行线性划分已经很难将复杂indel与其他结构变异划分开来，因此借助机器学习的SVM分类模型来解决。

综上所述，由于设计了变异分数函数计算策略，并对变异分数函数进行卡尔曼滤波或高斯滤波将“伪变异”平滑掉，从而使得本发明可以降低测序错误在微观上对检测的影响；按照阈值提取出结构变异区域，初步完成了结构变异与“伪变异”的区分，进而根据在结构变异中区分复杂indel和其他结构变异是二元分类这一问题本质，从而对结构变异区域进行特征提取并训练出SVM模型进行复杂indel和其他结构变异的分类，很大程度上解决了在测序错误的情况下复杂indel和其他结构变异的特征界限变得模糊这一问题，使得本发明在不同覆盖度和测序错误率下效果良好，并且在较高测序错误率的情况下仍能保持良好的召回率和精确率，这些优势都是现有相关方法所无法达到的。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明整体流程图；

图2为变异分数的计算和卡尔曼滤波的流程图；

图3为复杂indel的形成过程图；

图4为复杂indel类型示意图，其中，(a)为3’F型，(b)为5’R型，(c)为3’F&RefR型，(d)为3’F&5’R型；

图5为测序错误对比对结果的影响示意图；

图6为变异分数计算示意图；

图7为变异分数计算的时序性示意图；

图8为结构变异区域内变异分数的特征提取图。

具体实施方式

本发明提供了一种基于滤波策略的高容错基因组复杂结构变异检测方法，基于滤波策略的高容错检测算法(CIDDⅡ)，是基于卡尔曼滤波与支持向量机(SVM)的容错机制提出的复杂indel检测算法。

请参阅图1，本发明一种基于滤波策略的高容错基因组复杂结构变异检测方法，包括以下步骤：

S1、对SAM文件预处理；

遍历SAM文件中所有的比对结果，SAM文件已经给出了读段的比对质量并以此排序，每条读段只处理比对质量最优的结果，即遍历最优质量比对读段中的CIGAR字段；

S2、根据比对后的CIGAR字段和变异分数计算准则，计算出各个位点在当前读段对应的变异分数，并将其预先保存在每个位点的变异分数集合中；

涉及到变异分数的计算，CIDDⅡ根据真正的结构变异区域的变异密度一定比含有测序错误的正常区域要高这一特点，定义了变异分数，如公式(1)所示。

具体的，S_i为位点i的变异分数，因为每个位点会被多条读段所覆盖，从而计算过程中会多次计算位点i的变异分数，分别记作{S_i1,S_i2,...,S_ij}，其中，j表示位点i的测序覆盖度，k为窗口半径，位点i变异分数将依赖于(i-k,i+k)范围的变异程度；C_j(i)为逻辑函数，表示在第j个读段中位点i在某个具体的比对结果中的变异情况，如果位点i发生了变异，则C_j(i)的值为1，反之为0。

S3、统计每个位点的变异分数集合中的平均数当作该位点最终的变异分数并得到此样本的变异分数函数；

变异分数函数即是计算出的所有位点的变异分数，其显然是个离散函数，描述了染色体位点对于参考序列的变异程度，值域为[0,1]。

S4、对变异分数函数进行卡尔曼滤波，以得到滤波降噪后的变异分数函数；

请参阅图2，卡尔曼滤波的过程噪声可以预先通过对变异分数进行随机抽样统计均值和方差得出，具体的，变异分数函数卡尔曼滤波过程如下：

S4011、通过位点i-1的最优预估值S′_i-1确定位点i的变异分数预测值

如公式(2)所示：

其中，A、B是系统参数，U是系统控制量。

S4012、根据上一次计算的误差方差P_i-1和预测过程噪声ω，预测当前位点i的误差方差

如公式(3)所示：

S4013、结合测量系统参数H和测量噪声ε，计算当前位点i的卡尔曼增益K_i，如公式(4)所示：

S4014、结合当前位点i的测量值y_i计算当前位点i的最优预估值S_i'，如公式(5)所示：

S4015、根据公式(6)更新误差方差，并重复迭代以上步骤。

其中，P_i为误差方差，H为测量系统参数，I为单位矩阵，ε为测量噪声。

卡尔曼滤波可以用高斯滤波代替，具体的，变异分数函数高斯滤波过程如下：

S4021、取变异位点i，变异位点i坐标代表其本身变异与否，若变异取1，否则取0；取距离i点最近的k个坐标点，即与位点i相邻的k个位点，根据数据量和计算能力确定k的取值，对于普通计算机建议取值为k＝8，具体地，若变异取值为1，否则为0；

S4022、设定μ和σ，由公式(7)计算出这9个位点对应的高斯模板取值G_σj，j＝1,2,...,9，即位点i的邻域点的变异与否对位点i是否为“伪变异”的影响权重；

其中，G_σ是标准差为σ的高斯核。

S4023、根据高斯模板的特性，为了确保这9个位点的高斯模板取值相加为1，对高斯模板归一化，即

得到最终的高斯模板；

S4024、利用公式(8)计算9个位点的高斯滤波I_σj，并相加得到位点i的高斯滤波的值。

I_σ＝I*G_σ (8)

其中，*是卷积操作。

经过以上滤波步骤，如果某个位点的变异是一个“伪变异”，依照测序错误的特点，它周围的变异点密度小且数量少，那么这个点及其周围的分数将会通过滤波尽量平滑掉。而如果某个位点是结构变异中的位点，它周围的变异点密度会增大，导致周围附近位点的变异分数普遍升高，滤波过程会使得增益提高，变异区域的分数会缓慢提升，从而使得结构变异与“伪变异”区分开。

S5、依照滤波后的变异分数函数，设定阈值并分离出其中的结构变异区域，提取特征；

设定阈值并分离结构变异区域具体如下：

将滤波后变异分数的均值和方差分别记作μ'和σ'，以μ'+3σ'为阈值，将变异分数函数和滤波后的变异分数函数进行上下分割并取出大于μ'+3σ'的区间，分别记作：

SV＝{(a₁,b₁),(a₂,b₂),...,(a_n,b_n)}

SV'＝{(A₁,B₁),(A₂,B₂),...,(A_n,B_n)}

将左端点A当作对SV集合的筛选的标准：对于(a,b)∈SV，如果

且A∈(a,b)，保留，否则丢弃，筛选后的SV集合中的每个元素就代表着一个结构变异区域。

提取特征具体如下：

CIDDⅡ选取结构变异区域的插入比率、删除比率、替换比率、变异分数最大值以及四等分位点对应的变异分数作为特征值，以上七个特征值分别记作insRatio、delRatio、subRatio、score₀、score₁、score₂、score₃，下面分别给出区域(a,b)内的计算公式：

其中，T_ij、P_ij、Q_ij分别表示第j个读段内目标序列，比对信息和读段序列的第i个字符的值，c为insert-size正常区域半径，I为指示函数，S_i是位点i的变异分数。

S6、训练SVM模型，再用训练好的SVM模型对结构变异区域分类并最终得到复杂indel结果集。

训练SVM模型具体如下：

S601、将结构变异集合中每个变异看作高维空间中的一个点，其中结构变异的每个特征代表一个维度；

S602、用高斯核函数将结构变异集合中的元素映射到结构变异特征个数组成的高维空间中，然后在高维空间中训练得到复杂indel和其他结构变异的最优分离超平面；

S603、最后把得到的模型用于结构变异的复杂indel分类中。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Pindel-C、SV-Bay和INDELseek是现有的拥有复杂indel检测能力的工具的代表，但是它们均只能处理第二代测序数据，而第二代测序数据的数据读段短、比对信息少等特点极大限制了可检测的复杂indel的范围，使得现有的复杂indel检测算法普遍只对小于80bps的复杂indel有敏感度，且召回率只能达到50％左右，大量的复杂indel被忽略，从而严重影响了复杂indel的下游分析。

与上述的现有复杂indel检测工具相比较，由于设计了变异分数函数计算策略，并对变异分数函数进行卡尔曼滤波或高斯滤波将“伪变异”平滑掉，从而使得本发明可以降低测序错误在微观上对检测的影响；在按照阈值提取出结构变异区域后，弄清了从结构变异中区分复杂indel和其他结构变异是二元分类这一问题本质，从而对结构变异区域进行特征提取并采用SVM模型进行复杂indel和其他结构变异的分类，很大程度上克服了在有测序错误的情况下复杂indel和其他结构变异的特征界限变得模糊这一问题，使得本发明在不同覆盖度和测序错误率下效果良好，并且在测序错误率15％的情况下仍能保持良好的召回率和精确率。

请参阅图3，展示了一种复杂indel的形成过程：正常的DNA由于某种原因产生了两个断裂点P1和P2，通常称P1和P2分别为左断点和右断点，断裂部分缺失形成了一个删除变异。紧接着DNA启动修复机制，在空间邻近的区域选取了两块与删除长度相近DNA片段a和片段b进行插入。这时，在删除变异的基础上又产生插入变异。最后，在插入过程中，片段a是反向互补后插入到对应位置，此时又可以看作片段a发生了倒置变异。由于修复时采用了双片段插入，那么在修复片段之间还存在着一个中间断点P3。图中的虚线箭头用来表示复杂indel形成过程中对DNA分子的改变。可以看出，虽然复杂indel的表现形式上是一种替换变异，但实际上复杂indel的形成过程本质上是多种简单变异(删除、插入、倒置)的互相叠加。

请参阅图4，展示了四种复杂indel的类型，为了便于表示，图中将DNA分子进行简化省略了碱基和配对。每个小图表示一种类型的复杂indel，下方黑线表示正常情况下的参考序列，上方黑线表示变异发生后的目标序列。图a表示的是一种单片段插入的复杂indel的类型，插入片段相对变异发生位置属于3’方向且为正向插入，因此该类复杂indel的类型被定义为3’F型。图b也是一种单片段插入的复杂indel，插入片段相对变异发生位置属于5’方向且为倒位后插入，因此该类复杂indel的类型被定义为5’R型。图c和图d数据多片段插入的复杂indel，它们的定义命名与单片段的类似。值得一提的是Ref表示的是插入片段的来源属于两个左右断点之间缺失序列的一部分。由于单片段的Ref型本质上属于简单的倒置变异，因此在复杂indel的分类中Ref型只会出现在多片段的复杂indel中。

请参阅图5，展示了测序错误对比对结果的影响。由于测序过程中随机出现的测序过快或者过慢导致了测序错误，其中测序过快会使连续的碱基未被检测到，测序过慢会使碱基重复测序多次。

请参阅图6，展示了单点变异分数的计算过程，0坐标对应的位点就是待计算变异分数的位点，由图可见一个位点的变异分数是周围位点共同作用的结果。

请参阅图7，展示了变异分数计算的时序性，当前T时刻计算出位点i的变异分数为S_i，代表了周围2k个位点对此位点的变异支持，而下一时刻窗口仍保存着上一窗口2k-2个位点的变异信息。

请参阅图8，展示了结构变异区域内变异分数的特征提取，选取了区域内变异分数的最大值以及四等分位点对应的变异分数。

下面结合具体的技术特征来说明本发明的优势。

1)CIDDⅡ设计的变异分数计算策略，公式1，实现了变异分数大小与该点是否变异无关且只与窗口内其他变异点有关，从而反映了“伪变异”与结构变异的区别，将进而将两者区分开来。该计算策略还实现了窗口内的其他变异位点对该点的支持度并不等价，距离越近支持度越大，距离越远支持度越小。总之，这样计算出来的变异分数不需要进行归一化更便于特征提取且位点变异分数高低和其他变异位点距离相关，更为合理。

2)CIDDⅡ采用的滤波算法可以对测序错误率加大时变异分数函数不能准确描述真正的变异情况进行降噪处理，得以将“伪变异”最大程度地平滑掉。其中，卡尔曼滤波要求序列具有时序性，而染色体的位点变异分数恰好可以看作一个时序序列；高斯滤波将中心点邻域内不同位置的位点赋予不同的权值，这恰好与周围位点对变异位点的“伪变异”与结构变异区分的贡献度大小按照距离远近而不同这一特点相符。因此，卡尔曼滤波与高斯滤波都是符合变异位点特征的滤波算法，可以有效对“伪变异”进行过滤平滑。

3)CIDDⅡ采用的SVM分类模型并且制定的特征值选取规则，充分考虑到高测序错误的影响下再简单地通过设定阈值来判断结构变异是否为复杂indel可行性已不高这一问题，将该判断问题正确划分为二分类问题，并且结合结构变异的特征来制定规则进行特征提取，进而将复杂indel和其他结构变异的区分精度大大提高。值得一说的是，在特征值选取时，考虑到结构变异区域内的各个变异比率并不能完全体现结构变异区域内的具体变异密度等信息，因此不仅选取了结构变异区域最直观的三个宏观特征(插入比率、删除比率、替换比率)作为特征值，还选取了区域内变异分数的最大值和四等分位点对应的变异分数来当作结构变异的具体特征。这样使得后期的SVM分类模型训练更加精确。

实验结果表明，在不损失召回率和精确率的情况下，CIDDⅡ可以在存在测序错误的情形下完成复杂indel的检测，并且召回率和精确率仍然能保持在70％和85％。这一结果大大优于其他复杂indel识别算法。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (4)

1.一种基于滤波策略的高容错基因组复杂结构变异检测方法，其特征在于，包括以下步骤：

S1、对SAM格式的输入文件进行预处理，遍历最优质量比对读段中的CIGAR字段；

S2、根据比对后的CIGAR字段和变异分数计算准则，计算出各个位点在当前读段对应的变异分数，并将其预先保存在每个位点的变异分数集合中，位点i的变异分数S_i为：

其中，多次计算位点i的变异分数并分别记作{S_i1,S_i2,...,S_ij}，j表示覆盖位点i的一条条读段，k为窗口半径，位点i变异分数依赖于(i-k,i+k)范围的变异程度；C_j(i)为逻辑函数，表示在第j个读段中位点i在某个具体的比对结果中的变异情况；

S3、统计每个位点的变异分数集合中的平均数当作该位点最终的变异分数并得到此样本的变异分数函数；

S4、对变异分数函数进行卡尔曼或高斯滤波，得到滤波降噪后的变异分数函数，变异分数函数卡尔曼滤波过程如下：

S4011、通过位点i-1的最优预估值S_i'_-1确定位点i的变异分数预测值

变异分数预测值

如下：

误差方差

卡尔曼增益K_i：

最优预估值S_i'：

其中，A，B是系统参数，U是系统控制量，P_i为误差方差，H为测量系统参数，I为单位矩阵，ε为测量噪声；

S4012、根据上一次计算的误差方差P_i-1和预测过程噪声ω，预测当前位点i的误差方差

S4013、结合测量系统参数H和测量噪声ε，计算当前位点i的卡尔曼增益K_i；

S4014、结合当前位点i的测量值y_i计算当前位点i的最优预估值S_i'；

S4015、更新误差方差，并重复迭代以上步骤；

变异分数函数高斯滤波过程如下：

S4021、取变异位点i，变异位点i坐标代表其本身变异与否，若变异取1，否则取0；取距离i点最近的k个坐标点，即与位点i相邻的k个位点，根据数据量和计算能力确定k的取值，对于普通计算机建议取值为k＝8；

S4022、设定μ和σ，计算出9个位点对应的高斯模板取值G_σj，j＝1,2,...,9，即位点i的邻域点的变异与否对位点i是否为伪变异的影响权重；

其中，G_σ是标准差为σ的高斯核；

S4023、对高斯模板归一化，即

得到最终的高斯模板；

S4024、计算9个位点的高斯滤波I_σj，并相加得到位点i的高斯滤波的值：

I_σ＝I*G_σ

其中，*是卷积操作；

S5、依照滤波后的变异分数函数，设定阈值并分离出结构变异区域，提取特征，设定阈值并分离结构变异区域具体如下：

将滤波后变异分数的均值和方差分别记作μ'和σ'，以μ'+3σ'为阈值，将变异分数函数和滤波后的变异分数函数进行上下分割并取出大于μ'+3σ'的区间，分别记作：

SV＝{(a₁,b₁),(a₂,b₂),...,(a_n,b_n)}

SV'＝{(A₁,B₁),(A₂,B₂),...,(A_n,B_n)}

将左端点A当作对SV集合的筛选的标准，对于(a,b)∈SV，如果

且A∈(a,b)，保留，反之则丢弃，筛选后SV集合中的每个元素代表一个结构变异区域；

S6、训练支持向量机SVM模型，再用训练好的SVM模型对结构变异区域分类并得到复杂indel结果集，训练SVM模型如下：

S601、将结构变异集合中每个变异看作高维空间中的一个点，结构变异的每个特征代表一个维度；

S602、用高斯核函数将结构变异集合中的元素映射到结构变异特征个数组成的高维空间中，然后在高维空间中训练得到复杂indel和其他结构变异的最优分离超平面；

S603、把得到的模型用于结构变异的复杂indel分类中。

2.根据权利要求1所述的基于滤波策略的高容错基因组复杂结构变异检测方法，其特征在于，步骤S3中，变异分数函数即计算出的所有位点的变异分数为离散函数，能够描述染色体位点对于参考序列的变异程度，值域为[0,1]。

3.根据权利要求1所述的基于滤波策略的高容错基因组复杂结构变异检测方法，其特征在于，步骤S4中，如果某个位点的变异为伪变异，伪变异及其周围的分数将通过滤波平滑掉；如果某个位点是结构变异中的位点，结构变异与伪变异区分开。

4.根据权利要求1所述的基于滤波策略的高容错基因组复杂结构变异检测方法，其特征在于，步骤S5中，提取特征具体如下：

选取结构变异区域的插入比率、删除比率、替换比率、变异分数最大值以及四等分位点对应的变异分数作为特征值，以上七个特征值分别记作insRatio、delRatio、subRatio、score₀、score₁、score₂、score₃，分别给出区域(a,b)内的计算公式：

其中，T_ij、P_ij、Q_ij分别是第j个读段内目标序列，比对信息和读段序列的第i个字符的值，c为insert-size正常区域半径，I为指示函数，S_i是位点i的变异分数。

Images