基于群体样本深度信息检测基因拷贝数变异的方法及系统与流程

1.本发明涉及生物技术领域，更具体而言涉及一种基于群体样本深度信息检测基因拷贝数变异的方法及系统。


背景技术：

2.基因拷贝数变异(cnv)检出一般基于已有的几种测序信息：1，基于序列组装信息，将组装序列和参考基因组序列比较；2，基于测序深度信息，将目标区域的测序深度与周围区域及对照人群该区域对比看是否存在差异；3，基于测序读长信息，基于双端测序(pe)读长比对结果的插入片段长度以及读长可以剪切开比对到不同位置的情况。现有技术文章《pscc:sensitiveand reliable population-scale copy number variation detection method basedon low coverage sequencing》中提出了cnv检出方法(简称pscc方法)，该方法基于测序深度计算cnv。
3.然而，现有技术的方法只能检出罕见cnv，无法检出常见频率的cnv，原因在于测序深度的方法检出cnv的算法为了保证准确性会把群体样本作为正常对照即当成没有变异的情况，这样虽然提高了准确性但是牺牲了敏感性，尤其是常见频率cnv；或者检出常见频率cnv但是准确性差；现有技术的方法不能准确性针对群体遗传性cnv正确分型；而且，现有技术的方法针对拷贝数扩增(dup)不能细化具体拷贝数目。


技术实现要素：

4.针对现有技术中存在的问题，本发明提出了一种检测基因拷贝数变异的方法及系统，有效解决常见频率的cnv检测、群体遗传性cnv分型和细化dup拷贝数的问题。
5.因此，在一方面，本发明提供了一种基于群体样本深度信息检测基因拷贝数变异的方法，所述方法包括：
6.(1)获取多个样本的测序数据；
7.(2)对于每个样本，对于待检测区域，计算该待检测区域的平均测序深度，并计算全基因组的平均测序深度，将所述待检测区域平均深度除以所述全基因组平均深度得到该区域相对于全基因组的拷贝比值；
8.(3)将所述多个样本的拷贝比值从小到大排序，并且从最小值开始每个值依次与后一值比较，如果前后两个值的差小于第一阈值，则将所述样本合并为一组，所有样本进行分组；
9.(4)对于每个分组，根据其拷贝比值与拷贝数参考比值的距离确定最近的拷贝数参考比值，确定为每个分组的拷贝数。
10.在一个实施方案中，在(1)中，所述多个样本多于50个样本。
11.在一个实施方案中，在(2)中，该待检测区域的平均测序深度经过gc矫正。
12.在一个实施方案中，在(3)中，所述第一阈值小于0.15，更优选小于0.12，最优选小于0.1。
13.在一个实施方案中，在(3)中，去除样本数目小于第二阈值的分组，以及样本数目大于第二阈值、但分布不符合正态分布的分组。
14.在一个实施方案中，判断分组是否符合正态分布包括，计算每个分组的拷贝比值的平均值、最大值、最小值和方差。
15.在一个实施方案中，所述第二阈值大于25，优选大于30。
16.在一个实施方案中，在(4)中，所述拷贝数参考比值为一组从0开始的等差数列，等差数为1/n，n为物种倍体数值，优选地，所述样本来自2倍体物种，所述离散数从0开始、以0.5进阶。例如，所述拷贝数参考比值包括0、0.5、1、1.5、2。
17.在一个实施方案中，在(4)中，根据分组的拷贝比值的平均值与拷贝数参考比值的距离确定每个分组所属的拷贝数，优选地，所述确定每个分组所属的拷贝数的公式为：拷贝数＝n*距离最近的拷贝数参考比值，n为物种倍体数值。
18.在另一方面，本发明提供了一种基于群体样本深度信息检测基因拷贝数变异的系统，所述系统包括：
19.测序数据获取模块，用于获取样本的测序数据；
20.拷贝比值计算模块，用于基于所述样本的测序数据计算待检测区域的平均测序深度，计算全基因组的平均测序深度，并将所述待检测区域平均深度除以所述全基因组平均深度得到该区域相对于全基因组的拷贝比值；
21.样本拷贝比值分组模块，用于将多个样本的拷贝比值从小到大排序，并且从最小值开始每个值依次与后一值比较，如果前后两个值的差小于第一阈值，则将所述样本合并为一组，所有样本进行分组；
22.拷贝数确定模块，用于根据分组的拷贝比值与拷贝数参考比值的距离确定最近的拷贝数参考比值，确定为所述分组的拷贝数，优选地，所述确定每个分组所属的拷贝数的公式为：拷贝数＝n*距离最近的拷贝数参考比值，n为物种倍体数值。
23.在一个实施方案中，所述系统还包括分组校验模块，用于针对样本数目大于第二阈值的分组进行正态分布检验。
24.在一个实施方案中，所述测序数据获取模块包括测序仪器。
25.在一个实施方案中，拷贝比值计算模块还用于，对该待检测区域的平均测序深度进行gc矫正。
26.在一个实施方案中，所述第一阈值小于0.15，更优选小于0.12，最优选小于0.1。
27.在一个实施方案中，所述分组校验模块用于计算分组的拷贝比值的平均值、最大值、最小值和方差。
28.在一个实施方案中，所述第二阈值大于25，优选大于30，更优选大于50。
29.在一个实施方案中，所述拷贝数参考比值为一组从0开始的等差数列，等差数为1/n，n为物种倍体数值，优选地，所述样本来自2倍体物种，所述离散数从0开始、以0.5进阶。例如，所述拷贝数参考比值包括0、0.5、1、1.5、2。
30.在一个实施方案中，拷贝数确定模块根据分组的拷贝比值的平均值与拷贝数参考比值的距离确定每个分组所属的拷贝数，优选地，所述确定每个分组所属的拷贝数的公式为：拷贝数＝n*距离最近的拷贝数参考比值，n为物种倍体数值。
31.利用本发明的方法可以更接近真实情况，基于群体样本深度信息准确检出cnv不
同型别，计算拷贝数目值。本发明的方法可以准确检出常见频率的cnv，并正常计算包括dup在内的拷贝数目具体数值。本发明的方法可以检出群体样本不同类型的cnv并正确分型。
附图说明
32.通过以下附图对本发明进行说明。
33.图1示例性示出了790例样本该区域拷贝比例分布。
34.图2示例性示出了样本szch0056的chr8:39226335-39388919区域进行深度及拷贝比例分析。
具体实施方式
35.不希望拘囿于任何理论，发明人假设目标区域在检测样本和群体样本可能存在不同型别的拷贝数目，比如cn0、cn1、cn2等等，对于同一个型别如cn2，该型别的样本在该目标区域计算的测序深度比应该服从正态分布，不同拷贝数目的样本的之间分布存在明显差异。
36.首先，计算目标区域测序深度总和，用所述目标区域测序深度总和除以该目标区域长度得到该区域平均测序深度，并优选所述区域平均测序深度经过gc矫正。目标区域长度至少500bp，可以是数万bp的长度，优选已知具有拷贝数变异的区域，预测有拷贝数变异的区域，或者候选具有拷贝数变异的区域。对于长的目标区域，可以对目标区域进行分段处理，例如每段500bp或1k bp，对分段的目标区域分别运用本发明的方法或系统，在计算结果中将具有相同拷贝数变异的相邻段进行合并，获得拷贝数变异的区段全长。gc矫正的方法可以参见pscc的矫正方法，简述如下：划分窗口统计每个窗口的深度值及gc碱基比例，针对gc比例相同的窗口，取这些窗口的深度值的中位数值作为该gc比例下的深度值。针对所有gc比例及对应的深度值进行lasso回归，得到gc比例与对应的回归深度值，实际测序深度基于该区域gc比例的回归深度进行标准化。同时，计算每个样本的全基因组平均测序深度，全基因组平均测序深度可以是一个样本的所有测序序列，而不要求必须涵盖物种的全基因组，事实上基因组存在大量重复序列，这部分序列并不适合进行测序。然后，以待检测区域平均测序深度除以全基因组平均测序深度得到该目标区域相对于全基因组的拷贝比值。
37.然后，将所有样本的拷贝比值进行从小到大排序，从最小值的样本拷贝比值开始，如果第二个样本拷贝比值与第一个样本拷贝比值的差不超过0.1，则将两个样本拷贝比值合并为一组；然后比较第三样本拷贝比值与第二样本拷贝比值，如果大小不超过0.1，则合并为一组，如此类推，逐次合并；当后一样本拷贝比值与前一样本拷贝比值相比差值大于0.1时，则当前分组截止，重新开始下一个分组，直至对所有样本拷贝比值都进行了分组。
38.再后，针对所有样本数目大于30的分组进行正态分布检验质控；计算每组样本拷贝比值的平均值、最大值、最小值和方差。还计算了每组样本拷贝比值的平均值之间距离和每组平均值与参考点(例如，对于二倍体物种而言，0、0.5、1、1.5和2等)的距离。对于n倍体物种而言，离散数的进阶可以于2倍体物种不同，例如按0/n、1/n、2/n、3/n.....这种方式进阶。n可以为1。对于精子或卵子，可以将其基因组理解为单倍体，这时可以将n看作为1，其参考点可以包括0、1、2、3、4
……
。
39.最后，根据每组平均值最接近的与参考点的值判断该每组内样本该目标区域的拷
贝数目；对于n倍体物种，可以根据其染色体的倍数和平均值的归属确定拷贝数目。优选地，所述确定每个分组所属的拷贝数的公式为：拷贝数＝n*距离最近的拷贝数参考比值，n为物种倍体数值。例如，对于二倍体物种而言，如果一组的平均值与0接近则拷贝数目就是0，如果与0.5接近则拷贝数目就是1，以此类推。
40.在另一方面，本发明可以以系统的方式呈现。例如，本发明涉及一种基于群体样本深度信息检测基因拷贝数变异的系统，所述系统包括测序数据获取模块、拷贝比值计算模块、样本拷贝比值分组模块、拷贝数确定模块和优选的分组校验模块。
41.在本发明中，所述测序数据获取模块用于获取样本的测序数据。所述测序数据获取模块包括测序仪器，所述测序仪器可以对样本进行测序获得测序数据。所述测序数据获取模块也可以从他处获取样本的测序数据，例如存储在计算机本地或远程服务器端的测序数据，或者存储在介质例如光盘、软盘、硬盘上的测序数据。
42.在本发明中，所述拷贝比值计算模块，用于基于所述样本的测序数据计算待检测区域的平均测序深度，计算全基因组的平均测序深度，并将所述待检测区域平均深度除以所述全基因组平均深度得到该区域相对于全基因组的拷贝比值。测序深度可以理解为，测序得到的碱基总量(bp)与基因组大小(genome)的比值，它是评价测序量的指标之一。可以对基因组上的一段序列计算测序深度，即对于一段序列测序得到的碱基量(bp)与序列大小的比值。
43.在本发明中，所述样本拷贝比值分组模块，用于将多个样本的拷贝比值从小到大排序，并且从最小值开始每个值依次与后一值比较，如果前后两个值的差小于第一阈值，则将所述样本合并为一组，所有样本进行分组。将样本拷贝比值分组实际上是将目标区域具有不同拷贝数样本分开，使得具有相同拷贝数的样本在同一个组中。阈值的选择很重要，阈值过小会分成过多的组，将本来在一个组里的样本分成不同组里；阈值过大会将本来在不同组里的样本分到一个组里，无论哪种情况都会使得检测的基因拷贝数变异出现偏差。在本发明的一个实施方案中，所述第一阈值小于0.15，更优选小于0.12，最优选小于0.1；优选地，上述阈值适用于2倍体物种。对于精子或卵子，可以将其基因组理解为单倍体，第一阈值可以为小于0.2。对于n倍体物种而言，例如n大于2，所述第一阈值与2倍体物种不同，可以为小于0.05，甚至小于0.03。
44.在本发明中，所述分组校验模块，用于针对样本数目大于第二阈值的分组进行正态分布检验。在本发明中，发明人发现深度计算的群体样本的拷贝比在同一基因型下的正态分布、不同基因型下明显存在差异，基于此可以检出或判断拷贝数目变异。本发明的方法可以基于拷贝比进行计算，也可以基于拷贝数目计算。本发明的方法可以用来检出cnv，也可以用来鉴定cnv的真假。正态分布需要较多的样本量，个位数的样本量很难观察出来分组的正态分布。在一个实施方案中，所述第二阈值大于25，优选大于30，更优选大于50。正态分布检验的指标用于辅助判断结果的准确性，若分组小于25则不对分组进行检验，该辅助指标输出为空。
45.在本发明中，所述拷贝数确定模块，用于根据分组的拷贝比值与拷贝数参考比值的距离确定最近的拷贝数参考比值，确定为所述分组的拷贝数。计算每组样本拷贝比值的平均值、最大值、最小值和方差。还计算了每组样本拷贝比值的平均值之间距离和每组平均值与参考点(例如，对于二倍体物种而言，0、0.5、1、1.5和2等)的距离。最后，根据每组平均
值最接近的与参考点的值判断该每组内样本该目标区域的拷贝数目；对于n倍体物种，可以根据其染色体的倍数和平均值的归属确定拷贝数目。例如，对于二倍体物种而言，如果一组的平均值与0接近则拷贝数目就是0，如果与0.5接近则拷贝数目就是1，以此类推。
46.所述拷贝比值计算模块、所述样本拷贝比值分组模块、所述分组校验模块和所述拷贝数确定模块可以通过计算机程序实现，例如通过编程语音，根据测序深度的计算公式编写计算机可执行程序，然后借助计算机硬件实现。
47.本领域的技术人员应当理解，本发明的基于群体样本深度信息检测基因拷贝数变异的方法中各个步骤划分以及次序仅仅是示意性地而非限定性的，本领域的技术人员可以在不背离在附加的权利要求和其等价技术方案阐述的本发明精神和范围的情况下，进行删减、增加、替换、修改和变化。
48.本发明可以实现为一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序在由处理器执行时导致本发明的方法的步骤被执行。在一个实施例中，所述计算机程序被分布在网络耦合的多个计算机设备或处理器上，以使得所述计算机程序由一个或多个计算机设备或处理器以分布式方式存储、访问和执行。单个方法步骤/操作，或者两个或更多个方法步骤/操作，可以由单个计算机设备或处理器或由两个或更多个计算机设备或处理器执行。一个或多个方法步骤/操作可以由一个或多个计算机设备或处理器执行，并且一个或多个其他方法步骤/操作可以由一个或多个其他计算机设备或处理器执行。一个或多个计算机设备或处理器可以执行单个方法步骤/操作，或执行两个或更多个方法步骤/操作。
49.实施例1
50.本发明以人类chr20染色体1561287-1594229作为目标区域为例对本发明进行示例性描述。
51.首先，针对790例样本的测序数据，对于每个样本，计算该目标区域的测序深度总和，用所述目标区域测序深度总和除以该目标区域长度得到该区域平均测序深度，并所述区域平均测序深度经过gc矫正。同时，计算每个样本的全基因组平均测序深度，待检测目标区域平均测序深度除以全基因组平均测序深度得到该目标区域相对于全基因组的拷贝比值，790例样本该目标区域拷贝比例分布如图1所示；
52.然后，所有样本的拷贝比值经过从小到大排序，从最小值的样本拷贝比值开始，如果第二个样本拷贝比值与第一个样本拷贝比值的差不超过0.1，则将两个样本拷贝比值合并为一组；然后比较第三样本拷贝比值与第二样本拷贝比值，如果大小不超过0.1，则合并为一组，如此类推，逐次合并；当后一样本拷贝比值与前一样本拷贝比值的相比差值大于0.1时，则当前分组截止，重新开始下一个分组，直至对所有样本拷贝比值都进行了分组；
53.再后，针对所有样本数目大于30的分组进行正态分布检验质控；计算每组样本拷贝比值的平均值、最大值、最小值和方差。还计算了每组样本拷贝比值的平均值之间距离和每组平均值与参考点(0、0.5、1、1.5和2等)的距离；
54.最后，790样本的该目标区域拷贝比值被分成了三组；第一组平均值0.067，标准差0.013，与第二组距离0.561，与参考点中的0比较接近，所以该组所有样本该目标区域为0拷贝；第二组平均值0.627，标准差0.031，与第一组距离0.561，与第三组距离0.560，与参考点中的0.5比较接近，所以该组所有样本该区域为1拷贝；第三组平均值1.187，标准差0.569，
与第二组距离0.560，与参考点中的1比较接近，所以该组所有样本该区域为2拷贝。
55.实施例2
56.对样本szch0056的chr8:39226335-39388919区域进行测序深度及拷贝比例分析，与群体样本分布如图2所示。该样本拷贝比例为0.022，群体中分布在0参考系左右，所以拷贝数目为0。本发明对该区域利用pscc算法进行检测但未检出，经affymetrix公司cytoscan750k的cma芯片验证为真。
57.以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
58.尽管已经参考示例性实施例描述了本发明，但是应理解，本发明并不限于上述实施例的构造和方法。相反，本发明意在覆盖各种修改例和等同配置。另外，尽管在各种示例性结合体和构造中示出了本发明所公开的各种元件和方法步骤，但是包括更多、更少的元件或方法的其它组合也落在本发明的范围之内。