基于外显子测序数据的拷贝数变异检测方法及系统、终端和存储介质与流程

本发明涉及外显子测序技术领域，尤其是涉及一种基于外显子测序数据的拷贝数变异检测方法及系统、终端和存储介质。

背景技术：

人类是二倍体，有两套遗传信息，分别来自父母，通常认为dna分子有两个拷贝。拷贝数变异(copynumbervariants,cnvs)是dna序列结构变异的一种，通常定义为dna片段(>1kb)的扩增或删除。

因为发生cnvs的区域一般比较大，大概率会包含若干基因，所以它们在人类疾病和药物反应中具有重要的作用。若cnvs包含影响关键的发育基因会引发很严重的疾病，比如自闭症、精神分裂症、先天性心脏畸形、帕金森、老年痴呆症等等。而且cnvs与癌症的发生也密切相关，导致致癌基因的激活归因于拷贝数的扩增和抑癌基因的失活归因于拷贝数的减少。所以，对cnvs的研究可以在疾病、癌症的治疗和预后发挥重要的作用，也有助于发现一些药物的靶点，应用于临床治疗。

目前检测cnvs主要有两种技术：

1)微阵列比较基因组杂交技术。这种技术有其局限性,检测的结果的灵敏度与精度相对比较低。

2)二代测序技术。随着测序技术快速发展，基于测序技术的cnvs检测结果有更高的分辨率和准确度，而且测序成本不断降低。全基因组测序技术与外显子测序技术成为二代测序技术中检测cnvs和研究人类疾病的主要策略。因为外显子只占了人类基因组的1％左右，外显子测序技术比全基因组测序技术更便宜、更快。

dna序列由a、g、c、t四种碱基组成，测序实际上就是通过仪器得到由这四种碱基组成的一条很长的字符串。外显子测序技术就是将位于外显子区域的dna捕获并富集，通过测序得到外显子区域的字符串。所以，使用外显子测序数据检测病人的cnvs成为研究人类疾病主要方式之一。

与全基因组测序技术相比，在外显子数据中检测cnvs存在更高的难度。外显子测序技术的文库准备的时候，杂交这一过程引入了偏差；一些外显子区域内的reads数比较低；全基因组测序数据中reads是连续的，而外显子的reads比对到参考基因组却是孤立的。这些因素使得在外显子测序数据中检测拷贝数变异带来了额外的挑战。

目前，国内外已有的使用外显子测序数据检测cnvs的工具比较好的有以下三个工具：

1)canoes：利用负二项分布模型拟合数据，去噪声的算法使用的是广义相加模型，去掉gc含量对数据的影响。然后使用hmm和viterbi算法将拷贝数划分成删除、正常和扩增，但只对比较小的cnvs有很高的敏感度，而且没有考虑外显子区域信号分布不均匀、以及杂交等因素导致的偏差。

2)xhmm：利用高斯近似模型拟合数据，并且用pca标准化数据，去掉前k个潜在的影响噪声的因子，然而假定所有噪声都是随机的高斯噪声，但gc含量对数据的影响并不是随机的。

3)convex：同样使用高斯近似模型来拟合数据，但去噪声的算法采用的是通过比对正常样本，设置对照组来减少一些特定的噪声，然而没有考虑杂交等因素导致的偏差和噪声。

技术实现要素：

为了解决现有技术中外显子测序数据中检测不足的问题，本发明提供了一种基于外显子测序数据的拷贝数变异检测方法及系统、终端和存储介质。采用该方法可以降低成本，实现外显子测序数据拷贝数变异检测的准确性和有效性。

为达到上述目的，本发明采用以下技术方案：

一种基于外显子测序数据的拷贝数变异检测方法，包括以下步骤：

s100，获取正常样本和测试样本，对正常样本的外显子测序数据进行数据清理，然后对数据进行标准化处理，生成标准化后的正常样本集数据矩阵；

s200，根据正常样本集数据矩阵中每个外显子区域在所有样本中的离散程度，将外显子区域划分成稳定与不稳定的区域；在稳定的区域中，处理测序技术中批次效应因素对正常样本集与测试样本比对的影响，以此构建参考数据矩阵；

s300，使用pca方法对参考数据矩阵进行处理，通过用主成分重构原始数据，将参考数据矩阵转换到其他空间并得到新的参数；

s400，将测试数据变换到参考数据矩阵使用pca转换后的空间中，然后使用z-score方法得到测试数据与参考数据矩阵在当前空间中的差异程度，将外显子区域划分成扩增、正常、删除三种状态，根据各个外显子区域的状态使用贪心算法进行合并，完成对测试样本的拷贝数变异的检测。

步骤s100具体包括：

s101，根据样本集文件自身提供的质量信息对正常样本集的外显子测序数据进行质量控制，获得高质量的测序数据，形成正常样本集外显子测序数据对应的readscount数据矩阵；

s102，根据外显子测序数据自身的特点进行进一步清理,过滤掉数据矩阵中不符合参数值要求的外显子区域；

s103，对数据进行标准化处理，将每个样本每个外显子区域的reads数标准化后的值构建形成正常样本集数据矩阵。

s102中，所述参数值要求包括设定最低的外显子长度、外显子上可接受的最低gc含量、外显子上可接受的最高gc含量、所有样本在一个外显子上可接受的最低的reads数的中位数。

s103中，每个样本每个外显子区域的reads数标准化后的值为外显子上的reads数除以对应样本测序数据中的数据清理后的reads数总和。

步骤s200具体包括：

s201，使用变异系数判断正常样本集数据中每个外显子区域在所有样本中的离散程度，根据设定的变异系数阈值将外显子区域划分成稳定与不稳定的区域；

s202，将测试样本的外显子测序数据按照正常样本集数据矩阵的构建流程生成对应的测试样本数据矩阵，然后在稳定的外显子区域中计算测试样本与正常样本集之间的距离，进而生成参考数据矩阵。

一种基于外显子测序数据的拷贝数变异检测系统，包括：

正常样本集数据矩阵生产模块，用于获取正常样本和测试样本，对正常样本的外显子测序数据进行数据清理，然后对数据进行标准化处理，生成标准化后的正常样本集数据矩阵；

参考数据矩阵构建模块，用于根据正常样本集数据矩阵中每个外显子区域在所有样本中的离散程度，将外显子区域划分成稳定与不稳定的区域；在稳定的区域中，处理测序技术中批次效应因素对正常样本集与测试样本比对的影响，以此构建参考数据矩阵；

参考数据矩阵处理模块，使用pca方法对参考数据矩阵进行处理，通过用主成分重构原始数据，将参考数据矩阵转换到其他空间并得到新的参数；

拷贝数变异的检测模块，将测试数据变换到参考数据矩阵使用pca转换后的空间中，然后使用z-score方法得到测试数据与参考数据矩阵在当前空间中的差异程度，将外显子区域划分成扩增、正常、删除三种状态，根据各个外显子区域的状态使用贪心算法进行合并，完成对测试样本的拷贝数变异的检测。

所述正常样本集数据矩阵生产模块具体用于：

根据样本集文件自身提供的质量信息对正常样本集的外显子测序数据进行质量控制，获得高质量的测序数据，形成正常样本集外显子测序数据对应的readscount数据矩阵；

根据外显子测序数据自身的特点进行进一步清理,过滤掉数据矩阵中不符合参数值要求的外显子区域；

对数据进行标准化处理，将每个样本每个外显子区域的reads数标准化后的值构建形成正常样本集数据矩阵。

所述参考数据矩阵构建模块具体用于：

使用变异系数判断正常样本集数据中每个外显子区域在所有样本中的离散程度，根据设定的变异系数阈值将外显子区域划分成稳定与不稳定的区域；

将测试样本的外显子测序数据按照正常样本集数据矩阵的构建流程生成对应的测试样本数据矩阵，然后在稳定的外显子区域中计算测试样本与正常样本集之间的距离，进而生成参考数据矩阵。

一种终端，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现所述的基于外显子测序数据的拷贝数变异检测方法。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行如权利要求1至5任一项所述的一种基于外显子测序数据的拷贝数变异检测方法。

本发明相比于现有技术，具有以下优点：

本发明以外显子测序数据为目标对象，对测试样本检测其拷贝数变异，采用该方法可以降低成本，实现外显子测序数据拷贝数变异检测的准确性和有效性。

具体优点如下：

第一：数据清理时考虑的更加周全，首先对单个样本的测序数据进行清洗，去掉低质量、重复的reads，然后通过多样本间的比对，去掉数据量非常小测序时无法捕获到序列的外显子区域，最后考虑到gc含量对测序数据的影响，若外显子区域存在极高的gc含量，其reads数的分布完全不均匀，清理掉对应的区域可以有效防止模型过于复杂导致过拟合。从源头上确保拷贝数变异检测系统的可靠性。

第二：将基因组划分成稳定与不稳定的区域，这符合种群内存在基因组多态性的特点，进而对不稳定的区域进行进一步的分析：做基因本体论功能富集分析，从中发现样本的特异性，样本特异性的区域，这些区域内样本特异的特征，以及这些区域的功能是什么。

第三：充分考虑了样本测序数据之间的差异。由于实验批次、实验条件的不同会导致样本间的测序数据有较大的差异，通过计算样本与样本间的相似系数，过滤掉与测试样本相似系数较低的样本，只将距离较近的样本作为参考，降低了去除噪声的难度。

第四：本发明可以快速检测出测试样本发生的拷贝数变异，参考样本的测序数据只需处理一次产生对应的矩阵，测试样本的测序数据只需按照对应的流程、参数转换成对应的数据即可。处理测序数据时，只关注落在外显子区域中的reads数，线性时间内就可以将测序数据转换成所需的readscount数据矩阵。

第五：本发明提供了一个可交互界面供临床研究人员使用，只需简单地几步既可以得到测试样本的拷贝数变异的区域，查找感兴趣的基因是否发生拷贝数变异。

附图说明

图1为基于外显子测序数据的拷贝数变异检测系统的流程图；

图2为参考样本使用pca前数据矩阵的热图；

图3为参考样本使用pca后数据矩阵的热图；

图4为基于外显子测序数据的拷贝数变异检测系统的参数设置界面；

图5为基于外显子测序数据的拷贝数变异检测系统的导入输入文件界面；

图6为基于外显子测序数据的拷贝数变异检测系统的拷贝数变异的染色体选择界面；

图7为基于外显子测序数据的拷贝数变异检测系统的拷贝数变异结果搜索展示界面。

具体实施方式

本发明一种基于外显子测序数据的拷贝数变异检测方法，包括以下步骤：

步骤1：获取正常样本和测试样本，对正常样本的外显子测序数据进行数据清理，然后对数据进行标准化处理，生成标准化后的正常样本集数据矩阵。

第一步，根据文件自身提供的质量信息对正常样本集的外显子测序数据进行质量控制，去掉重复、低质量的reads(高通量测序仪测序所得到的碱基序列)，对数据进行初步清理，获得高质量的测序数据，然后得到正常样本集外显子测序数据对应的readscount数据矩阵。

第二步，根据外显子测序数据自身的特点进行进一步清理。设定最低的外显子长度、外显子上可接受的最低gc含量、外显子上可接受的最高gc含量、所有样本在一个外显子上可接受的最低的reads数的中位数这四个参数值，过滤掉数据矩阵中不符合上述参数值的外显子区域，保证后续的分析是在可靠的数据上进行的。

第三步，对数据进行标准化处理。在第一步、第二步对数据的清理的基础上进行标准化处理。每个样本每个外显子区域的reads数标准化后的值为外显子上的reads数除以对应样本测序数据中的数据清理后的reads数总和，每个样本每个外显子区域的reads数标准化后的值构建形成正常样本集数据矩阵。

步骤2：根据正常样本集数据矩阵中每个外显子区域在所有样本中的离散程度，进而将外显子区域划分成稳定与不稳定的区域。在不稳定的区域，可以用来做基因本体论功能富集分析；在稳定的区域中，处理测序技术中批次效应这一因素对正常样本集与测试样本比对的影响，构建参考数据矩阵。

第一步，使用变异系数判断正常样本集数据中每个外显子区域在所有样本中的离散程度，根据设定的变异系数阈值将外显子区域划分成稳定与不稳定的区域。

第二步，将测试样本的外显子测序数据按照正常样本集数据矩阵的构建流程生成对应的测试样本数据矩阵，然后在稳定的外显子区域中计算测试样本与正常样本集之间的距离，进而生成参考数据矩阵。

步骤3：使用pca(主成分分析方法)对参考数据矩阵进行处理，用主成分重构原始数据，将参考数据矩阵转换到其他空间。

步骤4：利用步骤3中使用pca转换到其他空间训练出来的参数，将测试数据变换到参考数据矩阵使用pca转换后的空间中。然后使用z-score方法得到测试数据与参考数据矩阵在当前空间中的差异程度。进而将外显子区域划分成扩增、正常、删除三种状态，根据各个外显子区域的状态使用贪心算法进行合并，完成对测试样本的拷贝数变异的检测。

本发明所述的基于外显子测序数据的拷贝数变异检测系统，主要包括五个部分，基于外显子测序数据的拷贝数变异检测方法四个部分，及利用可交互界面，模拟使用基于外显子测序数据的拷贝数变异检测系统一个部分。各部分具体如下。

正常样本集数据矩阵生产模块，用于获取正常样本和测试样本，对正常样本的外显子测序数据进行数据清理，然后对数据进行标准化处理，生成标准化后的正常样本集数据矩阵；

参考数据矩阵构建模块，用于根据正常样本集数据矩阵中每个外显子区域在所有样本中的离散程度，将外显子区域划分成稳定与不稳定的区域。在稳定的区域中，处理测序技术中批次效应因素对正常样本集与测试样本比对的影响，以此构建参考数据矩阵；

参考数据矩阵处理模块，使用pca方法对参考数据矩阵进行处理，通过用主成分重构原始数据，将参考数据矩阵转换到其他空间并得到新的参数；

拷贝数变异的检测模块，将测试数据变换到参考数据矩阵使用pca转换后的空间中，然后使用z-score方法得到测试数据与参考数据矩阵在当前空间中的差异程度，将外显子区域划分成扩增、正常、删除三种状态，根据各个外显子区域的状态使用贪心算法进行合并，完成对测试样本的拷贝数变异的检测。

所述正常样本集数据矩阵生产模块具体用于：

根据样本集文件自身提供的质量信息对正常样本集的外显子测序数据进行质量控制，获得高质量的测序数据，形成正常样本集外显子测序数据对应的readscount数据矩阵；

根据外显子测序数据自身的特点进行进一步清理,过滤掉数据矩阵中不符合参数值要求的外显子区域；

对数据进行标准化处理，将每个样本每个外显子区域的reads数标准化后的值构建形成正常样本集数据矩阵。

所述参考数据矩阵构建模块具体用于：

使用变异系数判断正常样本集数据中每个外显子区域在所有样本中的离散程度，根据设定的变异系数阈值将外显子区域划分成稳定与不稳定的区域；

将测试样本的外显子测序数据按照正常样本集数据矩阵的构建流程生成对应的测试样本数据矩阵，然后在稳定的外显子区域中计算测试样本与正常样本集之间的距离，进而生成参考数据矩阵。

1、正常样本集测序数据的数据预处理

对测序数据进行数据清洗，原始测序数据中存在次优的比对结果的reads、没有通过质量控制的reads、由于pcr或测序错误产生的reads、补充匹配的reads等低质量的数据，这些存在异常的reads数据甚至可能影响到拷贝数变异检测的结果，去掉这些reads以及mappingquality低于设定值的reads，从而获得高质量可靠的数据。

区域过小的外显子上的reads分布完全不均匀，根据gc含量对reads分布的影响发现，具有过低以或过高的gc含量的外显子上的reads分布差异很大，以及一些外显子上的reads数比较低，这些外显子上的数据是不可靠的。由于样本间的测序深度不一致，必须对数据进行标准化处理，将数据放在同一个量纲上进行分析。

首先通过所有样本间的比对，去掉那些测序时无法捕获到reads的外显子区域，其次去掉外显子长度很低的区域，然后清理掉gc含量极高或极低的外显子区域，进而有效防止模型过于复杂导致过拟合，从源头上确保拷贝数变异检测系统的可靠性。最后对数据进行标准化处理，计算所有样本中可靠的外显子区域的reads数之和，数据矩阵标准化后的值为原始的reads数除以对应自身样本的reads数总和。

从而完成对原始数据的预处理，得到正常样本外显子测序数据对应的readscount数据矩阵。

2、基因组稳定与不稳定区域的划分以及构建参考数据矩阵

种群内存在多样性，其内在是基因组的多态性。种群内基因组的一些区域因为外在、内在的因素可能存在多个拷贝数的状态，即拷贝数在这些区域内不稳定。因此将基因组划分成稳定与不稳定的区域，对不稳定的区域做基因本体论功能富集分析，从中发现样本的特异性，样本特异性的区域，这些区域内样本特异的特征，以及这些区域的功能是什么。在样本稳定的区域中分析拷贝数的状态，扩增、正常、删除。使用变异系数划分稳定与不稳定的区域，其公式如下：

其中，std和mean表示当前外显子在所有样本中的标准化后的reads数那一列的标准差、均值，cv表示当前外显子的变异系数。

将测试样本测序数据按照正常样本集数据矩阵同样的流程、参数转换成对应的测试样本数据矩阵，由于实验批次、实验条件等批次效应的影响，样本与样本间的测序数据有较大的差异。可以使用jensen-shannon距离、相关系数、皮尔逊相关系数等计算测试样本与正常样本间的距离，这里采用的是jensen-shannon距离，进而得到正常样本集与测试样本测序数据间的距离矩阵，过滤掉与测序样本测序数据距离较远的样本，进而生成了参考数据矩阵。降低了去除噪声的难度，也防止结果出现偏差。

3、pca变换处理参考数据矩阵

使用pca对参考数据矩阵进行处理。虽然readscount与拷贝数变异在一个区域内是相关的，但gc含量、一些外显子区域的reads分布不均匀、杂交导致的偏差、样本间的差异等因素，会导致偏差和噪声。将这些因素看出影响reads数的次要因素，任何一个主要成分对测序数据reads数的影响远远大于这些噪声的影响，所以相对于这些次要因素，主要成分相对不受影响，即可用主要成分重构原始数据的主要信息，使用pca将数据矩阵转换到其他空间。附图2为使用pca前数据矩阵的热图，每行代表一个样本，每列代表一个外显子区域，每列的颜色越相近说明对应的样本集测序数据在该列对应的外显子区域内差异较小，使用pca的目标是使每列的差异比较小。附图3位使用pca后数据矩阵的热图。通过两图的对比，说明使用pca对数据进行处理效果良好。

4、拷贝数变异检测

测试样本的数据矩阵使用上一步参考样本数据使用pca转换到其他空间训练出来的参数，转换到参考样本变换后的空间。使用z-score判断测试样本的外显子区域是否发生拷贝数扩增或删除变异。通常一个cnv区域大于1kb，根据各个外显子区域的状态、外显子的相邻情况，使用贪心算法将异常外显子区域的合并，确定拷贝数变异的状态及其长度范围，完成对测试样本的拷贝数变异的检测。

5、利用可交互界面对拷贝数变异检测系统的模拟使用

本发明基于使用外显子测序数据的拷贝数变异检测系统，提供一种基于pyqt5实现的可交互界面，作为一个拷贝数变异检测系统供用户使用。

第一，附图4为参数设置模块。mappingquality设为20，最低外显子长度设为20，gc含量下限(％)设为10，gc含量上限设为90，最低readscount设为20，z-score阈值设为3，js散度阈值设为0.9，合并外显子设为3，进程数设为2。点击下一步或选择文件。

第二，附图5为选择文件模块。导入正常样本集外显子测序数据的存储路径；导入测试样本的bam文件；输入外显子的区域或对应的外显子bed文件；导入上一步中外显子的区域中每个外显子对应的gc含量的文件。点击运行，等待对正常样本集测序数据的处理以及对测试样本的拷贝数变异检测。

第三，附图6为对测试样本的染色体进行选择查看。这里点击8号染色体，出现一个默认的区域结果展示图，如附图7所示。可以通过在搜索框中输入基因组区域多基因名，或在染色体上用鼠标选取区域查看感兴趣的区域是否发生变异。

输入的区域为chr8:66454059-66468352，在8号染色体中的位置在附图7中染色体单体示意图中红色的线所在的位置，包含adhfe1基因，一共包含5个外显子，这5个外显子都发生了拷贝数扩增的变异。

本发明还提供一种终端，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如所述的基于外显子测序数据的拷贝数变异检测方法。

具体地显示系统如图4至图7所示。

本发明还提供一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行如所述的一种基于外显子测序数据的拷贝数变异检测方法。

以上内容是对本发明所作的进一步详细说明，不能认定本发明的具体实施方式仅限于此，对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单的推演或替换，都应当视为属于本发明由所提交的权利要求书确定保护范围。