1.本发明属于生物信息检测技术领域，具体涉及一种基因拷贝数变异检测方法及装置、设备、存储介质。
背景技术：2.宏基因临床检测技术是利用二代高通量测序从遗传物质角度对微生物感染进行鉴定和诊断的新型临床技术。其中高通量测序包含的微生物序列(即非人源性序列)相对有限，超过90％的序列内容为人源性序列。但至今的临床应用中却仅使用了量较少的微生物序列进行感染学鉴定，量较多的人源性序列依然缺少应用和分析。宏基因技术利用游离核酸进行微生物鉴定，而肿瘤细胞因其高代谢率容易在游离核酸中检测到，游离核酸已经成为肿瘤筛查和检查的重要实验样本。如果通过宏基因在感染学检测的同时进行肿瘤早筛检测，不但能够解决当前疾病问题，还能对肿瘤早筛阳性结果提前介入和治疗，对人类健康有着巨大意义。但至今却极少联合宏基因和肿瘤基因拷贝数变异(copy number variations，cnv)检测的算法和软件。因此通过宏基因感染的检测手段能够同时进行肿瘤筛查的工作受技术限制无法进行。
3.目前针对cnv检测的主流理论包括read-pair(rp)法、split-read(sr)法、read-depth(rd)法和assembly(as)法等4种方法。其中rp是最早出现的算法，利用双端测序插入片段长度分布来检测cnv，也称之为pair end mapping(pem)方法。当插入片段长度过长或者过短时，都代表着基因组发生了结构变异。sr方法利用一端能够比对，另外一端比对不上的reads来识别cnv。另外一端比对不上，可能是存在cnv,通过将单独的reads进行拆分，使其能够正确比对到参考基因组上，拆分的点就是cnv的断裂点。rd方法利用拷贝数和对应区域测序深度的相关性来进行分析，基本模型是缺失区域的测序深度相对低，而插入区域的测序深度相对高。as方法利用测序得到的短序列进行组装，将组装的contig与参考基因组进行比较，从而确定发生了结构变异的区域。
4.由于宏基因检测为50-75bp的二代单端测序技术，与当前主流软件均不相符，同时地区人群位点、不同实验条件引入的位点差异均无法得到很好的矫正，使得其分辨率和结果在目前的软件中均表现不佳，至今仍无可用于宏基因cnv检测的算法。针对cnv检测主流理论，rp和sr两种方法依赖于双端测序技术，不适用于宏基因数据，且相对算法不够精确。as法依赖通量和测序覆盖度技术，宏基因的覆盖度与基因组组装技术相差甚远，该方法无法运行于宏基因数据。这三种理论方法均与宏基因测序数据不符合，仅rd法下的少量软件可以应用。
5.但是rd法依赖通量和测序深度技术，需要较高较稳定的深度变化才能用于识别cnv，因此在宏基因这种低测序深度的数据中应用该方法会引入极多的假阳性位点。且传统rd法分析模型使用相同的rd分辨率和cnv分辨率，过小的分辨率会导致数据离散型太强，假阳性过高，过大的分辨率会导致cnv平均化，产生假阴性结果，同时cnv的边缘位置也可能因为rd计算时的区间覆盖形成过渡型，影响cnv的识别和判断，因此导致cnv检测的准确性、灵
敏度均不够。
技术实现要素：6.本发明的目的在于提供一种基因拷贝数变异检测方法及装置、设备、存储介质，可以针对rd位点和cnv片段使用不同的分辨率，在提高rd位点分辨率的同时保证cnv片段的稳定性，提高cnv检测的准确性和灵敏度。
7.本发明第一方面公开一种基因拷贝数变异检测方法，包括：
8.从待测样本的测序数据中确定出多个指定位点的测序深度值；其中，每一所述指定位点对应有置信区间及rd设定均值；
9.若任一指定位点的测序深度值未位于该指定位点对应的置信区间内，确定该指定位点为拷贝异常位点；
10.根据各个所述拷贝异常位点对应的rd设定均值以及测序深度值，计算各个所述拷贝异常位点的拷贝数值；其中，所述拷贝数值与所述测序深度值呈正相关关系；
11.将所有所述拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，所述拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类；
12.将归属于所述高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于所述低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点；
13.将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。
14.本发明第二方面公开一种基因拷贝数变异检测装置，包括：
15.深度确定单元，用于从待测样本的测序数据中确定出多个指定位点的测序深度值；其中，每一所述指定位点对应有置信区间及rd设定均值；
16.异常检测单元，用于在任一指定位点的测序深度值未位于该指定位点对应的置信区间内时，确定该指定位点为拷贝异常位点；
17.计算单元，用于根据各个所述拷贝异常位点对应的rd设定均值以及测序深度值，计算各个所述拷贝异常位点的拷贝数值；其中，所述拷贝数值与所述测序深度值呈正相关关系；
18.聚类单元，用于将所有所述拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，所述拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类；
19.变异确定单元，用于将归属于所述高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于所述低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点；
20.合并单元，用于将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。
21.本发明第三方面公开一种电子设备，包括存储有可执行程序代码的存储器以及与所述存储器耦合的处理器；所述处理器调用所述存储器中存储的所述可执行程序代码，用
于执行第一方面公开的基因拷贝数变异检测方法。
22.本发明第四方面公开一种计算机可读存储介质，所述计算机可读存储介质存储计算机程序，其中，所述计算机程序使得计算机执行第一方面公开的基因拷贝数变异检测方法。
23.本发明的有益效果在于，所提供的基因拷贝数变异检测方法及装置、设备、存储介质，方法主要包括：利用每个指定位点的置信区间确定出拷贝异常位点，计算拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，其中拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类，此时拷贝异常位点能够利用空间聚类得到很好的分类，进而可以实现对高拷贝变异位点和低拷贝变异位点的自动识别，以及过滤掉归属于正常拷贝类的拷贝异常位点，减少误判率；然后，将归属于高拷贝变异子类且测序深度值大于相应的置信区间、以及归属于低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点确定为拷贝变异位点；接着将位置相邻且同变异类型的拷贝变异位点合并得到更加精确的cnv片段。
24.因此本发明通过分离rd位点的区间大小和cnv片段的区间大小，针对rd位点和cnv片段使用不同的分辨率，利用提高rd位点的稳定性连锁提高cnv片段的稳定性，这样可以在提高rd位点分辨率的同时保证cnv片段的稳定性，防止大片段深度不稳定、单位点离散度太大等数据结构问题，同时使用无监督机器学习方法对拷贝异常位点进行聚类分类和与置信区间比对的方法来双重识别拷贝变异位点，可以进一步提高拷贝变异位点的识别精确性，进而提高cnv检测的准确性、灵敏度和稳定性。
附图说明
25.此处的附图，示出了本发明所述技术方案的具体实例，并与具体实施方式构成说明书的一部分，用于解释本发明的技术方案、原理及效果。
26.除非特别说明或另有定义，不同附图中，相同的附图标记代表相同或相似的技术特征，对于相同或相似的技术特征，也可能会采用不同的附图标记进行表示。
27.图1是本发明实施例公开的一种基因拷贝数变异检测方法的流程图；
28.图2是本发明实施例公开的一种宏基因测序数据的基本结构图；
29.图3是本发明实施例与其他cnv软件在相同分辨率下的检测效果比对图；
30.图4是本发明实施例与其他cnv软件在最佳分辨率下的检测效果比对图；
31.图5是本发明实施例与其他cnv软件针对阳性样本的检测效果比对图；
32.图6是本发明实施例公开的一种基因拷贝数变异检测装置的结构示意图；
33.图7是本发明实施例公开的一种电子设备的结构示意图。
34.附图标记说明：
35.601、深度确定单元；602、异常检测单元；603、计算单元；604、聚类单元；605、变异确定单元；606、合并单元；701、存储器；702、处理器。
具体实施方式
36.除非特别说明或另有定义，本文所使用的所有技术和科学术语与所属技术领域的技术人员通常理解的含义相同。在结合本发明的技术方案以现实的场景的情况下，本文所
使用的所有技术和科学术语也可以具有与实现本发明的技术方案的目的相对应的含义。本文所使用的“第一、第二
…”
仅仅是用于对名称的区分，不代表具体的数量或顺序。本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。
37.除非特别说明或另有定义，本文所使用的“所述”、“该”为相应位置之前所提及或描述的技术特征或技术内容，该技术特征或技术内容与其所提及的技术特征或技术内容可以是相同的，也可以是相似的。
38.毫无疑义，与本发明的目的相违背，或者明显矛盾的技术内容或技术特征，应被排除在外。为了便于理解本发明，下面将参照说明书附图对本发明的具体实施例进行更详细的描述。
39.需要说明的是，本发明公开的基因拷贝数变异检测方法及装置，可适用于如宏基因数据的低测序深度数据，亦可适用于正常测序的高测序深度高数据，比如全外显子测序数据、全基因测序数据、芯片测序数据。在本发明实施例中，以宏基因测序数据为例进行阐述，不应认为是对本发明的限定。
40.(一)基线的建立过程
41.优选的在本发明实施例中，在对待测样本的测序数据进行拷贝数变异检测之前，可利用大量同批次或条件、来源的临床样本数据结合机器学习的方法和方差指标过滤不稳定、不合理的噪音位点，利用k-mean算法识别位点中的着丝粒、线粒体和重复序列位点，并且结合方差剔除波动性位点，保留稳定可用于cnv检测的位点，然后建立各个稳定位点的基线。
42.具体的，在对待测样本的测序数据进行拷贝数变异检测之前，可以执行以下步骤s01～s07：
43.s01、获取n1个训练样本的测序数据；每个训练样本的测序数据均包括m1个候选位点的测序深度值。
44.在本发明中，位点实则为位点区间，指的是染色体区段，测序深度值则对应为染色体区段的read-depth(rd)数据。需要说明的是，在对肿瘤进行早筛的时候，与肿瘤相关的大部分cnv片段通常发生在常染色体上，本基因拷贝数变异检测方法主要是应用于常染色体上cnv片段的检测场景，但不限于此，在一些可能的实施例中也可以应用于性染色体上cnv片段的检测场景。
45.由于执行的是无监督的聚类方式，获取的测序数据可仅包括原始的fastq格式的测序数据，不需要临床验证和随访数据，例如训练样本是属于cnv检测阳性或阴性的结果数据。因此本发明的应用场景可以为双盲，对早筛意义更大。在获取n1个训练样本的原始的fastq格式的测序数据之后，还进一步根据预先设计的对比和质控流程对原始的fastq格式的测序数据进行过滤、去接头、比对、去重等预处理工作。
46.然后按照预设的block区间梯度参数，确定位点区间(即染色体区段)大小，根据该位点区间大小将测序数据中的染色体划分为m1个染色体区段，每一个染色体区段则对应一个候选位点。例如，假设预设的block区间梯度参数为10kb，则可得到该流程下稳定的位点区间大小为10kb，基于此将染色体划分为m1个染色体区段，每一染色体区段的大小为10kb。
47.在获取各个训练样本的m1个候选位点之后，优选的还可以对n1个训练样本的m1个候选位点的测序深度值进行数据标准化处理，使得同一候选位点下n1个训练样本中的测序
深度值相差较大的不同特性数据落在设定范围[0，1]内，从而消除不同基因组批次和样本测序深度不同的批次差异，减少测序数据中取值过大或者取值过小的奇异特性数据导致的不良影响。
[0048]
s02、对所有候选位点的测序深度值进行无监督聚类分类，获得多个分类类别。
[0049]
步骤s02中则对n1
×
m1个候选位点的测序深度值进行区间质控，主要通过利用k-mean算法识别分类树，以将n1
×
m1个候选位点的测序深度值进行无监督聚类分类成多个分类类别。
[0050]
s03、计算每一分类类别的第一rd均值，根据第一rd均值从多个分类类别中识别出噪音类别。
[0051]
步骤s03中对各个分类类别所包括的所有候选位点的测序深度值求平均，获得每一分类类别的第一rd均值，然后基于第一rd均值，可识别出均值异常的分类类别，从而将均值异常的分类类别确定为噪音类别剔除。从而可以通过机器学习聚类算法识别线粒体位点、重复序列位点、着丝粒位点等背景噪音位点，实现高效、细致化降噪。
[0052]
具体的，将第一rd均值小于第一指定阈值的分类类别确定为着丝粒类别；以及，将第一rd均值大于第二指定阈值的分类类别确定为重复序列类别；以及，将第一rd均值大于第三指定阈值的分类类别确定为线粒体类别；而将着丝粒类别、重复序列类别和线粒体类别均作为噪音类别。
[0053]
其中，第一指定阈值小于第二指定阈值，第二指定阈值小于第三指定阈值。第一指定阈值、第二指定阈值、第三指定阈值可由开发人员根据实际需求而预先设定，而在一些优选的实施例中，上述第一指定阈值、第二指定阈值、第三指定阈值亦可根据多个分类类别的整体测序深度平均值确定。
[0054]
比如，在计算每一分类类别的第一rd均值之后，还可以进一步对所有分类类别的第一rd均值求平均，获得所有分类类别整体性的第二rd均值，然后将第二rd均值的十分之一作为第一指定阈值；以及，将第二rd均值的五倍作为第二指定阈值；以及，将第二rd均值的二十倍作为第三指定阈值。也即是说，假设所有分类类别的第二rd均值为cu，则第一指定阈值为cu/10、第二指定阈值为cu*5、第三指定阈值为cu*20。
[0055]
s04、将噪音类别所包括的候选位点剔除，获得m2个设定位点。m2小于m1。
[0056]
在聚类结果中，着丝粒类别包括的候选位点因为测序难度大是典型的噪音位点，且一般与肿瘤cnv变异无关。重复序列类别包括的候选位点拷贝数较高且在基因组位置不稳定，不具有明确的诊断意义，也作为噪音位点去除；线粒体类别包括的候选位点在不同的细胞中差异较大，在低测序深度数据中会成为重要噪音点，因此在本发明中将这些候选位点剔除后，所保留的m2个设定位点，包括但不限于位于y染色体上的设定位点、位于x染色体上的设定位点和位于常染色体上的设定位点，作为稳定的cnv检测靶点。
[0057]
s05、计算每一设定位点在n1个训练样本下的rd设定均值及方差。
[0058]
具体通过计算n1个训练样本分别对应于每一个设定位点的测序深度值(优选为标准化后的rd数据)的平均值，作为该设定位点对应的rd设定均值，每个设定位点的rd设定均值之间存在差异性；以及计算n1个训练样本分别对应于每一个设定位点的测序深度值的方差。
[0059]
s06、从m2个设定位点中确定出全部或部分设定位点作为指定位点。
[0060]
在确定出m2个设定位点之后，在其它一些可能的实施例中，可以将全部的设定位点均作为指定位点，即用于后续对待测样本进行检测时参照的稳定位点。而在本发明实施例中，优选对m2个设定位点再过滤一遍，在样本维度上对方差较大的少量设定位点进行剔除。具体的，在计算得到每一设定位点在n1个训练样本下的方差之后，可以将方差小于指定方差阈值的设定位点确定为指定位点，以获得m3个指定位点，而将方差大于或等于指定方差阈值的设定位点进行剔除。
[0061]
需要说明的是，m3小于或等于m2。当m2个设定位点的方差均小于指定方差阈值时，m3个指定位点包括全部的m2个设定位点，而当m2个设定位点中存在部分设定位点的方差大于或等于指定方差阈值时，则m3个指定位点包括部分的设定位点。
[0062]
s07、根据rd设定均值及方差，计算每一指定位点的置信区间。
[0063]
确定m3个指定位点之后，则通过以下公式计算每一指定位点的置信区间：
[0064][0065]
式中，代表第i个指定位点对应的rd设定均值，s
i2
代表第i个指定位点对应的方差，由于数据符合正态分布，为1.96。
[0066]
因此，第i个指定位点的置信区间表达式为
[0067]
也即是说，利用n1个训练样本计算完毕后，每个指定位点都对应有基线，该基线的内容包括置信区间及rd设定均值，置信区间用于规定该指定位点的测序深度值的正常值范围，落在置信区间外的点将会被识别为异常值。通过设置ci＝95％为指定位点的置信区间，也即为每个指定位点建立专属的残差范围，且随着样本量累积，残差范围更加准确，且不依赖已知结局样本。
[0068]
实施步骤s01～s07，通过利用识别的相对稳定的设定位点建立基线，利用基线去除不同位点之间的差异，去除gc含量、人群基因组来源、实验方法和质控方法等差异引入的位点间深度不平均的噪音点，能够很好对背景位点进行异质性处理，更好的去除批次效应和地方性人源特点带来的噪音，保证识别的拷贝数变异是真实的变化。
[0069]
(二)对待测样本/训练样本的检测过程
[0070]
如图1所示，本发明实施例公开一种基因拷贝数变异检测方法，包括以下步骤s10～s60：
[0071]
s10、从待测样本的测序数据中确定出多个指定位点的测序深度值。
[0072]
在本发明中，指定位点指的是预先指定的染色体区段，在其它一些可能的实施例中具体可由开发人员根据实际需求指定，而在本实施例中则是参照上述基线的建立过程中确定出的m3个指定位点进行确定的。在实际测序中，待测样本不一定测得指定的m3个指定位点中的全部位点，因此在步骤s10中，所确定出的待测样本的多个指定位点包括上述m3个指定位点中的全部或部分位点。然后可根据确定出待测样本的各个指定位点，调取出与各个指定位点对应的基线内容，即对应的置信区间及rd设定均值。
[0073]
s20、若任一指定位点的测序深度值未位于该指定位点对应的置信区间内，确定该指定位点为拷贝异常位点。
[0074]
其中，将各个指定位点的测序深度值与其对应的置信区间进行比对，只要测序深度值未位于该指定位点对应的置信区间内，则确定该指定位点为拷贝异常位点，遍历多个指定位点之后，可确定出多个拷贝异常位点。具体的，若测序深度值小于相应的置信区间，确定该指定位点为低拷贝异常位点，从而获得低拷贝异常位点集合；若测序深度值大于相应的置信区间，确定该指定位点为高拷贝异常位点，从而获得高拷贝异常位点集合。
[0075]
举例来说，假设待测样本的第1个指定位点的测序深度值为rd1，若则识别第1个指定位点为低拷贝异常位点；若则识别第1个指定位点为高拷贝异常位点。
[0076]
上述高、低拷贝异常位点可能为拷贝数变异位点或者生态分析下的误差点(α＝0.05)，单个异常位点不太明确，为了进一步提高精确性，可计算上述识别出异常的拷贝异常位点的拷贝数值进一步分析，即执行步骤s30。
[0077]
s30、根据各个拷贝异常位点对应的rd设定均值以及测序深度值，计算各个拷贝异常位点的拷贝数值。其中，拷贝数值与测序深度值呈正相关关系。
[0078]
在本实施例中，通过以下公式计算各个拷贝异常位点的拷贝数值：
[0079][0080]
式中，cpj代表第j个拷贝异常位点的拷贝数值，rdj代表第j个拷贝异常位点的测序深度值，代表第j个拷贝异常位点对应的rd设定均值。
[0081]
需要说明的是，各个拷贝异常位点的拷贝数值不限于上式(2)所述的将拷贝异常位点的测序深度值与rd设定均值之间的比值作为该拷贝异常位点的拷贝数值的计算方式，在其它一些可能的实施例中，计算拷贝数值的方式亦可采用其它的与上式(2)等同的替代公式或明显变型公式，同样能保证拷贝数值与测序深度值呈正相关关系。例如等。
[0082]
s40、将所有拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，其中拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类。
[0083]
具体的在本实施例中，样本结论无预知性，利用无监督聚类进行识别，预选范围k值为从1到10，交叉验证法搜索最优k值，经过样本大量验证，先验聚类k＝4时最佳。因此在本实施例中，设置k＝4，利用k近邻法(k-nearest neighbors，knn)自然聚类所有拷贝异常位点的拷贝数值，可获得4个分类类别。然后计算各个分类类别的拷贝数值(cp)平均值，按照cp平均值对4个分类类别进行排序，按照cp平均值从高到低的顺序分别鉴定为高拷贝变异子类、高拷贝离散类、低拷贝离散类、低拷贝变异子类。其中，高拷贝离散类、低拷贝离散类均可视为正常拷贝类，用于过滤掉归属于正常拷贝类的拷贝异常位点，减少误判率。
[0084]
s50、将归属于高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点。
[0085]
接着可将识别出的高拷贝变异子类中包括的拷贝异常位点与上述与置信区间比对得到的高拷贝异常位点集合取第一交集，第一交集中的拷贝异常位点则确定为高拷贝变
异位点。同样的，可将识别出的低拷贝变异子类中包括的拷贝异常位点与上述与置信区间比对得到的低拷贝异常位点集合取第二交集，第二交集中的拷贝异常位点则确定为低拷贝变异位点。其中，高拷贝变异位点与低拷贝变异位点统称为拷贝变异位点。
[0086]
通过上述取交集的方式，即同时使用k近邻法数据聚类方法对拷贝异常位点进行聚类分类和与基线置信区间比对的方法，来双重识别拷贝变异位点，可以进一步提高拷贝变异位点的识别精确性，同时拷贝数变异位点能够利用空间的聚类得到很好的聚类，进而实现对拷贝数增加和减少等突变类型的自动识别，更精准识别高拷贝变异位点和低拷贝变异位点，进而可以提高cnv检测的准确性、灵敏度和稳定性。
[0087]
s60、将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。
[0088]
其中，将位置相邻且同属于高拷贝变异子类的高拷贝变异位点进行合并，获得高拷贝变异片段；以及，将位置相邻且同属于低拷贝变异子类的低拷贝变异位点进行合并，获得低拷贝变异片段。高拷贝变异片段与低拷贝变异片段统称为拷贝变异片段。
[0089]
具体的，对于相同变异类型(即同属于高拷贝变异子类或低拷贝变异子类)的拷贝变异位点在染色体内部进行连续性合并，即将位置相邻且变异类型相同的染色体区段进行合并。其中，采用滑动窗口在多个拷贝变异位点上进行扫描，滑动窗口的长度大小应当至少涵盖2个拷贝变异位点，本实施例中每个位点的区间大小为10kb，因此滑动窗口的长度大小可设为20kb，且该滑动窗口的移动步长为10kb，即滑动窗口每次移动1个拷贝变异位点的距离，窗口内始终有2个拷贝变异位点，在滑动窗口每次移动后，判断滑动窗口内的后一个拷贝变异位点是否与其相邻的前一个拷贝变异位点归属于同一个拷贝变异子类，若是，将该后一个拷贝变异位点与其相邻的前一个拷贝变异位点进行合并。在扫描过程中，取允许异常点为1/100kb步长值延伸，即每100kb步长中允许存在1个非拷贝变异位点依然被判定为连锁位点，与拷贝变异位点合并到同一拷贝变异片段中。
[0090]
举例来说，假设多个拷贝异常位点为a1、a2、a3、
…
、a8、a9、
…
、aj，滑动窗口首先选择{a1、a2}，若a2与a1归属于同一个拷贝变异子类，则将两者合并；接着滑动窗口移动10kb步长，移动至{a2、a3}，若a3与a2归属于同一个拷贝变异子类，则将两者合并；以此类推直至遍历所有拷贝异常位点。
[0091]
在此过程中，假设从a1合并至a8，接着滑动窗口移动至{a8、a9}，若判断出a9与a8归属于不同的拷贝变异子类，而此时滑动窗口自a1开始所移动的步长合计并未超过100kb，可将a9视为该100kb步长内第一个出现的非拷贝变异位点，依然会将a9视为连锁位点与a8进行连锁合并；接着若判断出a
10
与a9归属于不同的拷贝变异子类，则判定当前合并的拷贝变异片段中断，不会合并a
10
，而若判断出a
10
与a9归属于同一个拷贝变异子类，则继续将a
10
与a9进行连锁合并。在此例子基础上，若连续合并至a
j-1
，且判断出aj与a
j-1
归属于不同的拷贝变异子类，则需要判断aj与上一个出现的非拷贝变异位点a9之间的距离，若aj与a9之间的距离超过100kb，则可依然将aj视为连锁位点与a
j-1
进行连锁合并；否则，若aj与a9之间的距离不超过100kb，则判定当前合并的拷贝变异片段中断。
[0092]
可见实施本发明实施例，当多个同变异类型的拷贝变异位点的位置相邻时，利用连锁算法延伸位点，将位置相邻且同变异类型的拷贝变异位点合并，以得到更加精确的cnv片段，而且通过分离rd位点的区间大小和cnv片段的区间大小，针对rd位点和cnv片段使用
不同的分辨率，然后利用提高rd位点的稳定性连锁提高cnv片段的稳定性，这样可以在提高rd位点分辨率的同时保证cnv片段的稳定性，防止大片段深度不稳定、单位点离散度太大等数据结构问题，同时使用无监督机器学习方法对拷贝异常位点进行聚类分类和与基线置信区间比对的方法来双重识别拷贝变异位点，可以进一步提高拷贝变异位点的识别精确性，进而提高cnv检测的准确性、灵敏度和稳定。
[0093]
进一步的，对于计算后的所有拷贝变异片段，可设置拷贝变异片段阈值进行筛选，从而进一步提高精确性。具体优选的，在本发明实施例中，执行步骤s60之后，还可以执行以下步骤s70或者s80：
[0094]
s70、若拷贝变异片段的大小达到第一片段阈值，确定拷贝变异片段为目标拷贝变异片段。
[0095]
s80、若拷贝变异片段的大小达到第二片段阈值且小于第一片段阈值，确定拷贝变异片段为疑似拷贝变异片段。
[0096]
例如，设置第一片段阈值为10mb，第二片段阈值为1mb，也即大小达到1mb且小于10mb的拷贝变异片段为疑似拷贝变异片段，大小达到10mb的拷贝变异片段为目标拷贝变异片段。
[0097]
执行步骤s70之后，还可以执行以下步骤s71～s72：
[0098]
s71、获取目标拷贝变异片段包括的各个目标拷贝异常位点的拷贝数值，通过以下公式对目标拷贝异常位点的拷贝数值进行log转换：
[0099]
cp
′z＝log(cpz+0.001)
ꢀꢀꢀ
(3)
[0100]
式中，cp
′z代表作图引用的拷贝参数，cpz代表第z个目标拷贝异常位点的拷贝数值。
[0101]
s72、根据各个目标拷贝异常位点的作图引用的拷贝参数按照染色体分布的顺序进行绘图。
[0102]
通过执行步骤s71～s72，对目标拷贝异常位点的拷贝数值log转换后得到稳定、可读性高的作图引用的拷贝参数，在对识别的目标拷贝变异片段添加辅助线进行提示，同时借助作图工具实现对数转化后的位点拷贝数作图可视化。
[0103]
在执行步骤s01～s07之后，以上步骤s10～s80为针对未知的单个待测样本的检测步骤，而在实际应用中还可以采用未知的批次测试样本进行实施，具体在步骤s10中针对的样本及样本数量有所不同，在对测试样本的检测过程中，首先获取n2个测试样本的测序数据；从每个测试样本的测序数据中确定出多个指定位点的测序深度值，然后再按照以上步骤s20～s80实施，可实现测试功能，同时对步骤s01～s07确定出的每一指定位点及其基线进行优化更新。
[0104]
例如按照步骤s40对n2个测试样本实施之后，所获得的拷贝变异类中的拷贝异常位点均会得到保留并作为新的指定位点，同时每个新的指定位点会利用历史样本(至少包括n1个训练样本)数据和本批次的n2个测试数据再次计算新的基线，并将各个新的指定位点对应的新的基线以及本批次的测试数据存储至数据库。以便下次对未知的单个待测样本进行检测或者对下批次的测试样本数据进行指定位点及其基线优化。
[0105]
上述涉及的测试样本可以大量来源于临床检验数据，不需要额外的科研型训练集和临床信息，临床检验样本多、数据大的优势都可以转换为指定位点及其基线的稳定性，通
过持续送检对指定位点及其基线的优化，可以针对检验单位的检验方式、检验人群、分析方式、测序方法有更好的兼容性和稳定性。
[0106]
如图2所示，图2中展示了常染色体、x/y性染色体、线粒体基因组的gc含量和深度情况，其中可见染色体存在gc含量不稳定异常点，这些点是是分析中引入假信号的异常点。从深度数据可见基因组整体测序深度很低，在1到2之间，不满足现有技术中组装等cnv识别技术的应用。
[0107]
图3为相同分辨率下结果比对，上半图为现有其他cnv识别软件，下半图为本发明实施例，可见上半图在相同分辨率下位点离散严重，不同位点的差异性未能去除，同时上下部都有少量离群点，无法判断为cnv或背景噪音。下半图为本发明实施例的检测结果，可见染色体标记明显，可读性更强，而且所有的位点正态分布于正常二倍拷贝线(中间实线)上下，位于中间实线上下的两个虚线分别表示拷贝数值为三倍拷贝线和单倍拷贝线，可见无明显cnv变异，该样本为阴性样本。
[0108]
图4为现有其他cnv识别软件在最佳分辨率下的分辨情况。可见该软件针对同一样本不同分辨率下的结果差异明显，上半图结果在单倍拷贝线下方存在大量离散点，且每个点表示1m的大片段突变，即上半图表示为样本阳性，但人工识别可见，该异常位点大多为低深度异常点，来源于低深度数据或低质量位点，而经过该软件识别形成假阳性结果。下半图为本发明实施例的检测结果，可见单倍拷贝线下方无明显离散点，且所有位点分布稳定在正常区间(三倍拷贝线和单倍拷贝线之间)内，可见为明显的正常的阴性样本。
[0109]
图5为阳性样本的分辨情况。现有其他cnv识别软件使用最优参数进行分析，可见13号染色体存在拷贝数加倍变异，但不同位点存在差异，特别在变异区域边缘的位点呈现离散过度点，使得整个变异区间位置无法确认，变异倍数不够明确，14号染色体仅有个别点呈现缺失突变信号，但数量少，是否为变异尚无法判断。下图为本发明实施例的检测结果，可见13号染色体突变结果十分明确，大段染色体呈现拷贝数增加现象，同时拷贝数中心介于三倍拷贝线和二倍拷贝线之间，即识别到细胞间变异存在异质性，提示可能存在肿瘤而非遗传病导致的cnv信号，同时14号染色体的缺失突变异常明显，本发明的机器学习工具也识别了该位点的突变并辅助标记。
[0110]
综上可见，本发明与有限的其他方法明显有着更好的特异性和稳定性。本发明能够过滤低深度和重复序列位点，而且能够利用庞大的位点基线去除批次效应、实验方法、质控手段、地区人群特点等带来的位点深度差异，通过连锁计算能够从基因组位置维度提高cnv检测的稳定性并绘图，因此能够很好适用宏基因测序等低深度测序数据结构下的cnv识别，可用于感染患者中肿瘤和遗传病的探索和早筛，临床意义重大。
[0111]
而且，实施本发明实施例，通过分离rd位点的区间大小和cnv片段的区间大小，针对rd位点和cnv片段使用不同的分辨率，利用提高rd位点的稳定性连锁提高cnv片段的稳定性，这样可以在提高rd位点分辨率的同时保证cnv片段的稳定性，防止大片段深度不稳定、单位点离散度太大等数据结构问题，同时使用无监督机器学习方法对拷贝异常位点进行聚类分类和与置信区间比对的方法来双重识别拷贝变异位点，可以进一步提高拷贝变异位点的识别精确性，进而提高cnv检测的准确性、灵敏度和稳定性。
[0112]
如图6所示，本发明实施例公开一种基因拷贝数变异检测装置，包括深度确定单元601、异常检测单元602、计算单元603、聚类单元604、变异确定单元605和合并单元606，其
中，
[0113]
深度确定单元601，用于从待测样本的测序数据中确定出多个指定位点的测序深度值；其中，每一指定位点对应有置信区间及rd设定均值；
[0114]
异常检测单元602，用于在任一指定位点的测序深度值未位于该指定位点对应的置信区间内时，确定该指定位点为拷贝异常位点；
[0115]
计算单元603，用于根据各个拷贝异常位点对应的rd设定均值以及测序深度值，计算各个拷贝异常位点的拷贝数值；其中，拷贝数值与测序深度值呈正相关关系；
[0116]
聚类单元604，用于将所有拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类；
[0117]
变异确定单元605，用于将归属于高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点；
[0118]
合并单元606，用于将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。
[0119]
可选的，图6所示的基因拷贝数变异检测装置还可以包括以下未图示的单元：
[0120]
数据获取单元，用于获取多个训练样本的测序数据；测序数据包括m1个候选位点的测序深度值；
[0121]
分类单元，用于对所有候选位点的测序深度值进行无监督聚类分类，获得多个分类类别；
[0122]
求均单元，用于计算每一分类类别的第一rd均值；
[0123]
噪音识别单元，用于根据第一rd均值从多个分类类别中识别出噪音类别；
[0124]
噪音剔除单元，用于将噪音类别所包括的候选位点剔除，获得m2个设定位点；
[0125]
基线计算单元，用于计算每一设定位点在多个训练样本下的rd设定均值及方差；以及，从m2个设定位点中确定出全部或部分设定位点作为指定位点；以及，根据rd设定均值及方差，计算每一指定位点的置信区间。
[0126]
如图7所示，本发明实施例公开一种电子设备，包括存储有可执行程序代码的存储器701以及与存储器701耦合的处理器702；
[0127]
其中，处理器702调用存储器701中存储的可执行程序代码，执行上述各实施例中描述的基因拷贝数变异检测方法。
[0128]
本发明实施例还公开一种计算机可读存储介质，其存储计算机程序，其中，该计算机程序使得计算机执行上述各实施例中描述的基因拷贝数变异检测方法。
[0129]
以上实施例的目的，是对本发明的技术方案进行示例性的再现与推导，并以此完整的描述本发明的技术方案、目的及效果，其目的是使公众对本发明的公开内容的理解更加透彻、全面，并不以此限定本发明的保护范围。
[0130]
以上实施例也并非是基于本发明的穷尽性列举，在此之外，还可以存在多个未列出的其他实施方式。在不违反本发明构思的基础上所作的任何替换与改进，均属本发明的保护范围。
技术特征：1.基因拷贝数变异检测方法，其特征在于，包括：从待测样本的测序数据中确定出多个指定位点的测序深度值；其中，每一所述指定位点对应有置信区间及rd设定均值；若任一指定位点的测序深度值未位于该指定位点对应的置信区间内，确定该指定位点为拷贝异常位点；根据各个所述拷贝异常位点对应的rd设定均值以及测序深度值，计算各个所述拷贝异常位点的拷贝数值；其中，所述拷贝数值与所述测序深度值呈正相关关系；将所有所述拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，所述拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类；将归属于所述高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于所述低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点；将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。2.如权利要求1所述的基因拷贝数变异检测方法，其特征在于，所述指定位点对应的置信区间及rd设定均值通过以下步骤计算：获取多个训练样本的测序数据；所述测序数据包括m1个候选位点的测序深度值；对所有候选位点的测序深度值进行无监督聚类分类，获得多个分类类别；计算每一所述分类类别的第一rd均值；根据所述第一rd均值从多个分类类别中识别出噪音类别；将所述噪音类别所包括的候选位点剔除，获得m2个设定位点；计算每一所述设定位点在多个训练样本下的rd设定均值及方差；从m2个设定位点中确定出全部或部分设定位点作为指定位点；根据所述rd设定均值及方差，计算每一所述指定位点的置信区间。3.如权利要求2所述的基因拷贝数变异检测方法，其特征在于，所述根据所述第一rd均值从多个分类类别中识别出噪音类别，包括：将所述第一rd均值小于第一指定阈值的分类类别确定为着丝粒类别；将所述第一rd均值大于第二指定阈值的分类类别确定为重复序列类别；将所述第一rd均值大于第三指定阈值的分类类别确定为线粒体类别；将所述着丝粒类别、所述重复序列类别和所述线粒体类别作为噪音类别。4.如权利要求3所述的基因拷贝数变异检测方法，其特征在于，所述计算每一所述分类类别的第一rd均值之后，所述方法还包括：根据每一所述分类类别的第一rd均值，计算多个所述分类类别的第二rd均值；将所述第二rd均值的十分之一作为所述第一指定阈值；将所述第二rd均值的五倍作为所述第二指定阈值；将所述第二rd均值的二十倍作为所述第三指定阈值。5.如权利要求2至4任一项所述的基因拷贝数变异检测方法，其特征在于，所述根据所述rd设定均值及方差，计算每一所述指定位点的置信区间，包括:每一所述指定位点的置信区间通过以下公式表示：
式中，代表第i个指定位点对应的rd设定均值，s
i2
代表第i个指定位点对应的方差。6.如权利要求1至4任一项所述的基因拷贝数变异检测方法，其特征在于，所述根据各个所述拷贝异常位点对应的rd设定均值以及测序深度值，计算各个所述拷贝异常位点的拷贝数值，包括：通过以下公式计算各个所述拷贝异常位点的拷贝数值：式中，cp
j
代表第j个拷贝异常位点的拷贝数值，rd
j
代表第j个拷贝异常位点的测序深度值，代表第j个拷贝异常位点对应的rd设定均值。7.如权利要求1至4任一项所述的基因拷贝数变异检测方法，其特征在于，所述将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并获得拷贝变异片段之后，所述方法还包括：若所述拷贝变异片段的大小达到第一片段阈值，确定所述拷贝变异片段为目标拷贝变异片段；若所述拷贝变异片段的大小达到第二片段阈值且小于所述第一片段阈值，确定所述拷贝变异片段为疑似拷贝变异片段。8.基因拷贝数变异检测装置，其特征在于，包括：深度确定单元，用于从待测样本的测序数据中确定出多个指定位点的测序深度值；其中，每一所述指定位点对应有置信区间及rd设定均值；异常检测单元，用于在任一指定位点的测序深度值未位于该指定位点对应的置信区间内时，确定该指定位点为拷贝异常位点；计算单元，用于根据各个所述拷贝异常位点对应的rd设定均值以及测序深度值，计算各个所述拷贝异常位点的拷贝数值；其中，所述拷贝数值与所述测序深度值呈正相关关系；聚类单元，用于将所有所述拷贝异常位点的拷贝数值进行空间聚类分类，获得正常拷贝类和拷贝变异类，所述拷贝变异类包括两个拷贝变异子类，分别是高拷贝变异子类和低拷贝变异子类；变异确定单元，用于将归属于所述高拷贝变异子类且测序深度值大于相应的置信区间的拷贝异常位点、以及归属于所述低拷贝变异子类且测序深度值小于相应的置信区间的拷贝异常位点，分别确定为拷贝变异位点；合并单元，用于将位置相邻且属于同一拷贝变异子类的拷贝变异位点进行合并，获得拷贝变异片段。9.电子设备，其特征在于，包括存储有可执行程序代码的存储器以及与所述存储器耦合的处理器；所述处理器调用所述存储器中存储的所述可执行程序代码，用于执行权利要求1至7任一项所述的基因拷贝数变异检测方法。10.计算机可读存储介质，其特征在于，所述计算机可读存储介质存储计算机程序，其中，所述计算机程序使得计算机执行权利要求1至7任一项所述的基因拷贝数变异检测方法。
技术总结本发明属于生物信息检测技术领域，公开了一种基因拷贝数变异检测方法，根据置信区间确定拷贝异常位点后聚类，将归属于高拷贝变异子类且RD值大于相应置信区间以及归属于低拷贝变异子类且RD值小于相应置信区间的拷贝异常位点确定为拷贝变异位点，将位置相邻且同变异类型的拷贝变异位点合并得到更精确的CNV片段，通过分离RD位点的区间大小和CNV片段的区间大小，这样可以在提高RD位点分辨率的同时保证CNV片段的稳定性，同时使用无监督机器学习方法对拷贝异常位点进行聚类分类和与置信区间比对的方法来双重识别拷贝变异位点，进一步提高拷贝变异位点的识别精确性，进而提高CNV检测的准确性、灵敏度和稳定性。灵敏度和稳定性。灵敏度和稳定性。
技术研发人员：赵哲 韩雪莹 刘玉霞 李敏 马楠
受保护的技术使用者：郑州金域临床检验中心有限公司
技术研发日：2022.07.20
技术公布日：2022/11/1
