CN108875311B - 基于高通量测序和高斯混合模型的拷贝数变异检测方法 - Google Patents

Abstract

本发明涉及基因工程技术领域，具体涉及一种基于高通量测序和高斯混合模型的拷贝数变异检测方法，包括以下步骤：数据的生产与预处理，利用最大期望算法估计高斯混合模型的参数，窗口的拷贝数估计，合并拷贝数一致率高的相邻窗口，确定最终的拷贝数和过滤。本发明方法通过高通量测序和按窗口的检测方法，提高了分辨率，即可以检测到长度更短的CNV；通过GMM的方法，不引入参考样本，可以检测常见的拷贝数变异；通过GMM的概率过滤，提高了准确率；通过确定拷贝数变异区域后再最终确定每个实验样本的拷贝数，使群体样本的结果保持一致性。

Description

基于高通量测序和高斯混合模型的拷贝数变异检测方法

技术领域

本发明涉及基因工程技术领域，具体涉及一种基于高通量测序和高斯混合模型的拷贝数变异检测方法。

背景技术

拷贝数变异(Copy Number Variations，CNVs)是指与基因组参考序列相比，基因组中长度大于等于1000碱基(1kb)且以不同拷贝数存在的DNA片段，其形式包括插入、缺失、扩增，及其相互组合衍生出的复杂变异。Redon等根据拷贝数变异的遗传和组成形式将拷贝数变异分为5类：(a)缺失；(b)扩增；(c)同一位点并发的缺失与扩增；(d)多等位基因位点(multiple alleles)；(e)复杂难以描述的位点。通常，扩增比缺失更为常见，且覆盖更大的范围。具体来说，拷贝数变化可以通过破坏基因编码蛋白的活性部分、改变基因的表达、或者破坏基因组控制基因活性的调节区域等影响基因的活性。寻找拷贝数变异有助于在有遗传可能性的区域里寻找关键基因。

目前已实现检测拷贝数变异的方法主要有通过在一张芯片上用标记不同荧光素的样品(病例样品和对照样品)进行共杂交来检测样本基因组相对于对照基因组的拷贝数变异的比较基因组杂交芯片，通过芯片上荧光信号的强弱检测拷贝数变异(荧光信号强度与拷贝数成正相关)的单核苷酸多态检测芯片，通过测序读段深度、双端测序的插入片段长度或劈开的测序读段检测拷贝数变异的高通量测序等。

但是，比较基因组杂交芯片的分辨率较低，无法检测长度较短的拷贝数变异，同时由于需要对照样本，因此不适合检测常见的拷贝数变异。单核苷酸多态检测芯片依赖于单核苷酸多态性位点在基因组的分布，若位点较少，则分辨率低；若分布不均，则会在分布稀疏的基因组区域遗漏拷贝数变异，另外该技术通常需要对照样本，因此不适合检测常见的拷贝数变异。基于高通量测序的拷贝数变异检测技术在处理外显子组等目标区域捕获测序数据时，由于捕获效率波动带来的测序深度偏向性，其准确率较低，另外通常需要对照样本，因此不适合检测常见的拷贝数变异。现有技术多基于单个样本得到拷贝数变异，使得群体的拷贝数变异缺乏一致性。

发明内容

本发明的目的是解决检测拷贝数变异时分辨率低、准确率低、群体缺乏一致性、不适合检测常见的拷贝数变异的问题，提供一种基于高通量测序和高斯混合模型的拷贝数变异检测方法。

本发明是通过以下技术方案实现的：

基于高通量测序和高斯混合模型的拷贝数变异检测方法，包括以下步骤：

S1、数据的生产与预处理：对实验样本的基因组DNA进行高通量测序，比对测序读段至参考基因组，计算实验样本的平均测序深度，且将参考基因组按预设的长度划分为窗口，计算窗口的平均测序深度和窗口的归一化的平均测序深度；

S2、利用最大期望算法估计高斯混合模型的参数：高斯混合模型是指具有如下形式的概率分布模型：

其中，K是分模型的数目，k是第k个分模型的代号，θ是所有分模型的参数，θ_k是第k个分模型的参数且

μ_k是第k个分模型的期望，

是第k个分模型的方差，α_k是第k个分模型的系数，α_k≥0且

是高斯分布密度，

称为第k个分模型；

S3、窗口的拷贝数估计：将所有实验样本的某窗口的归一化的平均测序深度作为观测数据，带入最大期望算法，估计高斯混合模型的参数，利用朴素贝叶斯的方法计算每个实验样本在该窗口属于各分模型的概率；

S4、合并拷贝数一致率高的相邻窗口：若相邻窗口拷贝数一致的实验样本数超过90％，则认为这两个窗口属于同一个拷贝数变异，合并为一个窗口；循环合并窗口的过程，直到不再有相邻窗口可以合并；合并后的窗口即为最终的拷贝数变异区域；

S5、确定最终的拷贝数和过滤：对步骤S4得到的拷贝数变异区域，按步骤S3的方法得到最终的拷贝数，利用高斯混合模型的概率过滤，若某拷贝数变异区域内90％以上的实验样本属于所属分模型的概率大于90％，则保留此拷贝数变异区域，否则去除此拷贝数变异区域。

优选地，步骤S1中所述实验样本的平均测序深度＝比对上的测序读段数目*测序读段长度/参考基因组长度。

优选地，步骤S1中所述窗口的平均测序深度＝窗口内比对上的测序读段数目*测序读段长度/窗口的长度。

优选地，步骤S1中所述窗口的归一化的平均测序深度＝窗口的平均测序深度/实验样本的平均测序深度。

优选地，步骤S2中所述最大期望算法的步骤包括：

S21、对高斯混合模型的参数(μ_k,σ_k,α_k)取初始值；

S22、E步：依据当前模型参数，计算分模型k对观测数据y_j的响应度

S23、M步：计算新一轮迭代的高斯混合模型参数，包括

S24、判断是否收敛：若是，停止；若否，重复步骤S22～S24。

优选地，步骤S3中所述利用朴素贝叶斯的方法是使用如下公式：

其中，取使得P(k|y_j)最大的k作为j样本在该窗口的所属分模型。

本发明的有益效果在于：

(1)通过高通量测序和按窗口的检测方法，提高了分辨率，即可以检测到长度更短的CNV。

(2)通过GMM的方法，不引入参考样本，可以检测常见的拷贝数变异。

(3)通过GMM的概率过滤，提高了准确率。

(4)通过确定拷贝数变异区域后再最终确定每个实验样本的拷贝数，使群体样本的结果保持一致性。

附图说明

图1为本发明最大期望算法求解高斯混合模型参数的流程图；

图2为本发明实施例中188个样本的归一化的平均测序深度分布图。

具体实施方式

为更好理解本发明，下面结合实施例及附图对本发明作进一步描述，以下实施例仅是对本发明进行说明而非对其加以限定。

实施例HLA-DRB5基因拷贝数变异在188个实验样本的检测

1.取188份来自不同人的血液样品，提取DNA。使用MHC芯片捕获，然后基于Illumina Hiseq2000平台进行高通量测序。使用针对HLA-DRB3/4/5基因的特异性引物进行PCR，对PCR产物做Sanger测序。

2.对高通量测序数据，按本发明的方法，检测拷贝数变异，即包括以下步骤：

(1)数据的生产与预处理。对实验样本的基因组DNA进行高通量测序。比对测序读段至参考基因组。计算实验样本的平均测序深度(实验样本的平均测序深度＝比对上的测序读段数目*测序读段长度/参考基因组长度)。将参考基因组按预设的长度划分为窗口。计算窗口的平均测序深度(窗口的平均测序深度＝窗口内比对上的测序读段数目*测序读段长度/窗口的长度)。计算窗口的归一化的平均测序深度(窗口的归一化的平均测序深度＝窗口的平均测序深度/实验样本的平均测序深度)。

(2)高斯混合模型(GMM)与最大期望算法(EM算法)。GMM是指具有如下形式的概率分布模型：

其中，α_k是系数，α_k≥0，

是高斯分布密度，

称为第k个分模型。EM算法可以用于依据观测数据y₁,y₂,…,y_N估计GMM的参数。EM算法的步骤和公式见附图1。

(3)窗口的拷贝数估计。将所有实验样本的某窗口的归一化的平均测序深度作为观测数据，带入EM算法，可以估计GMM的参数。利用朴素贝叶斯的方法可以计算每个实验样本在该窗口属于各分模型的概率，公式如下：

取使得P(k|y_j)最大的k作为j样本在该窗口的所属分模型。在使用EM算法时，GMM的分模型数目K的取值和参数(μ_k,σ_k,α_k)的初始取值至关重要，这些值与真实值越接近，拷贝数估计越准确。我们取K＝2,3,4或5，分别带入EM算法，最后计算有效性指标S_Dbw(Halkidi,M.andM.Vazirgiannis.Clustering validity assessment:Finding the optimalpartitioning of a data set.in Data Mining,2001.ICDM 2001,Proceedings IEEEInternational Conference on.2001.IEEE.)，选取使得指标取得最大值的K。参数(μ_k,σ_k,α_k)的初始取值如下：

α_k＝1/K

GMM的每个分模型代表一个拷贝数状态，但还需要其他信息确定具体的拷贝数。我们使用分模型的参数μ_k来确定具体的拷贝数：

(4)合并拷贝数一致率高的相邻窗口。若相邻窗口拷贝数一致的实验样本数超过90％，则认为这两个窗口属于同一个拷贝数变异，合并为一个窗口。循环合并窗口的过程，直到不再有相邻窗口可以合并。合并后的窗口即为最终的拷贝数变异区域。

(5)确定最终的拷贝数和过滤。对第(4)步得到的拷贝数变异区域，按第(3)步的方法得到最终的拷贝数。利用GMM的概率过滤，若某拷贝数变异区域内90％以上的实验样本属于所属分模型的概率大于90％，则保留此拷贝数变异区域，否则去除此拷贝数变异区域。

3.最终，188个样本的归一化的平均测序深度分布如图2所示。

对比例HLA-DRB5基因拷贝数变异在188个实验样本的检测(金标准)

使用实施例步骤1中的Sanger测序数据，人工判断HLA-DRB5基因的型别和拷贝数。

比较本发明方法在HLA-DRB5基因检测到的拷贝数变异和Sanger测序得到的金标准，一致率为97.14％。

以上所述实施方式仅仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案作出的各种变形和改进，均应落入本发明的权利要求书确定的保护范围内。

Claims (3)

1.基于高通量测序和高斯混合模型的拷贝数变异检测方法，其特征在于，包括以下步骤：

S1、数据的生产与预处理：对实验样本的基因组DNA进行高通量测序，比对测序读段至参考基因组，计算实验样本的平均测序深度，且将参考基因组按预设的长度划分为窗口，计算窗口的平均测序深度和窗口的归一化的平均测序深度；

所述实验样本的平均测序深度＝比对上的测序读段数目*测序读段长度/参考基因组长度；

所述窗口的平均测序深度＝窗口内比对上的测序读段数目*测序读段长度/窗口的长度；

所述窗口的归一化的平均测序深度＝窗口的平均测序深度/实验样本的平均测序深度；

S2、利用最大期望算法估计高斯混合模型的参数：高斯混合模型是指具有如下形式的概率分布模型：

其中，K是分模型的数目，k是第k个分模型的代号，θ是所有分模型的参数，θ_k是第k个分模型的参数且

μ_k是第k个分模型的期望，

是第k个分模型的方差，α_k是第k个分模型的系数，α_k≥0且

是高斯分布密度，

称为第k个分模型；

S3、窗口的拷贝数估计：将所有实验样本的某窗口的归一化的平均测序深度作为观测数据，带入最大期望算法，估计高斯混合模型的参数，利用朴素贝叶斯的方法计算每个实验样本在该窗口属于各分模型的概率；

S4、合并拷贝数一致率高的相邻窗口：若相邻窗口拷贝数一致的实验样本数超过90％，则认为这两个窗口属于同一个拷贝数变异，合并为一个窗口；循环合并窗口的过程，直到不再有相邻窗口可以合并；合并后的窗口即为最终的拷贝数变异区域；

S5、确定最终的拷贝数和过滤：对步骤S4得到的拷贝数变异区域，按步骤S3的方法得到最终的拷贝数，利用高斯混合模型的概率过滤，若某拷贝数变异区域内90％以上的实验样本属于所属分模型的概率大于90％，则保留此拷贝数变异区域，否则去除此拷贝数变异区域。

2.根据权利要求1所述的基于高通量测序和高斯混合模型的拷贝数变异检测方法，其特征在于，步骤S2中所述最大期望算法的步骤包括：

S21、对高斯混合模型的参数(μ_k,σ_k,α_k)取初始值；

S22、E步：依据当前模型参数，计算分模型k对观测数据y_j的响应度

S23、M步：计算新一轮迭代的高斯混合模型参数，包括

S24、判断是否收敛：若是，停止；若否，重复步骤S22～S24。

3.根据权利要求1所述的基于高通量测序和高斯混合模型的拷贝数变异检测方法，其特征在于，步骤S3中所述利用朴素贝叶斯的方法是使用如下公式：

其中，取使得P(k|y_j)最大的k作为j样本在该窗口的所属分模型。

Images