具体实施方式

下面参考附图详细描述本发明的实施例，需要说明的是，所描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

拷贝数异常(copy number variations,CNVs)是属于基因组结构变异(structural variation)，根据大小可分为两个层次：显微水平(microscopic)和亚显微水平(submicroscopic)。本文中所提到的CNV既包括显微水平上的变异，也包括亚显微水平上的变异。所提到的显微水平的基因组结构变异主要是指显微镜下可见的染色体畸变，包括整倍体或非整倍体、缺失、插入、倒位、易位、脆性位点等结构变异。所提到的亚显微水平的基因组结构变异是指DNA片段长度在1Kb-3Mb的基因组结构变异，包括缺失、插入、重复、重排、倒位、DNA拷贝数目变化等。

为了获得准确的CNV检测结果，发明人创造性地想到将机器学习应用于CNV的检测中。那么如何构建合适的机器学习模型，对于CNV检测结果的准确性至关重要。

参照图1所示，为了构建机器学习模型，首先需要对测序数据进行质控和过滤，获得合格数据，然后选取合适的CNV检测软件，进行CNV检测，初始CNV集合。同时需要有已知的金标准CNV样本。例如可以使用千人基因组的几个标准细胞系样本，这些样本具备全面的已经发表的CNV金标准集(参考文献Mark J.P.Chaisson et al.Multi-platform discoveryof haplotype-resolved structural variation in human genomes.NatureCommunication,2019)。然后根据已知的金标准CNV集和正常基因组区域，判断软件检测出来的原始CNV是真实的还是假阳性的，若为真实的则标记为1，若为假阳性的则标记为0。采用的判断的标准为：初始CNV集合与已知的金标准CNV集有50％以上的重叠率，即标记为1。这样，就得到了具有1/0标签的CNV集合(即图1中标记为真或假的CNV集)。当然，所采用的判断标准也可以更加严格，例如可以选择有70％以上的重叠率，才标记为1，剩余的即为0。然后从原始比对的BAM文件中提取所有1/0标签的CNV的分类特征信息，得到用于训练的信号矩阵。选取适当的机器学习算法，对提取的特征进行训练，通过参数调整，不断优化已有的模型。从而构建得到机器学习模型。

在本发明的一些优选实施方式中，所选取的分类特征分别如下表1所示。

表1分类特征信息

其中，表1中所提到的CNV类型是指缺失(deletion)和重复(duplication)，所提到的分裂比对读段是指将测序读段比对到基因组的时候，由于测序读段只有一部分比对到基因组上，所以将该测序读段称为分裂比对读段，其数目可以通过下述方法或者公式计算获得：

将测序数据进行基因组比对得到比对结果BAM文件，根据BAM文件的第6列(CIGAR)的比对标志信息“S”来判断该读段是否为分裂比对读段，若有“S”标志则为分裂比对读段，计数加1。使用pysam包实现统计某一比对区域内所有分裂比对读段的数目。

所提到的插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外(可以用英文insertsize>mean+/-3sd表示)的测序读段的数目，可以通过下述公式或者方法计算获得：

将测序数据进行基因组比对得到比对结果BAM文件，使用CollectInsertSizeMetrics.jar函数统计bam文件平均的插入片段大小和标准差；使用配对的读段基因组比对起始位置相减得到该配对读段的插入片段大小；若该配对读段的插入片段大小超过了平均的插入片段大小左右三倍标准差的范围，则计数。

采用PE进行测序，所获得的读段是配对的(read1和read2)，所以除去配对的测序读段即为非配对测序读段，非配对测序读段数目可以通过下述方法计算获得：

使用pysam包的is_pair()函数判断某个读段是否为配对的读段，若为非配对的读段，则计数。

本发明使用已有的真集样本训练得到最终的分类模型，类似于现有的基于经验阈值基于的过滤策略，不同的是现有的方法采用的阈值往往是人为根据经验主观定义的，这种定义方式往往会引入人工的误差，同时这种方法往往只考虑单一或很少因素来筛选CNV，如在候选的CNV断点处，要求两个以上测序读长支持，这种筛选方式是片面的。而本发明通过比对，最终确定了CNV的测序深度、GC含量、读段剪切、插入片段的大小等16个因素，使得结果更准确可靠，机器学习的训练也很好的排除了人工引入的误差。

为此，在本发明的一个方面，本发明提供了一种构建机器学习分类模型的方法，包括：

(1)基于已知CNV信息和正常基因组区域的样本，建立训练集和测试集；

(2)获取所述训练集和所述测试集的每一个样本的测序数据；

(3)基于所述测序数据，利用CNV检测软件，确定由多个初始CNV构成的初始CNV集合；

(3)针对每个所述初始CNV，确定分类特征；

(4)基于所述初始CNV与所述已知CNV信息和正常基因组区域的比对，对所述初始CNV在真实CNV和假阳性CNV之间进行分类，其中，将真实CNV赋值1，假阳性CNV赋值0；和

(5)利用步骤(3)中获得的所述分类特征作为输入特征，步骤(4)中得到的所述分类结果作为输出结果，采用所述训练集和测试集，进行机器学习训练，以便构建用于确定核酸样本中CNV的机器学习模型，

其中，

所述分类特征包括选自下列的至少之一：

所述初始CNV的测序深度、

所述初始CNV左侧两个测序读段长度区域的测序深度、

所述初始CNV右侧两个测序读段长度区域的测序深度、

所述初始CNV的GC含量、

所述初始CNV左侧两个测序读段长度区域的GC含量、

所述初始CNV右侧两个测序读段长度区域的GC含量、

所述初始CNV的类型、

所述初始CNV的长度、

所述初始CNV左侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV右侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV左侧两个测序读段长度区域的分裂比对读段的数目、

所述初始CNV右侧两个测序读段长度区域的分裂比对读段的数目、

所述初始CNV左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV左侧两个测序读段长度区域的非配对测序读段数目、和

所述初始CNV右侧两个测序读段长度区域的非配对测序读段数目。

根据本发明的实施例，以上所述构建机器学习分类模型的方法可以进一步包括如下技术特征：

根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。

根据本发明的实施例，所述机器学习模型为随机森林模型。随机森林模型具有良好的预测性能。

根据本发明的实施例，所述CNV检测软件可以是本领域检测CNV的常用软件，包括但不限于manta、CNVnator、lumpy、delly和breakdancer等。

根据本发明的实施例，预先对所述初始CNV集合进行过滤和合并处理，所述过滤和合并处理包括：

(a)除去与基因组序列gap区重叠的CNV；

(b)除去长度小于100bp的CNV；

(c)确定在不同软件中均存在的重复CNV，并对所述重复CNV进行合并；和

(d)如果所述重复CNV在不同软件中体现的类型不同，则删除所述重复CNV。

所提到的gap区也称为N区，可以使用UCSC数据库下载基因组的N区的基因组坐标文件，使用自定义脚本查看CNV的坐标与N区的基因组坐标是否有重叠的地方，若重叠区域>20bp则去除该CNV，以此除去与基因组序列gap区重叠的CNV。通过过滤和合并处理，可以过滤去除一些基因组未知区域，或者软件检测局限所导致的假阳性CNV，提高CNV的检测准确性。

在本发明的另一个方面，本发明提供了一种确定核酸样本中CNV的方法，包括：(1)获取所述核酸样本的测序数据；(2)基于所述测序数据，利用CNV检测软件，确定由多个初始CNV构成的初始CNV集合；(3)针对每个所述初始CNV，确定分类特征；(4)基于每个所述初始CNV的所述分类特征，利用预先构建的机器学习模型，对所述初始CNV集合进行筛选，以便获得最终CNV集合，

所述分类特征包括选自下列的至少之一：

所述初始CNV的测序深度、

所述初始CNV左侧两个测序读段长度区域的测序深度、

所述初始CNV右侧两个测序读段长度区域的测序深度、

所述初始CNV的GC含量、

所述初始CNV左侧两个测序读段长度区域的GC含量、

所述初始CNV右侧两个测序读段长度区域的GC含量、

所述初始CNV的类型、

所述初始CNV的长度、

所述初始CNV左侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV右侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV左侧两个测序读段长度区域的分裂比对读段的数目、

所述初始CNV右侧两个测序读段长度区域的分裂比对读段的数目、

所述初始CNV左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV左侧两个测序读段长度区域的非配对测序读段数目、和

所述初始CNV右侧两个测序读段长度区域的非配对测序读段数目。

本发明所提供的确定核酸样本中CNV的方法，使用单个软件训练模型，不需要多个软件同时检测，步骤简单。所提到的核酸样本可以是任何所需要检测的核酸样本。这些核酸样本可以通过本领域常用的一些方法获得，例如借助于试剂盒进行基因组DNA的提取等等。

根据本发明的实施例，以上所述确定核酸样本中CNV的方法可以进一步包括如下技术特征：

根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。

根据本发明的实施例，所述机器模型为随机森林模型。

根据本发明的实施例，所述CNV检测软件包括选自manta、CNVnator、lumpy、delly和breakdancer的至少之一。

根据本发明的实施例，在进行步骤(4)之前，预先对所述初始CNV集合进行过滤和合并处理，所述过滤和合并处理包括：

(a)除去与基因组序列gap区重叠的CNV；

(b)除去长度小于100bp的CNV；

(c)确定在不同软件中均存在的重复CNV，并对所述重复CNV进行合并；和

(d)如果所述重复CNV在不同软件中体现的类型不同，则删除所述重复CNV。

在本发明的又一个方面，本发明提供了一种确定核酸样本中CNV的装置，如图2所示，包括：获取模块，用于获取所述核酸样本的测序数据；初始CNV检测模块，用于基于所述测序数据，利用CNV检测软件，确定由多个初始CNV构成的初始CNV集合；分类特征确定模块，用于针对每个所述初始CNV，确定分类特征；筛选模块，用于基于每个所述初始CNV的所述分类特征，利用预先构建的机器学习模型，对所述初始CNV集合进行筛选，以便获得最终CNV集合，

所述分类特征包括选自下列的至少之一：

所述初始CNV的测序深度、

所述初始CNV左侧两个测序读段长度区域的测序深度、

所述初始CNV右侧两个测序读段长度区域的测序深度、

所述初始CNV的GC含量、

所述初始CNV左侧两个测序读段长度区域的GC含量、

所述初始CNV右侧两个测序读段长度区域的GC含量、

所述初始CNV的类型、

所述初始CNV的长度、

所述初始CNV左侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV右侧两个测序读段长度区域中质量值<10的测序读段数目、

所述初始CNV左侧两个测序读段长度区域的分裂比对读段数目、

所述初始CNV右侧两个测序读段长度区域的分裂比对读段数目、

所述初始CNV左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、

所述初始CNV左侧两个测序读段长度区域的非配对测序读段数目、和

所述初始CNV右侧两个测序读段长度区域的非配对测序读段数目。

根据本发明的实施例，以上所述确定核酸样本中CNV的装置可以进一步包括如下技术特征：

根据本发明的实施例，所述装置进一步包括：

标准化单元，用于对所述分类特征进行标准化处理，所述标准化处理采用覆盖度进行。

根据本发明的实施例，所述装置中所述机器模型为随机森林模型。

根据本发明的实施例，所述装置中所述CNV检测软件包括选自manta、CNVnator、lumpy、delly和breakdancer的至少之一。

根据本发明的实施例，所述装置进一步包括：

过滤合并单元，用于对所述初始CNV集合进行过滤和合并处理，所述过滤和合并处理包括：

(a)除去与基因组序列gap区重叠的CNV；

(b)除去长度小于100bp的CNV；

(c)确定在不同软件中均存在的重复CNV，并对所述重复CNV进行合并；和

(d)如果所述重复CNV在不同软件中体现的类型不同，则删除所述重复CNV。

下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。

实施例1

实施例1提供了一个对现有CNV检测软件的检测结果进行评估的分类模型，该模型是通过对6个已知真实CNV的全基因组测序数据进行训练得到的，并且在另外的3个样本中得到验证(实现了90％以上的精准度和查全率)。基于该模型可对未知CNV样本的CNV检测结果进行评估和过滤，以实现较高的准确度，为临床检测CNV提供很好的助力。具体包括如下步骤：

(1)测序数据获取、质控和过滤

从NCBI下载9个千人基因组标准细胞系的测序数据，其中六个样本(SRR7782690,SRR7782691，SRR7782669，SRR7782675，SRR7782676，SRR7782677)作为模型训练集样本，另外的三个样本(SRR7782672，SRR7782673，SRR7782674)作为验证集样本。

使用FastQC(v0.11.5)对下载的测序数据进行质控，使用Soapnuke(v1.5.0)进行低质量读段的过滤(其中质量值低于15的碱基数大于读段(reads)总长的50％的读段被过滤掉；未知碱基(N)数占到读段总碱基数的10％以上的读段被过滤掉)，得到最终合格的数据用于后续分析。

(2)比对，去重

使用BWA-men将测序数据比对到人类参考基因组上(hg19)，然后使用picard(v2.18.11)去除PCR重复序列，即去除由于PCR扩增导致的重复序列。获得bam文件。

(3)CNV检测

使用CNV检测软件(manta/Lumpy/Delly)对比对后的bam文件进行CNV的检测得到初始CNV集合。

(4)CNV初步过滤合并：

a)移除与基因组gap区(N区)重叠的CNVs：使用UCSC数据库下载基因组的N区的基因组坐标文件，使用自定义脚本查看CNV的坐标与N区的基因组坐标是否有重叠的地方，若重叠区域>20bp则去除该CNV，以此除去与基因组序列gap区重叠的CNV；

b)移除小于100bp的CNVs；

c)合并各软件内部重叠的CNVs；

d)移除软件内部矛盾的CNVs(两个CNV重叠，但CNV的类型不同)。

基于以上标准得到原始的候选CNV集合。

(5)标记CNV：

将原始的候选CNV集中与已知真集的CNV区域(参考文献Mark J.P.Chaisson etal.Multi-platform discovery of haplotype-resolved structural variation inhuman genomes.Nature Communication,2019，详见：dbVar:nstd152)重叠率>50％标记为1；否则标记为0，得到具有1/0标签的CNV集合。

(6)特征提取：

使用自定义的python脚本(Get_Feature.py)提取每个CNV的比对信息，形成特征矩阵。具体来说：对(5)中得到的带有1/0标签的CNV集合中的每个CNVs分别提取对应的比对信息(详细信息见上表1)，其中测序深度使用log2标准化，低比对质量(Q<10)的read数目、软剪切的read数目、异常的配对的read数目、非配对的read数目都使用整个样本平均的测序深度进行均一化。该信息的提取主要使用python包pysam对比对后的BAM文件进行处理得到，最终每一个CNV都有一行对应的特征值列表，合并即得到所有CNVs特征值的一个二维矩阵。

(7)模型训练：

使用随机森林的算法对特征矩阵进行训练，得到最终用于CNV分类的模型。该模型训练主要是构建(6)中提取的特征值(可看作X)与已知的类别标签(可看作Y)的关系函数。因此，当得到最终的函数/模型时，可根据提供的新样本的特征值(X)预测最终的类别标签(Y)。

(8)模型预测性能验证：

利用另外三个无关的样本对该分类模型进行验证，对这三个样本同样进行以上(1)(2)(3)(4)(6)，最终将(6)中得到的特征矩阵输入模型，得到每个CNV的评估和预测结果。

结果分析如下：

1、与现有技术预测结果进行比较：

(1)基于经验阈值的过滤策略：

由于该策略不唯一，我们暂且使用两个较常用的标准：

1)split reads数目>2或discordant read pair数目>2，使用“SR>2||PE>2”表示；

2)split reads数目>＝2且discordant read pair数目>＝2，使用“SR>＝2&&PE>＝2”表示。

(2)基于多软件合并的策略：

这里同时使用Delly和Lumpy进行检测，并取其共同检测到的CNVs。使用“Delly&&Lumpy”表示。

初始的检测CNVs使用“Row CNVs”和“Row(Delly+Lumpy)”表示。

本发明方法预测结果使用“Pre CNVs”和“Pre(Delly+Lumpy)”表示。

具体结果如图3所示。

可以看出使用本方法极大的降低了假阳性CNV的数目(图3中阴影部分所示)，而最大程度保留了真实的CNV(图3中非阴影部分所示)，显著优于现有的技术。

2、本发明方法和已有方法预测的具体精准度和召回率如下：

精准度(precise)＝PRE_FP/(PRE_FP+PRE_FP)

召回率(recall)＝PRE_TP/ROW_TP

ROW_TP:未使用该模型之前原始软件检测到的真阳性CNVs的个数；

ROW_FP:未使用该模型之前原始软件检测到的假阳性CNVs的个数；

PRE_TP:使用该模型之后软件检测到的真阳性CNVs的个数；

PRE_FP:使用该模型之后软件检测到的假阳性CNVs的个数。

其结果如图4所示。其结果表明：使用机器学习模型具有很好的精准率和召回率，并且优于现有方法。

另外，发明人同时研究了在进行预测时，采用不同的模型对于预测结果的影响，即在上述步骤(7)中除了随机森林(Random Forest，缩写为RF)模型之外，还利用梯度提升数(Gradient Boosting classifier，缩写为GBC)和支持向量机(Support Vector Machine，缩写为SVM)两种预测模型进行了预测。三种模型的预测结果如图5所示。

从图5不难看出，随机森林模型相较于梯度提升数和支持向量机模型来说，无论是精准率还是召回率，更优，表现出更优异的效果。

同时，发明人研究了在进行预测时，分类特征的多少对于预测精准度和召回率的影响。图6为分类特征组合效果图，其中1-feature是只选用1个分类特征(read_depth)的预测精准度和召回率；3-features是选用3个分类特征(read_depth，read_depth_l，read_depth_r)的预测精准度和召回率；All-features是采用所有的分类特征的预测精准度和召回率。可以看出即使选用少量分类特征也有一定的预测效果，例如，选择其中一个分类特征进行预测，其精准度在88％以上，召回率在88％以上。进一步地，当分类特征越多时，结果表现越好。

在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。