CN107423534B - 基因组拷贝数变异的检测方法和系统 - Google Patents

Abstract

本发明提供了用于对受测者基因组中的拷贝数变异进行分析的方法。本发明的方法创新性地对从测序数据中提取的变量进行数学变换，导出新的二维变量,即测序深度比SDR和交替等位基因频率AAF或者交替等位单倍体频率AHF，进一步提高了信号/噪音比，提高了检测灵敏度和精度，更充分地利用了高通量DNA测序数据所包含的信息。另外，本发明的方法用于处理受测者体液样品中游离DNA，从测序数据中提取基因组某位点的CNV的特征信号时引入临近位点的信息，极大的提高了信号强度，可以有效地检测尿液或血清中游离DNA中的拷贝数变异。本发明还提供了用于实现上述方法的系统。

Description

基因组拷贝数变异的检测方法和系统

技术领域

本发明涉及基因组研究和疾病诊断治疗领域。具体的，本发明提供了一种用于对受测者(例如哺乳动物，特别是人)样品基因组中的拷贝数变异进行分析的方法和相关系统。

背景技术

基因组异常通常与各种遗传疾病、退行性疾病以及癌症关联。例如，癌症中基因拷贝的缺失或增加与基因片段或特定区域的缺失或扩增屡见不鲜。因此肿瘤发生的研究与研发更好的诊断与预后方法都对关联癌症和各种遗传疾病的特定基因区域的鉴定与克隆感兴趣。

癌症的核心特征是基因组的体突变(somatic mutation)，这些突变可以是点突变、染色体结构突变或者是染色体区段的拷贝数突变(copy number variation)。癌症细胞凋亡后，其含有变异的DNA分子，也会被释放到血液中和尿液中。于是，通过检测血液或者尿液中的DNA分子，就可以探知是否有带有变异的DNA存在，从而诊断癌症。

癌症的早期诊断，是治愈或者控制癌症的关键。目前很多癌症(比如肝癌)并没有灵敏的早期诊断手段；或者，诊断手段费用较高(比如PET-CT)；或者，诊断手段有一定的副作用(比如PET-CT导致受试者接受一定剂量的辐射)；或者，诊断手段仅仅在大型医疗机构才能实施(比如PET-CT检测需要受试者亲自前往大型医院)；或者，检测结果不能给予确定性的诊断(比如PET-CT或B超的成像，即使有经验的医生也不能给予确定性的诊断定论)。

检测样品(例如体细胞、血液或者尿液)的核酸的方法在近年有很大的发展，已发展出对核酸的高通量测序方法，能够获得大量病人或同一病人的大量的和完整的基因组信息。这种方法需要分析600亿或更多的序列数据点以提供一个精确的基因组序列。在诊断性基因组测序中，临床诊断的精确度要求进一步地加剧了序列分析的计算复杂性。早期的测序方法中通过从数以千计的孤立的、非常长的DNA片段中产生序列数据，从而保留序列信息的语境完整性并且减少精确数据所需的冗余测试来处理这一复杂性。

拷贝数变异(copy-number variant，CNV)是指在人类基因组中存在的大量核酸片段多态，包括片段的插入、缺失、重复等。这种多态也被称为拷贝数多态(copy-numberpolymorphism，CNP)。CNV发生的频率高于染色体结构变异，而且在整个基因组中覆盖的核苷酸总数超过单核苷酸多态性(single nucleotide polymorphism，SNP)的总数。CNV可能和表型变异紧密关联，同时在物种的演化和发展中发挥着重要作用。

现在已研发了一些检测DNA序列中拷贝数变异的方法。但这些方法中还存在所用的数据相对地易于出错，具有人为偏差等问题。因此，本领域还需要一种更准确、更全面地检测受测者的DNA序列中拷贝数变异的方法。

发明内容

本发明针对以上问题，提供了一种新的更准确的对基因组核酸序列中拷贝数变异(CNV)进行检测和分析的方法，对从测序数据中提取的变量进行数学变换，导出新的二维变量,即LSDR和AAF(或者AHF),进一步提高了信号/噪音比，提高了检测灵敏度和精度，更充分地利用了高通量DNA测序数据所包含的信息。另外，本发明的方法在处理体液中游离DNA样品中从测序数据中提取基因组某位点的CNV的特征信号时，引入临近位点的信息，极大的提高了信号强度，可以有效地检测尿液或血清在游离DNA中的CNV。本发明还提供了用于实现上述方法的系统。

在本发明的其中一个方面，提供了一种用于对受测者(例如哺乳动物，特别是人)样品基因组中的拷贝数变异进行分析的方法，所述方法适合用于检测待测者或患者的各种组织的细胞，包括上皮细胞、血液中的正常白细胞等，所述方法包括以下步骤：

(a)采集参考样本(其中参考样本数量为N，N>50，优选N>100。样本数越大越好)，对其核酸进行测序，在M个单核苷酸多态性(SNP)位点上检测每个样本的基因型(genotype)，对于每一个SNP位点，参考等位基因(reference allele)记为A，交替等位基因(alternative allele)记为B，三种基因型为AA、AB和BB；

对第i个SNP(i为1，2，3，…，M)，获得其参考等位基因和交替等位基因的测序深度d_iA和d_iB，计算得到这个SNP位点的测序深度d_i，其中d_i＝d_iA+d_iB，其交替等位基因B占测序深度的比值为θ_i＝d_iB/(d_iA+d_iB),θ_i∈[0,1]；当θ_i靠近0、0.5和1时对应AA、AB和BB基因型；

以d为横坐标，θ为纵坐标做图，所述N个样本的每个SNP的d_i和θ_i值在图上形成相对于AA、AB和BB基因型的三个点簇，算出三个点簇的中心位置的D_g和Θ_g值，其中D_g和Θ_g值可为点簇上所有点的d_i或θ_i的平均值或中位数值；

(b)通过线性插值法计算待测样本的每一个SNP在个体的正常状态(在哺乳动物，例如人，中，正常状态即基因拷贝数为2，基因型为AA，AB，BB型)下的测序深度d的期望值d_e；

(c)对待测样本的核酸进行测序，获得其每个SNP测序深度；

(d)计算得到待测样本的每个SNP位点的测序深度比(Sequencing Depth Ratio，SDR)，其对数形式为LSDR；以及计算得到交替等位基因频率(Alternative AlleleFrequency，AAF)，其中

SDR＝d/d_e；

LSDR＝log₂(d/d_e)；

以及

(e)采用隐马尔可夫模型分析上述步骤得到的数据，检查目标基因组的拷贝数变异状况。

在本发明的其中一个方面，所述隐马尔可夫模型分析包括以下步骤：

在M个SNP位点的AAF和LSDR值记为x_i和y_i(i＝1,...,M)，构成隐马尔可夫模型的显层(observation layer)；记M个SNP位点的拷贝数变异CNV状态为z_i(i＝1,...,M)，构成隐马尔可夫模型的隐层(hidden layer)，其中每个z_i有6个可能的状态(如以下表1所示)：

表1

从隐层到显层的发射概率为：

p(x_i,y_i|z_i)＝p(x_i|z_i)p(y_i|z_i)。

对应AAF的发射概率为：

其中b(g；G(z)-1,p_B)为二项分布(binomial distribution)的概率密度函数(probability density function)：

φ(x；μ,σ)为正态分布(normal distribution)概率密度函数：

其中G(z)为该SNP位点CNV状态z对应的基因型数(即基因型可能的取值个数)，对应的基因型如上表所示；g＝0,…，G(z)-1，为这G(z)个基因型中B等位基因的个数；当g＝0或g＝G(z)，即CNV基因型为纯合基因型时，对应的正太概率密度函数替换成以0或1为边界的截断(truncated)正态分布概率密度函数；

p_B为受测者所属基因种群的群体B等位基因频率(population frequency of Ballele)；

μ_x,z,g为不同CNV状态基因型对应的AAF均值，按如下公式计算初始值：

σ_x,z,g为不同CNV状态基因型对应的AAF标准差，初始值通过待测样本的数据估计。

对应LSDR的发射概率为：

其中φ(y；μ,σ)为正太分布概率密度函数，参数μ_y,z和σ_y,z为不同CNV状态对应的LSDR均值和标准差；μ_y,z的初始值用如下公式计算：

其中C(z)为不同CNV状态对应的总拷贝数(如表1所示)；

σ_y,z的初始值通过待测样本的数据估计。

隐层相邻SNP位点之间的转移概率为：

其中l为相邻两个SNP位点之间的距离(单位bp)，对于正常状态(即z＝4)到其他状态的转移，L取值范围为10-100Mb，优选为100Mb；对于其他状态之间的转换L取0.1～1Mb，优选为0.1Mb；

p_s,t为基础状态转移概率，当s＝t＝4时，即相邻两个SNP为正常状态，p_4,4取值为非常接近于1的概率，如0.999995；当s＝t、s≠4且t≠4时，即相邻两个SNP为相同的拷贝数变异状态，p_s,t取值接近于1的概率，但小于p_4,4，如0.95；当s≠4且t＝4时，即从拷贝数变异状态转移到正常状态，p_s,t取值为较小的概率，如0.049996；当s和t为其他组合时，即从正常状态转移到拷贝数变异状态或不同拷贝数状态之间的转移，p_s,t取值为极小的概率，如0.000001。

在本发明的方法中，计算σ_x,z,g和σ_y,z的初始值的具体步骤可以为：

(1)首先定义用差分计算标准差的方法：

不失一般性，对于一列2n个(对于奇数个观测值的情况，可以不用最后一个观测值)SNP的观测值(a₁,...,a_2n)，它们按SNP在染色体上的位置排列,数值的分布符合以下性质，即

a_i＝μ_k+ε_i,i∈[i_k-1,i_k),k＝1,...,K,

这里i_k为变点(change point)且变点数K<<2n；ε_i为独立同分布噪声，均值为0，标准差为σ。定义Δ_i＝z_2i-1-z_2i，那么Δ_i(i＝1,...,n)为独立同分布，均值为0，标准差为

Δ_i的标准差可以用较为鲁棒的平均绝对偏差(median absolute deviation，MAD)估计，即

于是a_i的标准差估计为

(2)估计σ_x,z,g的初始值时，先取待测样本AAF值x_i(i＝1,...,M)的子集{x_i:0.2<x_i<0.8}，即选择那些杂合基因型SNP对应的AAF值；然后再对这个子集应用(1)的方法，得到的估计记为σ_x；于是对于不同CNV状态基因型的σ_x,z,g的初始估值为

(3)估计σ_y,z的初始值时，直接对于待测样本的LSDR值y_i(i＝1,...,M)应用(1)的计算方法。

上述公式中的参数(μ_x,z,g、σ_x,z,g、μ_y,z和σ_y,z)在得到初始值后，根据待测样本AAF和LSDR的观测数据，由Baum-Welch算法对待测样本的参数进一步准确估计。

由Viterbi算法估计隐层各SNP位点的CNV状态。由此得到整个目标基因组中各拷贝数变异位置(包括拷贝数变异的起始和结束位点)，拷贝数变异的长度和变异后的拷贝数。

在本发明的另一个方面，还提供了一种用于对受测者(例如哺乳动物，特别是人)样品(特别是体液样品，例如血清和尿液)基因组中的拷贝数变异进行分析的方法。所述方法特别适合用于检测待测者或患者的体液中游离DNA上的拷贝数变异。体液包括血清，尿液等。该方法包括以下步骤：

(a)分别采集目标个体的体液样品以及含有正常基因组的组织样品(如血沉棕黄层或口腔上皮细胞)，从体液样品分离获得体液游离核酸，从含有正常基因组的组织样品分离得到参考核酸，对所述体液核酸和参考核酸进行测序，在M个SNP位点上检测其基因型(genotype)，对于每一个SNP位点，参考等位基因记为A，变异等位基因记为B，三种基因型为AA、AB和BB；

(b)计算得到目标个体的每个杂合型SNP位点的测序深度比SDR，其对数形式为LSDR，其方法如下：

获得每个杂合型SNP位点附近2L宽度区域内(SNP上游宽度为L和下游宽度为L；例如L＝1000bp)每个碱基位置的测序深度，参考核酸和体液核酸在第j个碱基位置上的测序深度记为d_N,i,j和d_T,i,j；

对第i个杂合型SNP(i为1，2，3，…，M)，参考核酸和体液核酸的加权测序深度分别为：

其中w(l_i,j)为权重；

l_i,j为区域中第j个碱基距离SNP_i位点的距离(单位bp)；σ的取值范围为0.3L至L；

第i个杂合型SNP位点的LSDR值计算为：

(c)计算得到目标个体的每个杂合型SNP位点对应的交替等位单倍体频率(Alternative Haplotype Frequency，AHF)，其方法如下：

获得目标个体的体液游离核酸中第i个杂合型SNP周围宽度为2R范围内的杂合基因型SNP位点(上游宽度为R,下游宽度为R；例如R＝50Kb)的基因型；对基因型数据做单倍型定相(haplotype phasing)，判断临近的杂合型SNP位点的两个等位基因A和B各处在哪个单倍体(haplotype)上；在体液游离DNA中，在A和B单倍体上第i个杂合型SNP位点附近其他杂合型SNP位点的测序深度记为d_T,A,i,j和d_T,B,i,j；

计算第i个杂合型SNP在A和B单倍体上的加权测序深度：

其中权重可以设定为多种形式，比如高斯核(Gaussian Kernel)

l_i,j为区域中第j个其他杂合型SNP位点距离SNP_i位点的距离(单位bp)；σ的取值范围在0.2R至R；

体液游离DNA中SNP_i位点的AHF值为：

参考DNA中SNP_i位点的AHF值为：

进一步计算：

LmAHF_i＝log₂(mAHF_T,i/mAHF_N,i)

其中mAHF_N,i＝|AHF_N,i-0.5|+0.5，mAHF_T,i＝|AHF_T,i-0.5|+0.5；

(d)计算出M个杂合SNP位点的LmAHF值(简记为x)和LSDR值(简记为y)后，得到一组2维序列{(x_i,y_i),i＝1,...,M}；采用协同分段(joint segmentation)方法，来确定血液游离DNA中的CNV片段(segment)，即把{(x_i,y_i),i＝1,...,M}分为K个片段

{η_k＝{(x_i,y_i):i∈(t_k-1,t_k]}；k＝1,...,K}，0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M

其中每个片段(η_k)中的杂合型SNP位点({i:i∈(t_k-1,t_k]})对应同一个CNV状态，相邻的片段对应不同的CNV状态；

(e)对于所述K个片段，推断其CNV状态，例如哪些片段属于正常的拷贝数状态(normal)，哪些片段有拷贝数扩增(gain)，缺失(loss)或杂合性缺失(loss ofheterozygosity，LOH)。

本发明的用于对受测者(哺乳动物，如人)的样品基因组中的拷贝数变异进行分析的方法可用于与生物遗传性相关的研究，例如用于研究人种的遗传信息与表观、高度和颜色的关系等。另一些应用包括筛选为了增强或显示期望特征而繁殖的动物等。

用于对受测者(哺乳动物，人)的样品基因组中的拷贝数变异进行分析的方法可以用于鉴定SNP和遗传表型之间的物理连锁关系，由此用于研发遗传图谱，确定对表型(包括疾病表型)重要的基因组区域。另外，本领域已知各种与基因突变有关的疾病。这些疾病包括哮喘、癌症、自身免疫疾病等。癌症包括膀胱癌、脑癌、乳腺癌、结肠癌、食道癌、白血病、肝癌、肺癌、口腔癌、胃癌等。用于对受测者(哺乳动物，人)的样品基因组中的拷贝数变异进行分析的方法也可用于观察、判断和监视与基因突变有关的疾病。

本发明还提供了实现以上用于对受测者(哺乳动物，如人)的样品基因组中的拷贝数变异进行分析的方法的系统。本发明的系统包含计算机处理器以及与所述处理器连接的计算机可读的存储介质，所述存储介质具有明确呈现其上的指令，当由计算机处理器执行时，所述指令引起处理器进行实现以上用于对受测者(哺乳动物，如人)的样品基因组中的拷贝数变异进行分析的方法的各个步骤的操作。

在本发明的其中一个方面，本发明提供的系统包括以下装置中的一种或多种：

(a)一种装置，其获取采集的参考样本的核酸的测序结果，在其SNP位点上检测每个样本的基因型，对于每一个SNP位点，参考等位基因记为A，交替等位基因记为B，三种基因型为AA、AB和BB，获得SNP位点上的参考等位基因和交替等位基因的测序深度d_A和d_B，计算得到这个SNP位点的测序深度d，其中d＝d_A+d_B，记其交替等位基因B占这个SNP位点的测序深度的比值为θ，即θ＝d_B/(d_A+d_B),θ∈[0,1]；

参考样本的个数设为N，检测的SNP位点个数设为M，对第i个SNP(i为1，2，3，…，M)，获得其参考等位基因和交替等位基因的测序深度d_iA和d_iB，计算得到这个SNP位点的测序深度d_i，其中d_i＝d_iA+d_iB，其交替等位基因B占测序深度的比值为θ_i＝d_iB/(d_iA+d_iB),θ_i∈[0,1]；当θ_i靠近0、0.5和1时对应AA、AB和BB基因型；

以d为横坐标，θ为纵坐标做图，所述N个样本的每个SNP的d_i和θ_i值在图上形成相对于AA、AB和BB基因型的三个点簇，算出三个点簇的中心位置的D_g和Θ_g值，其中D_g和Θ_g值可为点簇上所有点的d_i或θ_i的平均值或中位数值：

D_g＝median(d_j,j∈{i:g_i＝g}),g＝AA,AB,BB；

Θ_g＝median(θ_j,j∈{i:g_i＝g}),g＝AA,AB,BB。

(b)一种装置，其通过线性插值法计算待测样本的每一个SNP在个体的正常状态下的测序深度d的期望值d_e；

(c)一种装置，其获取待测样本的核酸的测序结果，获得其每个SNP测序深度；

(d)一种装置，其计算得到待测样本的每个SNP位点的测序深度比，即SDR，其对数形式为LSDR；以及计算得到交替等位基因频率，即AAF，其中

SDR＝d/d_e；

LSDR＝log₂(d/d_e)；

以及

(e)一种装置，其采用隐马尔可夫模型分析上述步骤得到的数据，检查和展示目标基因组的拷贝数变异状况。

在本发明的又一个方面，本发明提供用于对受测者(哺乳动物，例如人)的体液样品(例如血清和尿液样品)基因组中的拷贝数变异进行分析的系统，其包括以下装置中的一种或多种：

(a)一种装置，其获得目标个体的体液游离核酸和含有正常基因组的组织样品的参考核酸的测序结果，在其SNP位点上检测其基因型，对于每一个SNP位点，参考等位基因记为A，变异等位基因记为B，三种基因型为AA、AB和BB；

(b)一种装置，其计算得到目标个体的每个杂合型SNP位点的测序深度比SDR，其对数形式为LSDR，其方法如下：

获得每个杂合型SNP位点附近2L宽度区域内(SNP上游宽度为L和下游宽度为L；L范围为300-2000bp，例如L＝1000bp)每个碱基位置的测序深度，参考核酸和体液游离核酸在第j个碱基位置上的测序深度记为d_N,i,j和d_T,i,j；

检测的SNP位点个数设为M，对第i个杂合型SNP(i为1，2，3，…，M)，参考核酸和体液游离核酸的加权测序深度分别为：

其中w(l_i,j)为权重；优选的，所述权重为高斯核(Gaussian Kernel)

σ的取值范围为0.3L至L；

l_i,j为区域中第j个碱基距离SNP_i位点的距离，

第i个杂合型SNP位点的LSDR值计算为：

(c)一种装置，其计算得到受测者的每个杂合型SNP位点对应的交替等位单倍体频率，即AHF，其方法如下：

获得受测者的体液游离核酸中第i个杂合型SNP周围宽度为2R范围内的杂合基因型SNP位点(上游宽度为R,下游宽度为R；R范围为10-100kb，例如R＝50Kb)的基因型；对基因型数据做单倍型定相，判断临近的杂合型SNP位点的两个等位基因A和B各处在哪个单倍体上；在血清游离DNA中，在A和B单倍体上第i个杂合型SNP位点附近其他杂合型SNP位点的测序深度记为d_T,A,i,j和d_T,B,i,j；

计算第i个杂合型SNP在A和B单倍体上的加权测序深度：

其中w(l_i,j)为权重；优选的，所述权重为高斯核(Gaussian Kernel)

σ的取值范围为0.2R至R；

l_i,j为区域中第j个其他杂合型SNP位点距离SNP_i位点的距离；

血清游离DNA中SNP_i位点的AHF值为：

参考DNA中SNP_i位点的AHF值为：

进一步计算：

LmAHF_i＝log₂(mAHF_T,i/mAHF_N,i)

其中mAHF_N,i＝|AHF_N,i-0.5|+0.5，mAHF_T,i＝|AHF_T,i-0.5|+0.5；

(d)一种装置，其计算出M个杂合SNP位点的LmAHF值(简记为x)和LSDR值(简记为y)后，得到一组2维序列{(x_i,y_i),i＝1,...,M}；采用协同分段方法，来确定体液游离DNA中的CNV片段，即把{(x_i,y_i),i＝1,...,M}分为K个片段

{η_k＝{(x_i,y_i):i∈(t_k-1,t_k]}；k＝1,...,K}，0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M；

其中每个片段(η_k)中的杂合型SNP位点({i:i∈(t_k-1,t_k]})对应同一个CNV状态，相邻的片段对应不同的CNV状态；

(e)一种装置，其对于所述K个片段，推断和显示其CNV状态，例如推断和显示哪些片段属于正常的拷贝数状态，以及哪些CNV片段有拷贝数扩增(gain)，缺失(loss)或杂合性缺失(loss of heterozygosity，LOH)。

附图说明

图1本发明的一种实施方式的对SNP位点的d和θ值进行正则化的示例图。图1A是从测序结果得到的d和θ的散点图；图1B是经过变换，基于d和θ计算出了LSDR和AAF的散点图。

图2本发明的一种实施方式的检测受试者正常细胞的拷贝数变异的医学设备系统的显示图。

图3本发明的一种实施方式的检测肝癌病人的血液中游离DNA的拷贝数变异的医学设备系统的显示图。图3A为杂和型SNP位点的LSDR和LmAHF值的图。图3B为mAHF-LSDR值的图。

具体实施方式

实施例1 获得测试者的组织样品并提取DNA

从测试者抽取10毫升静脉血，收集至含EDTA的收集管，使用从血液中提取核酸的试剂盒(Qiagen QIAamp DSP DNA Blood Mini Kit)提取血清中的游离DNA。该试剂盒的操作方法在网页https://www.qiagen.com/us/shop/sample-technologies/dna/dna-preparation/qiaamp-dsp-dna-blood-mini-kit#orderinginformation进行了详述。

另外，采用口腔拭子在口腔内壁来回挂拭多次，获得患者口腔上皮样本。使用提取核酸的试剂盒(Promega Maxwell RSC Buffy Coat DNA Kit试剂盒)提取口腔上皮的DNA。该试剂盒的操作方法在网页https://www.promega.com/products/dna-purification-quantitation/genomic-dna-purification/maxwell-rsc-system-dna-purification-kits/maxwell-rsc-buffy-coat-dna-kit/进行了详述。

实施例2 对测试者的DNA进行全基因组测序，对其SNP位点进行基因型测定

将5ug的DNA用超声法破碎成大约200bp的片段，使用illumina公司的Paired-endSequencing Sample Preparation Kit试剂盒建立文库，利用illumina HiSeq平台进行全基因组测序。测序得到4千万至4亿个200～300bpDNA片段的序列。

实施例3 对组织细胞，例如口腔上皮细胞DNA的测序和基因型测定结果进行分析

本发明的方法可用于检测待测者或患者的各种组织的细胞，包括上皮细胞、血液中的正常白细胞等。

步骤(a)采集一批正常人群的参考样本(reference samples)，并根据实施例1和实施例2的DNA采集、库准备、测序仪、测序深度相同的方法对参考样本的DNA进行测序和基因型测定。正常细胞DNA样品中，不含有癌症细胞的基因组信息。在正常细胞DNA中检测到的基因组拷贝数多态性(polymorphism)为患者遗传得到的(inherited)。在本发明的方法中，对参考样品所用的测序流程和待检测样品使用的测序流程(包括DNA采集、库准备、测序仪、测序深度)相同或基本相似，这样可以保持在后续所述的对于参考样品和测试样品进行数据变换(transformation)时的准确性。

设参考样品的样本量为N(N需要至少50个以上，优选为100个以上，能够获得越多参考样品，其后分析结果越准确)，在目标基因组内的所以已知SNP位点(设为M个)上测出每个样本的每个SNP的基因型(genotype)，并获得其参考等位基因(reference allele)和交替等位基因(alternative allele)的测序深度(read depth)，即每个SNP位点的测序次数。目标基因组可以是受测者的基因组全长，也可以是选择的基因组片段。

一个样品的测序深度是指测序得到的总碱基数与待测基因组大小的比值。在本方法中，每个SNP位点的测序深度是指对每个SNP位点的测序次数。对基因组上某个位点，包括SNP位点的测序深度受到测序方法和条件以及与该位点自身状况(包括该位点的碱基类型，周围碱基的类型，在基因组上的位置、拷贝数等)的影响。在测序方法和条件类似的情况下，对不同样品的同一个SNP位点的测序深度(即测序次数)主要由该位点的碱基类型，周围碱基的类型，在基因组上的位置和拷贝数决定。

核酸的单核苷酸多态性(single nucleotide polymorphism，SNP)是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。它是人类可遗传的变异中最常见的一种，占所有已知多态性的90％以上。SNP所表现的多态性只涉及到单个碱基的变异，这种变异可由单个碱基的转换或颠换所引起，也可由碱基的插入或缺失所致。SNP既可能是二等位多态性，也可能是3个或4个等位多态性，后两者非常少见。因此，通常所说的SNP都是二等位多态性的。这种变异可能是转换(C←→T，在其互补链上则为G←→A)，也可能是颠换(C←→A，G←→T，C←→G，A←→T)。转换的发生率总是明显高于其它几种变异，具有转换型变异的SNP约占2/3，其它几种变异的发生几率相似。

对于每一个SNP位点，参考等位基因记为A，交替等位基因记为B，可能的三种基因型为AA、AB和BB。对于参考样本i(i＝1,...,N)，其基因型记为g_i(g_i∈{AA,AB,BB})，根据测序的结果直接获得等位基因A和B的原始测序深度，分别记为

和

计算得到这个SNP位点的原始测序深度

其中

为了使不同样本之间的测序深度更准确地进行比较，在本发明的方法中，还可通过每个样本的总测序深度来正则化(normalize)根据测序的结果直接获得的原始测序深度，即通过以下方式：

其中C_i为第i个样本的总测序深度，即所有M个SNP位点的原始测序深度之和；

为N个样本C_i的平均值，即

在这个SNP位点的正则化后的测序深度为d_i＝d_iA+d_iB。

交替等位基因B占测序深度的比值为θ_i＝d_iB/(d_iA+d_iB)(θ_i∈[0,1])。当θ_i靠近0、0.5和1时对应AA、AB和BB基因型。

以d为横坐标，θ为纵坐标做图，所述N个样本的每个SNP的d_i和θ_i值在图上形成相对于AA、AB和BB基因型的三个点簇(cluster)(参见图1A)。分别计算出三个点簇的中心位置的D_g和Θ_g值，其中D_g和Θ_g值可为点簇上所有点的d_i或θ_i的平均值或中位数值。

在本实施方式的方法中，用中位数(median)来估计点簇的中心位置，即

D_g＝median(d_j,j∈{i:g_i＝g}),g＝AA,AB,BB

Θ_g＝median(θ_j,j∈{i:g_i＝g}),g＝AA,AB,BB

通过用中位数(median)来估计点簇的中心位置，更有利于避免异常值对估计的影响。

从参考样本得到的M个SNP的D_g和Θ_g值可以保存用作对于参考样本和新样本进行数据变换(transformation)的基准值。

步骤(b)通过线性插值法计算待测样本的每一个SNP在正常状态(对于人类基因组，正常状态即基因拷贝数为2，基因型为AA，AB，BB型)下的测序深度d的期望值d_e(expectation)。

例如，用以下线性插值的方法计算：

步骤(c)对待测样本的核酸进行测序，获得其每个SNP测序深度d，记为d_i，(i为1，2，3，…，M)；

步骤(d)计算得到待测样本的每个SNP位点的测序深度比(Sequencing DepthRatio，SDR)，其对数形式为LSDR；以及计算得到交替等位基因频率(Alternative AlleleFrequency，AAF)，其中

SDR＝d/d_e；

LSDR＝log₂(d/d_e)；

以及

图1为本发明的方法对SNP位点的d和θ值进行正则化的示例图。图1A是从测序结果得到的d和θ的散点图；经过以上所述变换，基于d和θ计算出了LSDR和AAF，其散点图为图1B。如图1所示，本发明的方法通过对于所有M个SNP位点经过以上数据变换，每个SNP位点的LSDR和AAF值都正则化到统一的尺度，以便于聚集(aggregate)邻近的SNP位点的总信号强度和等位基因比值，提高本发明的CNV检测方法的信噪比。

步骤(e)采用隐马尔可夫模型(Hidden Markov Model，HMM)分析上述步骤得到的测序深度比SDR和LSDR；以及交替等位基因频率AAF，检查目标基因组的拷贝数变异状况。

其中，在M个SNP位点的AAF和LSDR值记为x_i和y_i(i＝1,...,M)，构成了隐马尔可夫模型的显层(observation layer)；记M个SNP位点的拷贝数变异CNV状态为z_i(i＝1,...,M)，构成隐马尔可夫模型的隐层(hidden layer)，其中每个z_i有6个可能的状态(如下表)：

表1

从隐层到显层的发射概率为：

p(x_i,y_i|z_i)＝p(x_i|z_i)p(y_i|z_i)。

对应AAF的发射概率为：

其中b(g；G(z)-1,p_B)为二项分布(binomial distribution)的概率密度函数(probability density function)：

φ(x；μ,σ)为正态分布(normal distribution)概率密度函数：

其中G(z)为该SNP位点CNV状态z对应的基因型数(即基因型可能的取值个数)，对应的基因型如上表所示；g＝0，…，G(z)-1，为这G(z)个基因型中B等位基因的个数；当g＝0或g＝G(z)，即CNV基因型为纯合基因型时，对应的正太概率密度函数替换成以0或1为边界的截断(truncated)正太概率密度函数；

p_B为受测者所属基因种群的群体B等位基因频率(population frequency)。p_B可从已知的参考样本分析和统计后获得。p_B也可从已公开的研究数据获得。选择的基因种群与受测者越接近越合适。例如，受测者为人类时，可以选择其所属种族、家族或是属于感兴趣的基因型群体等。在本实施例中，从千人基因组计划提供的待测样本所属的种群数据得到，例如在其网页http://www.1000genomes.org/data#download得到。例如位于第17号染色体上的，编号为rs2305480的SNP,参考等位基因为G；替换等位基因为A，替换等位基因的频率在东亚人群中为25.8％，在非洲人群中为14.1％，欧洲白人中为45.3％。

μ_x,z,g为不同CNV状态基因型对应的AAF均值。例如，在本实施例的方法中，按如下公式计算初始值：

σ_x,z,g为不同CNV状态基因型对应的AAF标准差，初始值通过待测样本的数据估计，例如通过下面描述的方法估算。

对应LSDR的发射概率为：

同样φ(y；μ,σ)为正太分布概率密度函数，参数μ_y,z和σ_y,z为不同CNV状态对应的LSDR均值和标准差。μ_y,z的初始值可用如下公式计算：

其中C(z)为不同CNV状态对应的总拷贝数(如上表)；

σ_y,z的初始值通过待测样本的数据估计，例如通过下面描述的方法估算。

在本实施例中，计算σ_x,z,g和σ_y,z的初始值的具体步骤为：

(1)首先定义用差分计算标准差的方法：

不失一般性，对于一列2n个(对于奇数个观测值的情况，可以不用最后一个观测值)SNP的观测值(a₁,...,a_2n)，按SNP在染色体上的位置排列,数值的分布符合以下性质，即

a_i＝μ_k+ε_i,i∈[i_k-1,i_k),k＝1,...,K,

其中i_k为变点(change point)且变点数K<<2n；ε_i为独立同分布噪声，均值为0，标准差为σ。定义Δ_i＝z_2i-1-z_2i。那么Δ_i(i＝1,...,n)为独立同分布，均值为0，标准差为

Δ_i的标准差可以用较为鲁棒(robust)的平均绝对偏差(median absolutedeviation，MAD)估计，即

于是a_i的标准差估计为

(2)估计σ_x,z,g的初始值时，先取待测样本AAF值x_i(i＝1,...,M)的子集{x_i:0.2<x_i<0.8}，即选择那些杂合基因型的SNP对应的AAF值；然后再对这个子集应用(1)的方法，得到的估计记为σ_x；于是对于不同CNV状态基因型的σ_x,z,g的初始估值为

(3)估计σ_y,z的初始值时，直接对于待测样本的LSDR值y_i(i＝1,...,M)应用(1)的计算方法。

隐层相邻SNP位点之间的转移概率为：

其中l为相邻两个SNP位点之间的距离(单位bp)。对于正常状态(即z＝4)到其他状态的转移，L取值范围为10-100Mb，优选为100Mb；对于其他状态之间的转换L取0.1～1Mb，优选为0.1Mb。

p_s,t为基础状态转移概率。当s＝t＝4时，即相邻两个SNP为正常状态，p_4,4取值为非常接近于1的概率，如0.999995；当s＝t、s≠4且t≠4时，即相邻两个SNP为相同的拷贝数变异状态，p_s,t取值接近于1的概率，但小于p_4,4，如0.95；当s≠4且t＝4时，即从拷贝数变异状态转移到正常状态，p_s,t取值为较小的概率，如0.049996；当s和t为其他组合时，即从正常状态转移到拷贝数变异状态或不同拷贝数状态之间的转移，p_s,t取值为极小的概率，如0.000001。

上述公式中的参数(μ_x,z,g、σ_x,z,g、μ_y,z和σ_y,z)在得到初始值后，根据待测样本AAF和LSDR的观测数据，由Baum-Welch算法对待测样本的参数进一步准确估计。

然后，由Viterbi算法估计隐层各SNP位点的拷贝数变异状态，由此得到整个目标基因组中各拷贝数变异位置(包括拷贝数变异的起始和结束位点)，拷贝数变异的长度和变异后的拷贝数。

本发明的上述检测受试者正常细胞的拷贝数变异的方法可以通过带计算机系统的医学设备实现，其结果可以通过软件转化为直观的系统显示结果，包括目标基因组中各拷贝数变异位置(包括拷贝数变异的起始和结束位点)，拷贝数变异的长度和变异后的拷贝数等信息。图2为示例性的应用本发明的上述方法检测受试者正常细胞的拷贝数变异后的医学设备系统的显示图。图2显示受试者第15号染色体的一个区段。在15号染色体32Mb碱基位置附近，多个SNP位点的LSDR值偏离正常状态的预测值(LSDR的预测值为0)；同时多个SNP的AAF的值也偏离了预测值(杂合状态的SNP的AAF的预测值为0.5)。使用本发明的方法，可以确定该受试者的第15号染色体，32Mb碱基位置附近有一个大约1.2Mb碱基长的CNV。

实施例4 对游离DNA(例如血清和尿液等体液中的游离DNA)的测序和基因型测定结果进行分析

在癌症患者的血液或尿液中存在肿瘤细胞凋亡后释放的DNA片段，也就带有肿瘤细胞基因组中的体突变(somatic mutation)和体拷贝数变异(肿瘤细胞基因组中的CNV特别称为somatic copy number alteration,SCNA)的信息。但游离DNA样品中的主要成分是正常细胞凋亡后释放的DNA，SCNA的信号被严重稀释。因此需要进一步增强LSDR的AAF的信噪比(signal-to-noise ratio)。

本发明还提供了根据患者的血液游离DNA来检测肿瘤细胞的CNV状态的方法，其步骤如下：

步骤(a)从同一测试者提取的参考DNA(含有正常基因组)和血清游离DNA。从血清样品分离获得血清游离核酸，从含有正常基因组的组织样品(如血沉棕黄层或口腔上皮细胞)分离得到前述参考核酸，对所述血清核酸和参考核酸进行测序，在M个SNP位点上检测其基因型。对于每一个SNP位点，参考等位基因记为A，变异等位基因记为B，三种基因型为AA、AB和BB。

步骤(b)计算得到目标个体的每个杂合基因型(下文简称为杂合型)SNP位点的测序深度比(Sequencing Depth Ratio，SDR)，其对数形式为LSDR，其计算方法如下：

获得每个杂合型SNP位点附近2L宽度区域内(SNP上游宽度为L和下游宽度为L；L范围为300-2000bp，例如L＝1000bp)每个碱基位置的测序深度，参考核酸和血清核酸在第j个碱基位置上的测序深度记为d_N,i,j和d_T,i,j；

对第i个杂合型SNP(i为1，2，3，…，M)，参考核酸和血清核酸的加权测序深度分别为：

其中权重可以设定为多种形式，例如可使用高斯核，或是可以设定为1/2L，即各位点的影响都一样。

在本发明的实施例中，使用高斯核(Gaussian Kernel)

l_i,j为区域中第j个碱基距离SNP_i位点的距离(单位bp)；σ的取值范围为0.3L至L。

第i个杂合型SNP位点的LSDR值计算为：

与实施例3中的方法比较，LSDR的定义在逻辑上是相同的，都是用于描述实际测量到的测序深度与预计测序深度的比值，只是为了进一步增强LSDR的AAF的信噪比而做了以下修正(modification)。实施例3中的方法的每个SNP位点的LSDR只用到了该位点的测序深度，而实施例4的方法用到了当前位点周围其他碱基位置的测序深度。

(c)计算得到目标个体的每个杂合型SNP位点的交替等位单倍体频率(Alternative Haplotype Frequency，AHF)。

获得目标个体的血清游离核酸中第i个杂合型SNP周围宽度为2R范围内的杂合型SNP位点(上游宽度为R,下游宽度为R；R范围为10-100kb，例如R＝50Kb)的基因型。

通过SHAPEIT软件(https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html)，对基因型数据做单倍型定相(haplotype phasing)，判断临近的杂合型SNP位点的两个等位基因各处在哪个单倍体(haplotype)上。哺乳动物中，每个染色体，都有两条。一条来自父亲，一条来自母亲。单倍型，就是单倍体上多态性位点的基因型。对基因型数据做单倍型定相(haplotype phasing)，即为从某个个体的双倍体数据，估算其单倍型。

不失一般性，把两条单倍体分别记作A和B单倍体。在血清游离DNA中，在A和B单倍体上第i个杂合型SNP位点附近其他杂合型SNP位点的测序深度记为d_T,A,i,j和d_T,B,i,j。

计算第i个杂合型SNP在A和B单倍体上的加权测序深度：

其中权重可以设定为多种形式，比如高斯核(Gaussian Kernel)

l_i,j为区域中第j个其他杂合型SNP位点距离SNP_i位点的距离(单位bp)；σ的值在0.2R至R。

血清游离DNA中SNP_i位点的AHF值为：

同样，参考DNA中SNP_i位点的AHF值为：

AHF可视为实施例3中的定义AAF的扩展。AHF的定义和实施例3的方法中的AAF的定义在逻辑上是相同的，都是用于描述二倍体物种细胞，两套染色体(一套来自父本，另一套来自母本)在某个位置的测序深度的比值。实施例3中的AAF只用到单个杂合型SNP位点的等位基因测序深度信息，而实施例4中的AHF用到当前杂合型SNP位点附近其他杂合型SNP位点的等位基因测序深度信息。

进一步计算：

LmAAF＝log₂(mAAF_T,i/mAAF_N,i)

其中mAAF_N,i＝|AAF_N,i-0.5|+0.5，mAAF_T,i＝|AAF_T,i-0.5|+0.5；

通过如上进一步数据变换得到LmAHF，更符合正态分布，更适合应用在以下步骤(d)中的算法。

步骤(d)计算出M个杂合型SNP位点的LmAHF值(简记为x)和LSDR值(简记为y)后，得到一组2维序列{(x_i,y_i),i＝1,...,M}。采用协同分段(joint segmentation)方法，来确定血液游离DNA中的CNV片段(segment)，即把(x_i,y_i),i＝1,...,M分为K个片段：

{η_k＝{(x_i,y_i):i∈(t_k-1,t_k]}；k＝1,...,K}，0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M

其中每个片段(η_k)中的杂合型SNP位点({i:i∈(t_k-1,t_k]})对应同一个CNV状态，相邻的片段对应不同的CNV状态。

分段(segmentation)是统计中变点检测(change-point detection)问题的常用概念。协同分段就是对多维序列进行分段。可采用统计学中各种合适的方法进行上述协同分段。在本实施例的方法中，采用的协同分段的具体算法如下：

(d.1)对于2维序列{(x_i,y_i),i＝0,...,M}(注意：为了算法符号含义清晰简洁，在原有序列前加上辅助的数据点(x₀＝0,y₀＝0))，首先定义统计量：

其中0≤s<t≤M；

和

可由统计学上惯用的方法，例如差分计算法获得。在本实施例中，可以下述方法获得：

不失一般性，对于一列2n个(对于奇数个观测值的情况，可以不用最后一个观测值)SNP的观测值(a₁,...,a_2n)，它们按SNP在染色体上的位置排列，数值的分布符合以下性质，即

a_i＝μ_k+ε_i,i∈[i_k-1,i_k),k＝1,...,K,

其中i_k为变点(change point)且变点数K<<2n；ε_i为独立同分布噪声，均值为0，标准差为σ。定义Δ_i＝z_2i-1-z_2i。那么Δ_i(i＝1,...,n)为独立同分布，均值为0，标准差为

Δ_i的标准差可以用较为鲁棒的平均绝对偏差(median absolute deviation，MAD)估计，即

于是a_i的标准差估计为

然后定义扫瞄统计量(scan statistic)：

(d.2)扫瞄统计量的显著性由以下公式计算：

其中

为自由度为2的卡方分布(chi-square distribution)概率密度函数；v(x)≈[(2/x)(Φ(x/2)-1/2)]/[(x/2)Φ(x/2)+φ(x/2)]；φ(·)和Φ(·)为标准正态分布的概率密度函数和概率分布函数。

(d.3)寻找变点0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M的算法：

步骤1：给定显著性水平α，T₁＝0，T₂＝M；

步骤2：计算

其中(s*,t*]是使Z(s,t)达到最大值的片段；

步骤3：如果p_Z(Z_max)<α，那么片段(T₁,T₂]分为三个(T₁,s*]，(s*,t*]和(t*,T₂](注意若T₁＝s*或t*＝T₂，则实际分为两个片段)；对原有序列进行更新：

其中

步骤4：对于(T₁,T₂]∈{(T₁,s*],(s*,t*],(t*,T₂]}和更新的2维序列{(x_i',y_i'),i＝0,...,M}分别应用步骤2和步骤3，直到所有的片段对应的p_Z>α。

(e)对于这K个片段，可采用统计学上已知的合适方法来推断哪些片段属于正常的拷贝数状态，哪些片段有拷贝数扩增(gain)，缺失(loss)或杂合性缺失(loss ofheterozygosity，LOH)。

在本实施例中，采用以下推断方法：

对于(d)中得到的K个片段，对每个片段(t_k-1,t_k]分别计算：

其中n_k＝t_k-t_k-1；

然后采用自助法(bootstrap)，在指标集{0,1,...,M}中随机选取H个长度为n_k的片段

用相同公式计算上述三个统计量

然后计算

其中I(·)为指示(indicator)函数，即条件为真，则函数值为1，否则为0。

每个片段的CNV状态根据下表2给出。

表2推断基因组片段拷贝数状态的标准。

表中p_xy、p_x、p_y、

和

的片段索引k在这里忽略，因为是对每一个片段都应用上述计算过程；非平衡是指两个等位基因的拷贝数不等；平衡是指两个等位基因的拷贝数相等；通常θ＝0.05，u＝1.5。

图3为示例性的应用本发明的上述方法检测一个肝癌病人的血液中游离DNA的拷贝数变异后的医学设备系统的显示图。图3显示，本发明的算法计算出这例肝癌基因组为四倍体，并且定位和定量基因组上的大量CNV。本发明的算法计算出了每一个基因组区段的总拷贝数，以及父系和母系两条染色体的比例。图3A中，上下两个子图中的点分别表示杂和型SNP位点的LSDR和LmAHF值，数据点用交替的浅绿色和浅灰色来区分它们所在的染色体。每条线段表示一个CNV片段。线段用不同颜色对应不同的CNV类型(基准为四倍体)，例如灰色为正常(normal)、蓝色为拷贝数缺失(loss)、红色为拷贝数增加(gain)，绿色为杂和性缺失(LOH)，浅蓝色为不确定(undecided)。图3B中，在mAHF-LSDR坐标平面上(注意：在图3A中所示的协同分段完成后，后续步骤中的计算只需要用到血清游离DNA的mAHF值)，每个圆圈表示一个CNV片段，分别对应于图3A中的线段，其尺寸与实际CNV片段的尺寸成比例。圆圈的颜色所表示的意义与图3A相同。圆圈簇(cluster)附近的数字标注表示推断出的相应的绝对拷贝数(absolute copy number)/多数等位基因的拷贝数(major allele copy number)。灰色虚线表示修正过的在mAHF维度和LSDR维度的基准线。坐标平面旁边是对应于mAHF维度和LSDR维度的加权直方图，数值为每个CNV片段的mAHF中位数和LSDR中位数，权重为每个CNV片段所含的杂和SNP位点的数量。

除非另外指出，本发明的实践将使用生物技术、计算机科学等的常规技术，显然除在上述说明和实施例中所特别描述之外，还可以别的方式实现本发明。其它在本发明范围内的方面与改进将对本发明所属领域的技术人员显而易见。根据本发明的教导，许多改变和变化是可行的，因此其在本发明的范围之内。本文所提到的所有专利、专利申请与科技论文均据此通过引用结合到本文。

Claims (18)

1.一种用于对受测者样品基因组中的拷贝数变异进行分析的方法，所述方法包括以下步骤：

(a)采集参考样本，对其核酸进行测序，在其SNP位点上检测每个样本的基因型，对于每一个SNP位点，参考等位基因记为A，交替等位基因记为B，三种基因型为AA、AB和BB，获得SNP位点上的参考等位基因和交替等位基因的测序深度d_A和d_B，计算得到这个SNP位点的测序深度d，其中d＝d_A+d_B，记其交替等位基因B占这个SNP位点的测序深度的比值为θ，即θ＝d_B/(d_A+d_B)，θ∈[0,1]；

参考样本的个数设为N，检测的SNP位点个数设为M，对第i个SNP，其中i为1，2，3，…，M，获得其参考等位基因和交替等位基因的测序深度d_iA和d_iB，计算得到这个SNP位点的测序深度d_i，其中d_i＝d_iA+d_iB，其交替等位基因B占测序深度的比值为θ_i＝d_iB/(d_iA+d_iB)，θ_i∈[0,1]；当θ_i靠近0、0.5和1时对应AA、AB和BB基因型；

以d为横坐标，θ为纵坐标做图，所述N个样本的每个SNP的d_i和θ_i值在图上形成相对于AA、AB和BB基因型的三个点簇，算出三个点簇的中心位置的D_g和Θ_g值，其中D_g和Θ_g值可为点簇上所有点的d_i或θ_i的平均值或中位数值：

D_g＝median(d_j,j∈{i:g_i＝g}),g＝AA,AB,BB；

Θ_g＝median(θ_j,j∈{i:g_i＝g}),g＝AA,AB,BB，

(b)通过线性插值法计算待测样本的每一个SNP在个体的正常状态下的测序深度d的期望值d_e；

(c)对待测样本的核酸进行测序，获得其每个SNP测序深度；

(d)计算得到待测样本的每个SNP位点的测序深度比，即SDR，其对数形式为LSDR；以及计算得到交替等位基因频率，即AAF，其中

SDR＝d/d_e；

LSDR＝log₂(d/d_e)；

以及

(e)采用隐马尔可夫模型分析上述步骤得到的每个SNP位点基因型、AAF和LSDR值，检查目标基因组的拷贝数变异状况。

2.如权利要求1所述的方法，其中步骤(e)中所述采用隐马尔可夫模型分析包括以下步骤：

在M个SNP位点的AAF和LSDR值记为x_i和y_i(i＝1,...,M)，构成隐马尔可夫模型的显层；记M个SNP位点的拷贝数变异状态为z_i(i＝1,...,M)，构成隐马尔可夫模型的隐层，其中每个z_i有6个可能的状态：

从隐层到显层的发射概率为：

p(x_i,y_i|z_i)＝p(x_i|z_i)p(y_i|z_i)；

对应AAF的发射概率为：

其中b(g；G(z)-1,p_B)为二项分布的概率密度函数：

φ(x；μ,σ)为正态分布概率密度函数：

其中G(z)为该SNP位点CNV状态z对应的基因型数(即基因型可能的取值个数)，对应的基因型如上表所示；g＝0，…，G(z)-1，为这G(z)个基因型中B等位基因的个数；当g＝0或g＝G(z)，即CNV基因型为纯合基因型时，对应的正太概率密度函数替换成以0或1为边界的截断正态分布概率密度函数；

p_B为群体B等位基因频率；

μ_x,z,g为不同CNV状态基因型对应的AAF均值，按如下公式计算初始值：

σ_x,z,g为不同CNV状态基因型对应的AAF标准差，初始值通过待测样本的数据估计；

对应LSDR的发射概率为：

其中φ(y；μ,σ)为正态分布概率密度函数，参数μ_y,z和σ_y,z为不同CNV状态对应的LSDR均值和标准差；μ_y,z的初始值用如下公式计算：

其中C(z)为不同CNV状态对应的总拷贝数；

σ_y,z的初始值通过待测样本的数据估计；

隐层相邻SNP位点之间的转移概率为：

其中l为相邻两个SNP位点之间的距离，对于正常状态，即z＝4，到其他状态的转移，L取值范围为10-100Mb；对于其他状态之间的转换L取0.1～1Mb；

p_s,t为基础状态转移概率，当s＝t＝4时，即相邻两个SNP为正常状态，p_4,4取值为非常接近于1的概率；当s＝t、s≠4且t≠4时，即相邻两个SNP为相同的拷贝数变异状态，p_s,t取值接近于1的概率，但小于p_4,4；当s≠4且t＝4时，即从拷贝数变异状态转移到正常状态，p_s,t取值为较小的概率；当s和t为其他组合时，即从正常状态转移到拷贝数变异状态或不同拷贝数状态之间的转移，p_s,t取值为极小的概率；

由Viterbi算法估计隐层各SNP位点的拷贝数变异状态。

3.如权利要求2所述的方法，其中对于正常状态，即z＝4，到其他状态的转移，L取值为100Mb，对于其他状态之间的转换L取值为0.1Mb。

4.如权利要求2所述的方法，其中计算所述σ_x,z,g和σ_y,z的初始值的具体步骤为：

(1)定义用差分计算标准差的方法：

对于一列2n个SNP的观测值(a₁,...,a_2n)，它们按SNP在染色体上的位置排列,数值的分布符合以下性质，即

a_i＝μ_k+ε_i,i∈[i_k-1,i_k),k＝1,...,K,

i_k为变点且变点数K<<2n；ε_i为独立同分布噪声，均值为0，标准差为σ；定义Δ_i＝z_2i-1-z_2i，那么Δ_i，其中i＝1,...,n，为独立同分布，均值为0，标准差为

Δ_i的标准差

于是a_i的标准差估计为

(2)估计σ_x,z,g的初始值时，先取待测样本AAF值x_i，其中i＝1,...,M，的子集{x_i:0.2<x_i<0.8}，即选择那些杂合基因型的SNP对应的AAF值；然后再对这个子集应用(1)的方法，得到的估计为σ_x；于是对于不同CNV状态基因型的σ_x,z,g的初始估值为

(3)估计σ_y,z的初始值时，直接对于待测样本的LSDR值y_i，其中i＝1,...,M，应用(1)的计算方法。

5.如权利要求4所述的方法，其中参数μ_x,z,g、σ_x,z,g、μ_y,z和σ_y,z在得到初始值后，根据待测样本AAF和LSDR的观测数据，由Baum-Welch算法对待测样本的参数进一步准确估计。

6.如权利要求1所述的方法，其中计算第i个SNP位点的测序深度d_i时，对获得的原始测序深度进行正则化处理，即

对于样本i，其中i＝1,…,N，其基因型记为g_i，其中g_i∈{AA,AB,BB}，等位基因A和B的原始测序深度记为

和

对原始测序深度进行如下正则化，

其中C_i为第i个样本的总测序深度，即所有M个SNP位点的原始测序深度之和；

为N个样本C_i的平均值，即

该SNP位点的正则化后的测序深度为d_i＝d_iA+d_iB。

7.如权利要求1所述的方法，其中步骤(b)中通过以下线性插值法公式计算得到d_e：

8.如权利要求1所述的方法，其中受测者是人。

9.一种用于对受测者的样品基因组中的拷贝数变异进行分析的方法，所述方法包括以下步骤：

(a)分别采集目标个体的体液样品以及含有正常基因组的组织样品，从体液样品分离获得待测的体液游离核酸，从含有正常基因组的组织样品分离得到参考核酸，对所述体液游离核酸和参考核酸进行测序，在其SNP位点上检测其基因型，对于每一个SNP位点，参考等位基因记为A，变异等位基因记为B，三种基因型为AA、AB和BB；

(b)计算得到目标个体的每个杂合型SNP位点的测序深度比SDR，其对数形式为LSDR，其方法如下：

获得每个杂合型SNP位点附近2L宽度区域内，即SNP上游宽度为L和下游宽度为L，L范围为300-2000bp，每个碱基位置的测序深度，参考核酸和体液游离核酸在第j个碱基位置上的测序深度记为d_N,i,j和d_T,i,j；

检测的SNP位点个数设为M，对第i个杂合型SNP，其中i为1，2，3，…，M，参考核酸和体液游离核酸的加权测序深度分别为：

其中w(l_i,j)为权重，所述权重为高斯核

σ的取值范围为0.3L至L；

l_i,j为区域中第j个碱基距离SNP_i位点的距离，

第i个杂合型SNP位点的LSDR值计算为：

(c)计算得到受测者的每个杂合型SNP位点对应的交替等位单倍体频率，即AHF，其方法如下：

获得受测者的体液游离核酸中第i个杂合型SNP周围宽度为2R范围内的杂合基因型SNP位点，即上游宽度为R,下游宽度为R；R范围为10-100kb，的基因型；对基因型数据做单倍型定相，判断临近的杂合型SNP位点的两个等位基因A和B各处在哪个单倍体上；在血清游离DNA中，在A和B单倍体上第i个杂合型SNP位点附近其他杂合型SNP位点的测序深度记为d_T,A,i,j和d_T,B,i,j；

计算第i个杂合型SNP在A和B单倍体上的加权测序深度：

其中w(l_i,j)为权重；

血清游离DNA中SNP_i位点的AHF值为：

参考DNA中SNP_i位点的AHF值为：

进一步计算：

LmAHF_i＝log₂(mAHF_T,i/mAHF_N,i)

其中mAHF_N,i＝|AHF_N,i-0.5|+0.5，mAHF_T,i＝|AHF_T,i-0.5|+0.5；

(d)计算出M个杂合SNP位点的LmAHF值记为x和LSDR值记为y后，得到一组2维序列{(x_i,y_i),i＝1,...,M}；采用协同分段方法，来确定体液游离DNA中的CNV片段，即把{(x_i,y_i),i＝1,...,M}分为K个片段

{η_k＝{(x_i,y_i):i∈(t_k-1,t_k]}；k＝1,...,K}，0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M；

其中每个片段η_k中的杂合型SNP位点{i:i∈(t_k-1,t_k]}对应同一个CNV状态，相邻的片段对应不同的CNV状态；

(e)对于所述K个片段，推断其CNV状态，所述推断方法如下：

对于(d)中得到的K个片段，对每个片段(t_k-1,t_k]分别计算：

其中n_k＝t_k-t_k-1；

然后采用自助法，在指标集{0,1,...,M}中随机选取H个长度为n_k的片段

用相同公式计算上述三个统计量

然后计算

其中I(·)为指示函数，即条件为真，则函数值为1，否则为0，

其中CNV状态为哪些片段属于正常的拷贝数状态，以及哪些CNV片段有拷贝数扩增，缺失或杂合性缺失。

10.如权利要求9所述的方法，其中片段的CNV状态根据下表给出

表中p_xy、p_x、p_y、

和

的片段索引k在这里忽略，因为是对每一个片段都应用上述计算过程；非平衡是指两个等位基因的拷贝数不等；平衡是指两个等位基因的拷贝数相等；θ＝0.05，u＝1.5。

11.如权利要求9所述的方法，其中步骤(d)中所述协同分段的具体算法如下：

(d.1)对于2维序列{(x_i,y_i),i＝0,...,M}，在原有序列前加上辅助的数据点即x₀＝0,y₀＝0，定义统计量：

其中0≤s<t≤M；

和

由差分计算法得到；

然后定义扫瞄统计量：

(d.2)扫瞄统计量的显著性由以下公式计算：

其中

为自由度为2的卡方分布概率密度函数；v(x)≈[(2/x)(Φ(x/2)-1/2)]/[(x/2)Φ(x/2)+φ(x/2)]；φ(·)和Φ(·)为标准正态分布的概率密度函数和概率分布函数；

(d.3)寻找变点0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M的算法：

步骤1：给定显著性水平α，T₁＝0，T₂＝M；

步骤2：计算

其中(s*,t*]是使Z(s,t)达到最大值的片段；

步骤3：如果p_Z(Z_max)<α，那么片段(T₁,T₂]分为三个(T₁,s*]，(s*,t*]和(t*,T₂](注意若T₁＝s*或t*＝T₂，则实际分为两个片段)；对原有序列进行更新：

其中

步骤4：对于(T₁,T₂]∈{(T₁,s*],(s*,t*],(t*,T₂]}和更新的2维序列{(x_i',y_i'),i＝0,...,M}分别应用步骤2和步骤3，直到所有的片段对应的p_Z＞α。

12.如权利要求11所述的方法，其中步骤(d.1)用差分法计算标准差

和

方法如下：

对于一列2n个，其中对于奇数个观测值的情况，不用最后一个观测值，SNP的观测值(a₁,...,a_2n)，它们按SNP在染色体上的位置排列，数值的分布符合以下性质，即

a_i＝μ_k+ε_i,i∈[i_k-1,i_k),k＝1,…,K,

其中i_k为变点且变点数K<<2n；ε_i为独立同分布噪声，均值为0，标准差为σ，定义Δ_i＝z_2i-1-z_2i，那么Δ_i(i＝1,...,n)为独立同分布，均值为0，标准差为

Δ_i的标准差可以用较为鲁棒的平均绝对偏差估计，即

于是a_i的标准差估计为

13.如权利要求9所述的方法，其中所述w(l_i,j)为高斯核

σ的取值范围为0.3L至L；

l_i,j为区域中第j个碱基距离SNP_i位点的距离。

14.如权利要求9所述的方法，其中受测者是人。

15.如权利要求9所述的方法，其中样品是体液样品。

16.如权利要求9所述的方法，其中样品是体液样品。

17.一种用于实现如权利要求1-8中任一项所述的方法的系统，用于对受测者的样品基因组中的拷贝数变异进行分析，所述系统包括以下装置中的一种或多种：

(a)一种装置，其获取采集的参考样本的核酸的测序结果，在其SNP位点上检测每个样本的基因型，对于每一个SNP位点，参考等位基因记为A，交替等位基因记为B，三种基因型为AA、AB和BB，获得SNP位点上的参考等位基因和交替等位基因的测序深度d_A和d_B，计算得到这个SNP位点的测序深度d，其中d＝d_A+d_B，记其交替等位基因B占这个SNP位点的测序深度的比值为θ，即θ＝d_B/(d_A+d_B)，θ∈[0,1]；

参考样本的个数设为N，检测的SNP位点个数设为M，对第i个SNP，其中i为1，2，3，…，M，获得其参考等位基因和交替等位基因的测序深度d_iA和d_iB，计算得到这个SNP位点的测序深度d_i，其中d_i＝d_iA+d_iB，其交替等位基因B占测序深度的比值为θ_i＝d_iB/(d_iA+d_iB)，θ_i∈[0,1]；当θ_i靠近0、0.5和1时对应AA、AB和BB基因型；

以d为横坐标，θ为纵坐标做图，所述N个样本的每个SNP的d_i和θ_i值在图上形成相对于AA、AB和BB基因型的三个点簇，算出三个点簇的中心位置的D_g和Θ_g值，其中D_g和Θ_g值可为点簇上所有点的d_i或θ_i的平均值或中位数值：

D_g＝median(d_j,j∈{i:g_i＝g}),g＝AA,AB,BB；

Θ_g＝median(θ_j,j∈{i:g_i＝g}),g＝AA,AB,BB，

(b)一种装置，其通过线性插值法计算待测样本的每一个SNP在个体的正常状态下的测序深度d的期望值d_e；

(c)一种装置，其获取待测样本的核酸的测序结果，获得其每个SNP测序深度；

(d)一种装置，其计算得到待测样本的每个SNP位点的测序深度比，即SDR，其对数形式为LSDR；以及计算得到交替等位基因频率，即AAF，其中

SDR＝d/d_e；

LSDR＝log₂(d/d_e)；

以及

(e)一种装置，其采用隐马尔可夫模型分析上述步骤得到的数据，检查和展示目标基因组的拷贝数变异状况。

18.一种用于实现如权利要求9-16中任一项所述的方法的系统，用于对受测者的体液样品基因组中的拷贝数变异进行分析，其包括以下装置中的一种或多种：

(a)一种装置，其获得目标个体的体液游离核酸和含有正常基因组的组织样品的参考核酸的测序结果，在其SNP位点上检测其基因型，对于每一个SNP位点，参考等位基因记为A，变异等位基因记为B，三种基因型为AA、AB和BB；

(b)一种装置，其计算得到目标个体的每个杂合型SNP位点的测序深度比SDR，其对数形式为LSDR，其方法如下：

获得每个杂合型SNP位点附近2L宽度区域内，即SNP上游宽度为L和下游宽度为L；L范围为300-2000bp，每个碱基位置的测序深度，参考核酸和体液游离核酸在第j个碱基位置上的测序深度记为d_N,i,j和d_T,i,j；

检测的SNP位点个数设为M，对第i个杂合型SNP，其中i为1，2，3，…，M，参考核酸和体液游离核酸的加权测序深度分别为：

其中w(l_i,j)为权重；

第i个杂合型SNP位点的LSDR值计算为：

(c)一种装置，其计算得到受测者的每个杂合型SNP位点对应的交替等位单倍体频率，即AHF，其方法如下：

获得受测者的体液游离核酸中第i个杂合型SNP周围宽度为2R范围内的杂合基因型SNP位点，即上游宽度为R,下游宽度为R；R范围为10-100kb，的基因型；对基因型数据做单倍型定相，判断临近的杂合型SNP位点的两个等位基因A和B各处在哪个单倍体上；在血清游离DNA中，在A和B单倍体上第i个杂合型SNP位点附近其他杂合型SNP位点的测序深度记为d_T,A,i,j和d_T,B,i,j；

计算第i个杂合型SNP在A和B单倍体上的加权测序深度：

其中w(l_i,j)为权重；

血清游离DNA中SNP_i位点的AHF值为：

参考DNA中SNP_i位点的AHF值为：

进一步计算：

LmAHF_i＝log₂(mAHF_T,i/mAHF_N,i)

其中mAHF_N,i＝|AHF_N,i-0.5|+0.5，mAHF_T,i＝|AHF_T,i-0.5|+0.5；

(d)一种装置，其计算出M个杂合SNP位点的LmAHF值记为x和LSDR值记为y后，得到一组2维序列{(x_i,y_i),i＝1,...,M}；采用协同分段方法，来确定体液游离DNA中的CNV片段，即把{(x_i,y_i),i＝1,...,M}分为K个片段

{η_k＝{(x_i,y_i):i∈(t_k-1,t_k]}；k＝1,...,K}，0＝t₀<t₁<...<t_k-1<t_k<...<t_K＝M；

其中每个片段η_k中的杂合型SNP位点{i:i∈(t_k-1,t_k]}对应同一个CNV状态，相邻的片段对应不同的CNV状态；

(e)一种装置，其对于所述K个片段，推断其CNV状态，其中CNV状态为哪些片段属于正常的拷贝数状态，以及哪些CNV片段有拷贝数扩增，缺失或杂合性缺失。

Images