确定核酸样本中CNV的方法、装置及应用与流程

确定核酸样本中cnv的方法、装置及应用
技术领域
1.本发明涉及基因测序领域，具体涉及一种确定核酸样本中cnv的方法、装置及应用。


背景技术：

2.拷贝数变异(copy number variations,cnvs)是一种常见的基因组结构变异，和人类表型多态性、疾病易感性、药物反应程度及群体遗传学密切相关。传统上，cnv主要是通过基于微阵列的技术，例如arraycgh或者snparrays。然而，由于这些技术主要基于探针的密度或者已知的标记，因此在cnv的大小和分辨率上有很大的局限性。而随着二代技术的广泛应用，cnv的检测精度有了很大的提高。虽然目前已有大量的基于二代测序的cnv检测软件被开发出来，但普遍存在假阳性过高的问题。
3.为了减少假阳性的cnv，现有技术中常用的有三种策略：1)基于保守的比对策略：某些软件在使用时采用保守的比对策略，即仅仅考虑唯一比对的读段用于分析，模糊匹配的读段被过滤掉；2)基于经验阈值的过滤策略：对于检测到的cnv，定义一个过滤的阈值，例如其些软件会要求检测出来的cnv需要两个以上读段的支持或者其他的一些比对信息作为cnvs的证据；3)基于多软件合并的策略：采用多个软件同时检测cnv，最终以同时检测出来的cnv作为最终的cnv。
4.然而这些cnv评估策略未能得到理想的评估和过滤效果，针对样本中cnv的检测还需要进一步改进。


技术实现要素：

5.本发明旨在至少在一定程度上解决相关技术中的技术问题之一。
6.本发明的发明人在研究中发现：常用的针对cnv的三种评估策略的检测结果并不理想，表现为：会引入新的假阳性cnv和人工误差，步骤繁琐，且会遗漏大量真实的cnv，未能达到理想的评估和过滤效果。例如，基于保守的比对策略，在处理时，仅考虑唯一比对的读段，未比对上的读段也采用无错配比对的原则，而这会导致新的假阳性cnv被检出。基于经验阈值的过滤策略常常需要定义一个经验上的阈值，这种阈值定义很容易引入主观偏见，而且通常来说这种阈值是很难确定的。基于多软件合并的策略，需要多个软件同时进行检测，步骤繁琐，而且由于各软件的一致性较低，会导致大量的真实cnv被遗漏掉。
7.为了提高cnv检测的准确性，本发明的发明人提供了一种确定核酸样本中cnv的方法、确定核酸样本中cnv的装置、计算机可读介质和电子设备。所提供的方法和装置通过引入机器学习模型并结合能够表征不同cnv的分类特征，对初始cnv进行筛选，能够很好的排除人工筛选所引入的误差。而且只需要使用单个软件检测即可，最大程度保留了真实的cnv。
8.为此，本发明的一个目的在于提供一种确定核酸样本中cnv的方法、装置及应用。
9.具体而言，本发明提供了如下技术方案：
10.在本发明的第一方面，本发明提供了一种确定核酸样本中cnv的方法，包括：(1)获取所述核酸样本的测序数据；(2)基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；(3)针对每个所述初始cnv，确定分类特征；(4)基于每个所述初始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右侧两个测序读段长度区域的非配对测序读段数目。
11.这些分类特征经过证实，能够用于区分不同的cnv。基于这些分类特征，并利用预先构建的机器学习模型，能够实现初始cnv集合的筛选，并获得最终的cnv集合。在筛选时，所应用的分类特征越多，所获得的最终cnv集合的精准率也越高，所述分类特征包括上述所列出的16个分类特征中的至少1个，至少3个，至少5个，至少8个，至少10个，至少12个，至少14个。更优选包括上述所列出的16个分类特征的全部。由此所提供的确定核酸样本中cnv的方法，所获得的cnv结果更加真实可靠，且由于引入了机器学习模型进行筛选，能够很好的排除人工筛选所引入的误差。
12.根据本发明的实施例，以上所述确定核酸样本中cnv的方法可以进一步包括如下技术特征：
13.根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。例如，可以通过所获得的原始的特征值除以覆盖度获得。
14.根据本发明的实施例，所述机器模型为随机森林模型。采用该随机森林模型，能够获得高准确率的cnv检测结果。
15.根据本发明的实施例，所述cnv检测软件包括选自manta、cnvnator、lumpy、delly和breakdancer的至少之一。
16.根据本发明的实施例，在进行步骤(4)之前，预先对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
17.(a)除去与基因组序列gap区重叠的cnv；
18.(b)除去长度小于100bp的cnv；
19.(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
20.(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
21.在本发明的第二方面，本发明提供了一种确定核酸样本中cnv的装置，包括：获取
模块，用于获取所述核酸样本的测序数据；初始cnv检测模块，用于基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；分类特征确定模块，用于针对每个所述初始cnv，确定分类特征；筛选模块，用于基于每个所述初始cnv的所述分类特征，利用预先构建的机器学习模型，对所述初始cnv集合进行筛选，以便获得最终cnv集合，
22.所述分类特征包括选自下列的至少之一：
23.所述初始cnv的测序深度、
24.所述初始cnv左侧两个测序读段长度区域的测序深度、
25.所述初始cnv右侧两个测序读段长度区域的测序深度、
26.所述初始cnv的gc含量、
27.所述初始cnv左侧两个测序读段长度区域的gc含量、
28.所述初始cnv右侧两个测序读段长度区域的gc含量、
29.所述初始cnv的类型、
30.所述初始cnv的长度、
31.所述初始cnv左侧两个测序读段长度区域中质量值<10的测序读段数目、
32.所述初始cnv右侧两个测序读段长度区域中质量值<10的测序读段数目、
33.所述初始cnv左侧两个测序读段长度区域的分裂比对读段的数目、
34.所述初始cnv右侧两个测序读段长度区域的分裂比对读段的数目、
35.所述初始cnv左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
36.所述初始cnv右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
37.所述初始cnv左侧两个测序读段长度区域的非配对测序读段数目、和
38.所述初始cnv右侧两个测序读段长度区域的非配对测序读段数目。
39.根据本发明的实施例，以上所述确定核酸样本中cnv的装置可以进一步包括如下技术特征：
40.根据本发明的实施例，所述装置进一步包括：
41.标准化单元，用于对所述分类特征进行标准化处理，所述标准化处理采用覆盖度进行。
42.根据本发明的实施例，所述装置中所述机器模型为随机森林模型。
43.根据本发明的实施例，所述装置中所述cnv检测软件包括选自manta、cnvnator、lumpy、delly和breakdancer的至少之一。
44.根据本发明的实施例，所述装置进一步包括：
45.过滤合并单元，用于对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
46.(a)除去与基因组序列gap区重叠的cnv；
47.(b)除去长度小于100bp的cnv；
48.(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
49.(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
50.在本发明的第三方面，本发明提供了一种构建机器学习分类模型的方法，包括：
51.(1)基于已知cnv信息和正常基因组区域的样本，建立训练集和测试集；
52.(2)获取所述训练集和所述测试集的每一个样本的测序数据；
53.(3)基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；
54.(3)针对每个所述初始cnv，确定分类特征；
55.(4)基于所述初始cnv与所述已知cnv信息和正常基因组区域的比对，对所述初始cnv在真实cnv和假阳性cnv之间进行分类，其中，将真实cnv赋值1，假阳性cnv赋值0；和
56.(5)利用步骤(3)中获得的所述分类特征作为输入特征，步骤(4)中得到的所述分类结果作为输出结果，采用所述训练集和测试集，进行机器学习训练，以便构建用于确定核酸样本中cnv的机器学习模型，
57.其中，
58.所述分类特征包括选自下列的至少之一：
59.所述初始cnv的测序深度、
60.所述初始cnv左侧两个测序读段长度区域的测序深度、
61.所述初始cnv右侧两个测序读段长度区域的测序深度、
62.所述初始cnv的gc含量、
63.所述初始cnv左侧两个测序读段长度区域的gc含量、
64.所述初始cnv右侧两个测序读段长度区域的gc含量、
65.所述初始cnv的类型、
66.所述初始cnv的长度、
67.所述初始cnv左侧两个测序读段长度区域中质量值<10的测序读段数目、
68.所述初始cnv右侧两个测序读段长度区域中质量值<10的测序读段数目、
69.所述初始cnv左侧两个测序读段长度区域的分裂比对读段的数目、
70.所述初始cnv右侧两个测序读段长度区域的分裂比对读段的数目、
71.所述初始cnv左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
72.所述初始cnv右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
73.所述初始cnv左侧两个测序读段长度区域的非配对测序读段数目、和
74.所述初始cnv右侧两个测序读段长度区域的非配对测序读段数目。
75.根据本发明的实施例，以上所述构建机器学习分类模型的方法可以进一步包括如下技术特征：
76.根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。
77.根据本发明的实施例，所述机器学习模型为随机森林模型。
78.根据本发明的实施例，所述cnv检测软件包括选自manta、cnvnator、lumpy、delly和breakdancer的至少之一。
79.根据本发明的实施例，预先对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
80.(a)除去与基因组序列gap区重叠的cnv；
81.(b)除去长度小于100bp的cnv；
82.(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
83.(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
84.在本发明的第四方面，本发明提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现本发明第一方面或者本发明第三方面中任一项所述方法的步骤。
85.在本发明的第五方面，本发明提供了一种电子设备，包括：本发明第四方面所述的计算机可读存储介质；以及一个或者多个处理器，用于执行所述计算机可读存储介质中的程序。
86.本发明所取得的有益效果为：本发明所提供的确定核酸样本中cnv的方法和确定核酸样本中cnv的装置，其在获得初始cnv集合后，又基于机器学习模型对初始cnv集合中的初始cnv进行了筛选。这些机器学习模型能够基于不同的分类特征实现cnv的分类表征，由此经过筛选的cnv其结果更加真实可靠，不会引入新的假阳性cnv和人工误差，只需要使用单个软件检测即可，最大程度保留了真实的cnv。
附图说明
87.图1是根据本发明的实施例提供的构建机器学习模型的方法的示意图。
88.图2是根据本发明的实施例提供的确定核酸样本中cnv的装置的结构示意图。
89.图3是根据本发明的实施例提供的采用不同方法获得的真阳性cnv和假阳性cnv的个数。
90.图4是根据本发明的实施例提供的采用不同方法获得的cnv的精准度和召回率结果。
91.图5是根据本发明的实施例提供的采用不同预测模型所获得的cnv的精准度和召回率结果。
92.图6是根据本发明的实施例提供的采用不同分类特征组合效果图。
具体实施方式
93.下面参考附图详细描述本发明的实施例，需要说明的是，所描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
94.拷贝数异常(copy number variations,cnvs)是属于基因组结构变异(structural variation)，根据大小可分为两个层次：显微水平(microscopic)和亚显微水平(submicroscopic)。本文中所提到的cnv既包括显微水平上的变异，也包括亚显微水平上的变异。所提到的显微水平的基因组结构变异主要是指显微镜下可见的染色体畸变，包括整倍体或非整倍体、缺失、插入、倒位、易位、脆性位点等结构变异。所提到的亚显微水平的基因组结构变异是指dna片段长度在1kb-3mb的基因组结构变异，包括缺失、插入、重复、重排、倒位、dna拷贝数目变化等。
95.为了获得准确的cnv检测结果，发明人创造性地想到将机器学习应用于cnv的检测中。那么如何构建合适的机器学习模型，对于cnv检测结果的准确性至关重要。
96.参照图1所示，为了构建机器学习模型，首先需要对测序数据进行质控和过滤，获得合格数据，然后选取合适的cnv检测软件，进行cnv检测，初始cnv集合。同时需要有已知的金标准cnv样本。例如可以使用千人基因组的几个标准细胞系样本，这些样本具备全面的已经发表的cnv金标准集(参考文献mark j.p.chaisson et al.multi-platform discovery of haplotype-resolved structural variation in human genomes.nature communication,2019)。然后根据已知的金标准cnv集和正常基因组区域，判断软件检测出来的原始cnv是真实的还是假阳性的，若为真实的则标记为1，若为假阳性的则标记为0。采用的判断的标准为：初始cnv集合与已知的金标准cnv集有50％以上的重叠率，即标记为1。这样，就得到了具有1/0标签的cnv集合(即图1中标记为真或假的cnv集)。当然，所采用的判断标准也可以更加严格，例如可以选择有70％以上的重叠率，才标记为1，剩余的即为0。然后从原始比对的bam文件中提取所有1/0标签的cnv的分类特征信息，得到用于训练的信号矩阵。选取适当的机器学习算法，对提取的特征进行训练，通过参数调整，不断优化已有的模型。从而构建得到机器学习模型。
97.在本发明的一些优选实施方式中，所选取的分类特征分别如下表1所示。
98.表1分类特征信息
[0099][0100][0101]
其中，表1中所提到的cnv类型是指缺失(deletion)和重复(duplication)，所提到的分裂比对读段是指将测序读段比对到基因组的时候，由于测序读段只有一部分比对到基因组上，所以将该测序读段称为分裂比对读段，其数目可以通过下述方法或者公式计算获
得：
[0102]
将测序数据进行基因组比对得到比对结果bam文件，根据bam文件的第6列(cigar)的比对标志信息“s”来判断该读段是否为分裂比对读段，若有“s”标志则为分裂比对读段，计数加1。使用pysam包实现统计某一比对区域内所有分裂比对读段的数目。
[0103]
所提到的插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外(可以用英文insertsize>mean /-3sd表示)的测序读段的数目，可以通过下述公式或者方法计算获得：
[0104]
将测序数据进行基因组比对得到比对结果bam文件，使用collectinsertsizemetrics.jar函数统计bam文件平均的插入片段大小和标准差；使用配对的读段基因组比对起始位置相减得到该配对读段的插入片段大小；若该配对读段的插入片段大小超过了平均的插入片段大小左右三倍标准差的范围，则计数。
[0105]
采用pe进行测序，所获得的读段是配对的(read1和read2)，所以除去配对的测序读段即为非配对测序读段，非配对测序读段数目可以通过下述方法计算获得：
[0106]
使用pysam包的is_pair()函数判断某个读段是否为配对的读段，若为非配对的读段，则计数。
[0107]
本发明使用已有的真集样本训练得到最终的分类模型，类似于现有的基于经验阈值基于的过滤策略，不同的是现有的方法采用的阈值往往是人为根据经验主观定义的，这种定义方式往往会引入人工的误差，同时这种方法往往只考虑单一或很少因素来筛选cnv，如在候选的cnv断点处，要求两个以上测序读长支持，这种筛选方式是片面的。而本发明通过比对，最终确定了cnv的测序深度、gc含量、读段剪切、插入片段的大小等16个因素，使得结果更准确可靠，机器学习的训练也很好的排除了人工引入的误差。
[0108]
为此，在本发明的一个方面，本发明提供了一种构建机器学习分类模型的方法，包括：
[0109]
(1)基于已知cnv信息和正常基因组区域的样本，建立训练集和测试集；
[0110]
(2)获取所述训练集和所述测试集的每一个样本的测序数据；
[0111]
(3)基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；
[0112]
(3)针对每个所述初始cnv，确定分类特征；
[0113]
(4)基于所述初始cnv与所述已知cnv信息和正常基因组区域的比对，对所述初始cnv在真实cnv和假阳性cnv之间进行分类，其中，将真实cnv赋值1，假阳性cnv赋值0；和
[0114]
(5)利用步骤(3)中获得的所述分类特征作为输入特征，步骤(4)中得到的所述分类结果作为输出结果，采用所述训练集和测试集，进行机器学习训练，以便构建用于确定核酸样本中cnv的机器学习模型，
[0115]
其中，
[0116]
所述分类特征包括选自下列的至少之一：
[0117]
所述初始cnv的测序深度、
[0118]
所述初始cnv左侧两个测序读段长度区域的测序深度、
[0119]
所述初始cnv右侧两个测序读段长度区域的测序深度、
[0120]
所述初始cnv的gc含量、
[0121]
所述初始cnv左侧两个测序读段长度区域的gc含量、
[0122]
所述初始cnv右侧两个测序读段长度区域的gc含量、
[0123]
所述初始cnv的类型、
[0124]
所述初始cnv的长度、
[0125]
所述初始cnv左侧两个测序读段长度区域中质量值<10的测序读段数目、
[0126]
所述初始cnv右侧两个测序读段长度区域中质量值<10的测序读段数目、
[0127]
所述初始cnv左侧两个测序读段长度区域的分裂比对读段的数目、
[0128]
所述初始cnv右侧两个测序读段长度区域的分裂比对读段的数目、
[0129]
所述初始cnv左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0130]
所述初始cnv右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0131]
所述初始cnv左侧两个测序读段长度区域的非配对测序读段数目、和
[0132]
所述初始cnv右侧两个测序读段长度区域的非配对测序读段数目。
[0133]
根据本发明的实施例，以上所述构建机器学习分类模型的方法可以进一步包括如下技术特征：
[0134]
根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。
[0135]
根据本发明的实施例，所述机器学习模型为随机森林模型。随机森林模型具有良好的预测性能。
[0136]
根据本发明的实施例，所述cnv检测软件可以是本领域检测cnv的常用软件，包括但不限于manta、cnvnator、lumpy、delly和breakdancer等。
[0137]
根据本发明的实施例，预先对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
[0138]
(a)除去与基因组序列gap区重叠的cnv；
[0139]
(b)除去长度小于100bp的cnv；
[0140]
(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
[0141]
(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
[0142]
所提到的gap区也称为n区，可以使用ucsc数据库下载基因组的n区的基因组坐标文件，使用自定义脚本查看cnv的坐标与n区的基因组坐标是否有重叠的地方，若重叠区域>20bp则去除该cnv，以此除去与基因组序列gap区重叠的cnv。通过过滤和合并处理，可以过滤去除一些基因组未知区域，或者软件检测局限所导致的假阳性cnv，提高cnv的检测准确性。
[0143]
在本发明的另一个方面，本发明提供了一种确定核酸样本中cnv的方法，包括：(1)获取所述核酸样本的测序数据；(2)基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；(3)针对每个所述初始cnv，确定分类特征；(4)基于每个所述初始cnv的所述分类特征，利用预先构建的机器学习模型，对所述初始cnv集合进行筛选，以便获得最终cnv集合，
[0144]
所述分类特征包括选自下列的至少之一：
[0145]
所述初始cnv的测序深度、
[0146]
所述初始cnv左侧两个测序读段长度区域的测序深度、
[0147]
所述初始cnv右侧两个测序读段长度区域的测序深度、
[0148]
所述初始cnv的gc含量、
[0149]
所述初始cnv左侧两个测序读段长度区域的gc含量、
[0150]
所述初始cnv右侧两个测序读段长度区域的gc含量、
[0151]
所述初始cnv的类型、
[0152]
所述初始cnv的长度、
[0153]
所述初始cnv左侧两个测序读段长度区域中质量值<10的测序读段数目、
[0154]
所述初始cnv右侧两个测序读段长度区域中质量值<10的测序读段数目、
[0155]
所述初始cnv左侧两个测序读段长度区域的分裂比对读段的数目、
[0156]
所述初始cnv右侧两个测序读段长度区域的分裂比对读段的数目、
[0157]
所述初始cnv左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0158]
所述初始cnv右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0159]
所述初始cnv左侧两个测序读段长度区域的非配对测序读段数目、和
[0160]
所述初始cnv右侧两个测序读段长度区域的非配对测序读段数目。
[0161]
本发明所提供的确定核酸样本中cnv的方法，使用单个软件训练模型，不需要多个软件同时检测，步骤简单。所提到的核酸样本可以是任何所需要检测的核酸样本。这些核酸样本可以通过本领域常用的一些方法获得，例如借助于试剂盒进行基因组dna的提取等等。
[0162]
根据本发明的实施例，以上所述确定核酸样本中cnv的方法可以进一步包括如下技术特征：
[0163]
根据本发明的实施例，所述分类特征是经过标准化处理的，所述标准化处理采用覆盖度进行。
[0164]
根据本发明的实施例，所述机器模型为随机森林模型。
[0165]
根据本发明的实施例，所述cnv检测软件包括选自manta、cnvnator、lumpy、delly和breakdancer的至少之一。
[0166]
根据本发明的实施例，在进行步骤(4)之前，预先对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
[0167]
(a)除去与基因组序列gap区重叠的cnv；
[0168]
(b)除去长度小于100bp的cnv；
[0169]
(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
[0170]
(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
[0171]
在本发明的又一个方面，本发明提供了一种确定核酸样本中cnv的装置，如图2所示，包括：获取模块，用于获取所述核酸样本的测序数据；初始cnv检测模块，用于基于所述测序数据，利用cnv检测软件，确定由多个初始cnv构成的初始cnv集合；分类特征确定模块，用于针对每个所述初始cnv，确定分类特征；筛选模块，用于基于每个所述初始cnv的所述分类特征，利用预先构建的机器学习模型，对所述初始cnv集合进行筛选，以便获得最终cnv集
合，
[0172]
所述分类特征包括选自下列的至少之一：
[0173]
所述初始cnv的测序深度、
[0174]
所述初始cnv左侧两个测序读段长度区域的测序深度、
[0175]
所述初始cnv右侧两个测序读段长度区域的测序深度、
[0176]
所述初始cnv的gc含量、
[0177]
所述初始cnv左侧两个测序读段长度区域的gc含量、
[0178]
所述初始cnv右侧两个测序读段长度区域的gc含量、
[0179]
所述初始cnv的类型、
[0180]
所述初始cnv的长度、
[0181]
所述初始cnv左侧两个测序读段长度区域中质量值<10的测序读段数目、
[0182]
所述初始cnv右侧两个测序读段长度区域中质量值<10的测序读段数目、
[0183]
所述初始cnv左侧两个测序读段长度区域的分裂比对读段数目、
[0184]
所述初始cnv右侧两个测序读段长度区域的分裂比对读段数目、
[0185]
所述初始cnv左侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0186]
所述初始cnv右侧两个测序读段长度区域内插入片段的大小在整个文库插入片段大小平均值的左右3倍标准差的范围以外的测序读段的数目、
[0187]
所述初始cnv左侧两个测序读段长度区域的非配对测序读段数目、和
[0188]
所述初始cnv右侧两个测序读段长度区域的非配对测序读段数目。
[0189]
根据本发明的实施例，以上所述确定核酸样本中cnv的装置可以进一步包括如下技术特征：
[0190]
根据本发明的实施例，所述装置进一步包括：
[0191]
标准化单元，用于对所述分类特征进行标准化处理，所述标准化处理采用覆盖度进行。
[0192]
根据本发明的实施例，所述装置中所述机器模型为随机森林模型。
[0193]
根据本发明的实施例，所述装置中所述cnv检测软件包括选自manta、cnvnator、lumpy、delly和breakdancer的至少之一。
[0194]
根据本发明的实施例，所述装置进一步包括：
[0195]
过滤合并单元，用于对所述初始cnv集合进行过滤和合并处理，所述过滤和合并处理包括：
[0196]
(a)除去与基因组序列gap区重叠的cnv；
[0197]
(b)除去长度小于100bp的cnv；
[0198]
(c)确定在不同软件中均存在的重复cnv，并对所述重复cnv进行合并；和
[0199]
(d)如果所述重复cnv在不同软件中体现的类型不同，则删除所述重复cnv。
[0200]
下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。
[0201]
实施例1
[0202]
实施例1提供了一个对现有cnv检测软件的检测结果进行评估的分类模型，该模型是通过对6个已知真实cnv的全基因组测序数据进行训练得到的，并且在另外的3个样本中得到验证(实现了90％以上的精准度和查全率)。基于该模型可对未知cnv样本的cnv检测结果进行评估和过滤，以实现较高的准确度，为临床检测cnv提供很好的助力。具体包括如下步骤：
[0203]
(1)测序数据获取、质控和过滤
[0204]
从ncbi下载9个千人基因组标准细胞系的测序数据，其中六个样本(srr7782690,srr7782691，srr7782669，srr7782675，srr7782676，srr7782677)作为模型训练集样本，另外的三个样本(srr7782672，srr7782673，srr7782674)作为验证集样本。
[0205]
使用fastqc(v0.11.5)对下载的测序数据进行质控，使用soapnuke(v1.5.0)进行低质量读段的过滤(其中质量值低于15的碱基数大于读段(reads)总长的50％的读段被过滤掉；未知碱基(n)数占到读段总碱基数的10％以上的读段被过滤掉)，得到最终合格的数据用于后续分析。
[0206]
(2)比对，去重
[0207]
使用bwa-men将测序数据比对到人类参考基因组上(hg19)，然后使用picard(v 2.18.11)去除pcr重复序列，即去除由于pcr扩增导致的重复序列。获得bam文件。
[0208]
(3)cnv检测
[0209]
使用cnv检测软件(manta/lumpy/delly)对比对后的bam文件进行cnv的检测得到初始cnv集合。
[0210]
(4)cnv初步过滤合并：
[0211]
a)移除与基因组gap区(n区)重叠的cnvs：使用ucsc数据库下载基因组的n区的基因组坐标文件，使用自定义脚本查看cnv的坐标与n区的基因组坐标是否有重叠的地方，若重叠区域>20bp则去除该cnv，以此除去与基因组序列gap区重叠的cnv；
[0212]
b)移除小于100bp的cnvs；
[0213]
c)合并各软件内部重叠的cnvs；
[0214]
d)移除软件内部矛盾的cnvs(两个cnv重叠，但cnv的类型不同)。
[0215]
基于以上标准得到原始的候选cnv集合。
[0216]
(5)标记cnv：
[0217]
将原始的候选cnv集中与已知真集的cnv区域(参考文献mark j.p.chaisson et al.multi-platform discovery of haplotype-resolved structural variation in human genomes.nature communication,2019，详见：dbvar:nstd152)重叠率>50％标记为1；否则标记为0，得到具有1/0标签的cnv集合。
[0218]
(6)特征提取：
[0219]
使用自定义的python脚本(get_feature.py)提取每个cnv的比对信息，形成特征矩阵。具体来说：对(5)中得到的带有1/0标签的cnv集合中的每个cnvs分别提取对应的比对信息(详细信息见上表1)，其中测序深度使用log2标准化，低比对质量(q<10)的read数目、软剪切的read数目、异常的配对的read数目、非配对的read数目都使用整个样本平均的测序深度进行均一化。该信息的提取主要使用python包pysam对比对后的bam文件进行处理得
到，最终每一个cnv都有一行对应的特征值列表，合并即得到所有cnvs特征值的一个二维矩阵。
[0220]
(7)模型训练：
[0221]
使用随机森林的算法对特征矩阵进行训练，得到最终用于cnv分类的模型。该模型训练主要是构建(6)中提取的特征值(可看作x)与已知的类别标签(可看作y)的关系函数。因此，当得到最终的函数/模型时，可根据提供的新样本的特征值(x)预测最终的类别标签(y)。
[0222]
(8)模型预测性能验证：
[0223]
利用另外三个无关的样本对该分类模型进行验证，对这三个样本同样进行以上(1)(2)(3)(4)(6)，最终将(6)中得到的特征矩阵输入模型，得到每个cnv的评估和预测结果。
[0224]
结果分析如下：
[0225]
1、与现有技术预测结果进行比较：
[0226]
(1)基于经验阈值的过滤策略：
[0227]
由于该策略不唯一，我们暂且使用两个较常用的标准：
[0228]
1)split reads数目>2或discordant read pair数目>2，使用“sr>2||pe>2”表示；
[0229]
2)split reads数目>＝2且discordant read pair数目>＝2，使用“sr>＝2&&pe>＝2”表示。
[0230]
(2)基于多软件合并的策略：
[0231]
这里同时使用delly和lumpy进行检测，并取其共同检测到的cnvs。使用“delly&&lumpy”表示。
[0232]
初始的检测cnvs使用“row cnvs”和“row(delly lumpy)”表示。
[0233]
本发明方法预测结果使用“pre cnvs”和“pre(delly lumpy)”表示。
[0234]
具体结果如图3所示。
[0235]
可以看出使用本方法极大的降低了假阳性cnv的数目(图3中阴影部分所示)，而最大程度保留了真实的cnv(图3中非阴影部分所示)，显著优于现有的技术。
[0236]
2、本发明方法和已有方法预测的具体精准度和召回率如下：
[0237]
精准度(precise)＝pre_fp/(pre_fp pre_fp)
[0238]
召回率(recall)＝pre_tp/row_tp
[0239]
row_tp:未使用该模型之前原始软件检测到的真阳性cnvs的个数；
[0240]
row_fp:未使用该模型之前原始软件检测到的假阳性cnvs的个数；
[0241]
pre_tp:使用该模型之后软件检测到的真阳性cnvs的个数；
[0242]
pre_fp:使用该模型之后软件检测到的假阳性cnvs的个数。
[0243]
其结果如图4所示。其结果表明：使用机器学习模型具有很好的精准率和召回率，并且优于现有方法。
[0244]
另外，发明人同时研究了在进行预测时，采用不同的模型对于预测结果的影响，即在上述步骤(7)中除了随机森林(random forest，缩写为rf)模型之外，还利用梯度提升数(gradient boosting classifier，缩写为gbc)和支持向量机(support vector machine，缩写为svm)两种预测模型进行了预测。三种模型的预测结果如图5所示。
[0245]
从图5不难看出，随机森林模型相较于梯度提升数和支持向量机模型来说，无论是精准率还是召回率，更优，表现出更优异的效果。
[0246]
同时，发明人研究了在进行预测时，分类特征的多少对于预测精准度和召回率的影响。图6为分类特征组合效果图，其中1-feature是只选用1个分类特征(read_depth)的预测精准度和召回率；3-features是选用3个分类特征(read_depth，read_depth_l，read_depth_r)的预测精准度和召回率；all-features是采用所有的分类特征的预测精准度和召回率。可以看出即使选用少量分类特征也有一定的预测效果，例如，选择其中一个分类特征进行预测，其精准度在88％以上，召回率在88％以上。进一步地，当分类特征越多时，结果表现越好。
[0247]
在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0248]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
[0249]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。