一种拷贝数变异检测方法、系统、存储介质、计算机设备与流程

本发明属于基因测序技术领域，尤其涉及一种拷贝数变异检测方法、系统、存储介质、计算机设备。

背景技术：

目前：密度峰值算法是一种基于密度的算法，能够快速找到全局的异常点。整个算法最主要的是以下三个概念：dij：表示第i个点和第j个点之间的欧式距离。ρi：与第i个点的距离小于dc的点的个数，其中dc为一个阈值。σi：比第i个点密度大的点集合中，与第i个点的最小的欧式距离，若第i个点的密度最大，则σi等于与第i个点最大的欧式距离。密度峰值算法的主要思想：首先计算每个数据点的密度以及它们之间欧式距离，根据该距离矩阵找出每个数据点(以点a为例)周围的密度比它大且欧式距离最近的点，假设该点为点b，计算两点之间的欧式距离，作为点a的距离σ，密度最大的点的距离σ由离点a最远的点与点a的欧氏距离计算得到。这样得到所有点的密度和距离σ，根据这两个特征作高斯分布统计，异常点的特点是密度小距离大，由此特点并依据每个点的概率密度判断该点是否为异常点。

拷贝数变异：是dna结构变异的一种情况，包括长度为1kb以上的dna片段的扩增或丢失，主要表现为亚显微水平的缺失和重复。cnv是基因组结构变异(structuralvariation,sv)的重要组成部分。cnv位点的突变率远高于snp(singlenucleotidepolymorphism),是人类疾病的重要致病因素之一。cnv的形成机制有多种，可分为dna重组和dna错误复制两大类。cnv可引起单基因病与罕见疾病,同时与复杂疾病也相关。其致病的可能机制有基因剂量效应、基因断裂、基因融合和位置效应等。随着对cnv的深入研究,人们对人类基因组的构成、个体间的遗传差异以及遗传致病因素也有了新的认识。一般情况下，扩增较丢失更为多见、覆盖的范围也相对更大，这主要是因为染色体大片段丢失往往会造成比扩增更为严重的表型，使得后代更难在进化中生存下来。借助新一代测序技术和相应的实验策略可以对cnv进行发掘和精确定位。

目前，针对新一代测序技术的单样本拷贝数变异检测方法主要有以下几种：readdepth(基于读段深度的检测方法，简称rd)；通过窗口的读段深度来检测拷贝数变异；pair-end(双端比对的检测方法)；根据比较pair-end两端之间距离(插入片段的长度)与参考序列上的差异来确认拷贝数变异；splitread：将序列进行分割，比较该片段与参考序列的差异来判断拷贝数变异；

目前大多数使用readdepth方法来检测拷贝数变异，而pair-end、splitread方法大多用来检测其他结构体变异。从测序数据中获得的读段深度与拷贝数变异相关的，发生拷贝数扩增区域的读段深度高，而拷贝数丢失的区域读段深度较低，所以可以根据每个窗口的读段深度来检测拷贝数变异。现阶段使用该方法来检测拷贝数变异的主要思想可分为以下几种：提取数据特征并建立概率统计模型，通过计算每个区域变异的概率来确定拷贝数变异。使用某些算法来检测异常区域，进而判断该异常区域的拷贝数变异类型，如随机森林算法，凸包算法等，局部异常点算法，孤立森林算法。通过检测序列的snp来计算出拷贝数变异,通过对信号数据进行统计建模，算法可以正确检测出染色体基因型，并估计样本中正常细胞的比例。但当样本中正常细胞比例过高时，由于肿瘤细胞信号过小，伴随噪声的影响，这些方法检测结果会不准确。

现有技术虽然在拷贝数变异检测中有较好的性能，但总体存在以下不足之处：检测方法过于依赖原数据，导致低覆盖度数据的检测结果不准确；检测方法使用的阈值由人为确定，导致检测结果的灵敏度不高；检测方法涉及到迭代或递归的过程，导致计算复杂度较大。

通过上述分析，现有技术存在的问题及缺陷为：现有技术对拷贝数变异检测灵敏度低，拷贝数变异检测不准确，计算复杂度较大。

解决以上问题及缺陷的难度为：找到一种可以对数据进行分割的算法，不仅可以合并差异不大的相邻数据点，而且不会平滑掉数据的起伏特征，即不会将变异数据的特征弱化；找出数据之间存在的联系，并提取可以区分数据是否变异的有效特征；找一种简洁、合适的算法利用所提取的特征对数据进行分析，对这些数据作统计分析，通过简单的计算量化数据的变异程度。

解决以上问题及缺陷的意义为：因为拷贝数变异是以bin为单位进行检测的，相邻bin之间又存在内在联系，所以找到一种合适的分割方法对数据检测是很有必要的，之后我们的检测以分割之后的每一个片段为单位进行检测，因为每一个片段都是对多个bin的平均，所以提高了检测的准确率；检测拷贝数变异我们只提取了序列的readdepth作为特征，其次就是bin的位置特征，但是这两个特征没有联系，所以我们根据这两个特征得到每个片段的密度和距离，利用这两个特征对所有的数据执行densitypeak算法，将正常区域和变异区域分开；由于染色体序列数据庞大，如果我们的算法计算复杂度还高，无疑会损坏硬件设备，并且非常耗时，该算法大大减少了运行时间，我们再对每个片段作多元高斯分布统计，得到每个片段变异的概率密度，概率密度就表示该片段的变异程度，然后选取一个合适的概率密度阈值，快速得到最终结果。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种拷贝数变异检测方法、系统、存储介质、计算机设备。

本发明是这样实现的，一种拷贝数变异检测方法，所述拷贝数变异检测方法包括：

第一步，从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

第二步，对经过处理后的数据进行segment计算；

第三步，使用密度峰值检测算法对segment计算产生的数据进行计算；

第四步，得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

第五步，得到变异区域，确定变异是丢失还是增添。

进一步，所述第二步segment计算是采用cghflasso分割算法将数据平滑处理，将散乱的数据分割成长度不等的片段，为seg；将segment后的一维readdepth数据加上位置信息转化成二维数据。

进一步，所述第三步还包括：

(1)计算每个seg的与其他seg的距离，得到距离矩阵dis_matrix；

(2)通过dis_matrix和dc参数求出每个seg的密度density；

(3)通过dis_matrix和density求出簇中心的seg集合和每个seg对应的距离distance；若sega为全局密度最大的seg，则a的距离为a与离a最远的seg之间的的距离，否则，a的距离为a与密度比a大且离a最近的seg之间的距离；

(4)对每个seg执行(1)～(3)，得到所有seg的密度和距离特征。

进一步，所述第四步对二维数据进行多元正态分布；根据多元正态分布的概率密度函数求出每个seg的概率密度；设定一个可调整的概率密度阈值prob，根据每个seg的概率密度判断该seg是否发生变异，若某seg的概率密度小于prob并且其密度小于全局平均密度，认为该seg是一个异常片段；通过每个seg的密度和距离两个特征就筛选出变异区域。

进一步，所述第五步确定变异类型，簇中心的区域为正常区域，取簇中心的readdepth值的平均值作为标准值；若该变异区域的readdepth值大于标准值，则认为该区域为拷贝数扩增区域duplication，若该变异区域的readdepth值小于标准值，则认为该区域为拷贝数丢失区域deletion。

进一步，所述拷贝数变异检测方法还包括：通过位置特征和readdepth特征得到每个seg的密度和距离特征，通过密度和距离进行操作；使用余弦相似性特征代替位置特征来计算每个seg的密度和距离，余弦相似性是通过位置和readdepth二维数据通过计算得到的；有从左到右y，x，z三个点，当x是异常点时，向量y-x和z-x之间的夹角不会因为不同的y和z值而剧烈变化，用余弦相似性代替位置特征。

本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如下步骤：

第一步，从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

第二步，对经过处理后的数据进行segment计算；

第三步，使用密度峰值检测算法对segment计算产生的数据进行计算；

第四步，得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

第五步，得到变异区域，确定变异是丢失还是增添。

本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：

第一步，从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

第二步，对经过处理后的数据进行segment计算；

第三步，使用密度峰值检测算法对segment计算产生的数据进行计算；

第四步，得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

第五步，得到变异区域，确定变异是丢失还是增添。

本发明的另一目的在于提供一种运行所述拷贝数变异检测方法的拷贝数变异检测系统，所述拷贝数变异检测系统包括：

数据预处理模块，用于从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

segment计算模块，用于对经过处理后的数据进行segment计算；

数据计算模块，用于使用密度峰值检测算法对segment计算产生的数据进行计算；

筛选模块，用于得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

变异类型确定模块，用于得到变异区域，确定变异是丢失还是增添。

本发明的另一目的在于提供一种拷贝数变异检测终端，所述拷贝数变异检测终端搭载所述的拷贝数变异检测系统。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：

本发明的主要创新点：

1.充分考虑了拷贝数相邻bin之间存在的联系，提出一种数据分割方法，将我们的分析单位从bin转化为seg。

2.根据拷贝数变异的特性(变异点被孤立，而正常点聚集)，运用密度峰值算法提取出数据的密度和距离这两个与拷贝数密切相关的特征，然后进一步分析。

本发明还可以利用数据的余弦相似性特征代替位置特征，然后结合readdepth特征运用密度峰值算法提取数据的密度和距离特征，该余弦相似性特征也考虑到相邻bin之间的联系，对拷贝数变异的检测有很大的帮助，可作为改进方案对检测结果进一步优化。

本发明首先得到数据的密度和距离两个特征，对其做高斯分布统计，然后根据数据的概率密度来判断其是否发生变异，可以提高检测结果的准确率和灵敏度。本发明在技术上使用密度峰值算法得到数据的密度和距离两个特征，通过两个特征将变异区域与正常区域分隔开，可以更容易，更准确的检测出变异位置。本发明所使用的密度峰值算法，简洁有效，可以一次性地将变异区域与正常区域分开，计算所需要的时间和空间都不大。

本发明解决了因为测序覆盖度低导致变异区域的读段深度与正常区域的读段深度差异很小的问题。通过使用密度峰值算法将变异区域的特征凸显出来，并检测出变异的位点；本发明使用密度峰值算法得到数据的密度和距离特征，对两个特征做高斯分布统计得到数据的概率密度，通过概率密度进行判断，可以提高检测结果的各项指标；本发明的密度峰值检测算法可以快速将样本的簇中心、正常区域和变异区域分开，不需要迭代或者其他多余的处理。

本发明的方法简单易实现，资源开销少，通过给定的bam文件和参考序列文件即可检测出最后的结果，密度峰值算法可以快速对数据进行处理，然后本发明再对数据做高斯分布统计，根据概率密度判断从而得到变异区域。

表1是本发明与其他拷贝数变异检测方法对仿真数据检测的结果对比。

3.充分考虑了拷贝数相邻bin之间存在的联系，提出一种数据分割方法，将我们的分析单位从bin转化为seg。

4.根据拷贝数变异的特性(变异点被孤立，而正常点聚集)，运用密度峰值算法提取出数据的密度和距离这两个与拷贝数密切相关的特征，然后进一步分析。

5.本发明还可以利用数据的余弦相似性特征代替位置特征，然后结合readdepth特征运用密度峰值算法提取数据的密度和距离特征，该余弦相似性特征也考虑到相邻bin之间的联系，对拷贝数变异的检测有很大的帮助，可作为改进方案对检测结果进一步优化。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的拷贝数变异检测方法流程图。

图2是本发明实施例提供的拷贝数变异检测系统的结构示意图；

图中：1、数据预处理模块；2、segment计算模块；3、数据计算模块；4、筛选模块；5、变异类型确定模块。

图3是本发明实施例提供的拷贝数变异检测方法的实现流程图。

图4是本发明所使用到的数据是由经典数据仿真软件intsim得到的，纯度为0.2，0.4和0.6，覆盖度为4x和6x，一共六组数据，每组50个样本。为了直观的展现我们的算法对数据处理的结果，以一个样本为例的效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种拷贝数变异检测方法、系统、存储介质、计算机设备，下面结合附图对本发明作详细的描述。

如图1所示，本发明提供的拷贝数变异检测方法包括以下步骤：

s101：从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

s102：对经过处理后的数据进行segment计算；

s103：使用密度峰值检测算法对segment计算产生的数据进行计算；

s104：得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

s105：得到变异区域，确定变异是丢失还是增添。

如图2所示，本发明提供的拷贝数变异检测系统包括：

数据预处理模块1，用于从bam文件中提取出每个窗口的readdepth值，对数据进行gc校正；

segment计算模块2，用于对经过处理后的数据进行segment计算；

数据计算模块3，用于使用密度峰值检测算法对segment计算产生的数据进行计算；

筛选模块4，用于得到标准化以后的密度和距离的二维数据，通过每个seg的密度和距离两个特征筛选出变异区域；

变异类型确定模块5，用于得到变异区域，确定变异是丢失还是增添。

下面结合附图对本发明的技术方案作进一步的描述。

如图3所示，本发明提供的拷贝数变异检测方法具体包括以下步骤：

(1)数据预处理

从bam文件中提取出每个窗口的readdepth值，因为序列中gc的含量会影响到readdepth值，所以本发明需要对数据进行gc校正，然后就可以对数据进行segment计算。

(2)segment计算

对经过第(1)步处理后的数据进行segment计算，segment计算是采用cghflasso分割算法将数据平滑处理。通过此步计算，可将散乱的数据分割成长度不等的几个片段，本发明称之为seg，相对于原始的数据，可以更好地检测出变异的位点。然后将segment后的一维readdepth数据加上位置信息转化成二维数据，由于每个位点的位置差异为1，而readdepth的差异较小，所以将两个特征压缩到同一个长度空间中来减少检测出现的误差。

(3)使用密度峰值检测算法对数据进行计算

使用密度峰值检测算法对第(2)步产生的数据进行计算，首先计算每个seg的与其他seg的距离，得到距离矩阵dis_matrix，然后通过dis_matrix和dc参数求出每个seg的密度density，(密度density即以某seg为圆心，dc为半径范围内其他seg的个数)。紧接着通过dis_matrix和density求出簇中心的seg集合和每个seg对应的距离distance。倘若sega为全局密度最大的seg，则a的距离为a与离a最远的seg之间的的距离，否则，a的距离为a与密度比a大且离a最近的seg之间的距离。对每个seg进行这样的步骤，得到所有seg的密度和距离特征。通过密度峰值检测算法将第(2)步产生的一维readdepth数据转化为二维的特征矩阵，由于两个特征的值不在一个范围内，所以需要对两个数据集进行标准化处理。

(4)筛选结果

通过第(3)步的计算，得到标准化以后的密度和距离的二维数据，密度和距离都比较大的seg即为簇中心，而密度小但距离大的seg即为异常区域。由于每个样本不同，本发明不能人为设定两个阈值进行判断，所以首先对该二维数据进行多元正态分布，然后根据该多元正态分布的概率密度函数求出每个seg的概率密度，最后本发明设定一个可调整的概率密度阈值prob，根据每个seg的概率密度判断该seg是否发生变异，若某seg的概率密度小于prob并且其密度小于全局平均密度，本发明认为该seg是一个异常片段。这样本发明通过每个seg的密度和距离两个特征就可以筛选出变异区域。

(5)确定变异类型

通过第(4)步的计算，得到了变异区域，还需要确定该变异是丢失还是增添。而对于如何确定变异类型，该方法认为，簇中心的区域为正常区域，本发明取簇中心的readdepth值的平均值作为标准值(因为簇中心可能不只一个，所以取平均值)。若该变异区域的readdepth值大于标准值，则认为该区域为拷贝数扩增区域(duplication)，若该变异区域的readdepth值小于标准值，则认为该区域为拷贝数丢失区域(deletion)。

本发明通过位置特征和readdepth特征得到每个seg的密度和距离特征，然后通过密度和距离进行操作。这里本发明还可以使用余弦相似性特征代替位置特征来计算每个seg的密度和距离。余弦相似性是通过位置和readdepth二维数据通过计算得到的。假设有从左到右y，x，z三个点，当x是异常点时，向量y-x和z-x之间的夹角不会因为不同的y和z值而剧烈变化，而且夹角会随着点间距的增加而减少，由于夹角在0到180°的范围内，余弦值单调，所以不同异常程度的点会对应不同的余弦值，因此可以用余弦相似性代替位置特征。

应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、cd或dvd-rom的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。在本文最后的表1中还提供了本发明与其他几种检测方法对相同仿真数据的检测结果的对比，对于低覆盖度的数据，本发明检测的结果指标明显优于其他方法，并且在灵敏度指标方面表现优越。

本发明所使用到的数据是由经典数据仿真软件intsim得到的，纯度为0.2，0.4和0.6，覆盖度为4x和6x，一共六组数据，每组50个样本。为了直观的展现我们的算法对数据处理的结果，以一个样本为例，如图4所示。

该图中大部分的点都集中在右下方，这些都是正常区域。右上角的两个点，其密度是最大的，且距离也大，因为它们是簇中心，也是正常区域。而左边靠上的点是异常点的概率很大，属于异常区域，靠下的点不太好判断，我们需要根据点的概率密度进行判断。但是，显然异常区域比较稀疏，概率密度偏小，正常区域点很密集(由于簇中心所在区域也很稀疏，所以判断时我们还需要根据密度排除掉将簇中心判断为异常点的可能)，概率密度也相应偏高。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。