[发明专利]用于计算癌症样本纯度和染色体倍性的方法和装置

摘要

本发明提供了一种全自动、高效率、高准确性的计算癌症样本纯度和染色体倍性的方法和装置。通过本发明提供的层次混合高斯模型，实现了对癌症样本纯度和染色体倍性的快速和准确计算，节约了纯度估算的时间和经济成本，同时提高了计算结果的准确性。本发明在癌症样本纯度和染色体倍性计算上具有广阔的运用前景。

主权项

1.一种用于计算癌症样本中癌症细胞纯度和染色体倍性的方法，所述方法包括以下步骤：步骤A：获取配对的癌症组织样本和正常组织样本的全基因组测序数据，并将测序数据比对到参考基因组；步骤B：从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，计算公式如(1.1)所示：公式(1.1)中，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；步骤C：根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量；步骤D：使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC‑seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段：公式(1)中，和分别表示在癌症样本中覆盖片段s的read数量和在正常样本中覆盖片段s的read数量，N_t表示癌症样本总read数量，N_n表示相应正常样本总read数量，e_s为TRE值；步骤E：以步骤D中BIC‑seq处理后的基因组片段为单位，统计片段内所有window的TRE的均值、方差和该片段内window数量，根据均值和方差对基因组每个片段的window数量进行平滑化处理，使TRE的分布更均匀，然后将平滑化处理后所有片段的window分布汇总，得到基因组上window随TRE变化的分布结果；同时以片段为单位，计算片段中所有HGSNV的MAF的均值和方差；步骤F：使用如公式(12)、(13)所示的类自回归模型，计算相邻拷贝数片段内TRE的差值即P，其中，遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；变量P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；步骤G：根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际peak之前最多可能存在理论peak的数量N，最后当第一个实际peak之前存在n个理论peak时，计算Q的值，以Qn表示，其中步骤G包括：G1：根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的Xf作为第一个实际观测peak的TRE均值：公式(13.1)中，i表示第i个peak，C(Xf+P×i)表示在TRE为Xf+P×i的位点，对应的window数量，n表示Mt以内peak的最大数量，Mt表示TRE的最大值；G2：使用公式(13.2)，根据步骤F计算的P和步骤G1计算的Xf，计算在Xf之前最多可能存在的peak数量N：公式(13.2)中，Xf表示第一个peak的均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数；G3：利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Qn的值：Qn＝Xf‑n×P+2×P＝Xf+(2‑n)×P，n∈[0，N]   (13.3)公式(13.3)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Qn表示在Xf之前理论上存在n个peak时的Q值；步骤H：使用步骤F计算的P与步骤G计算的Qn，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，由此对(P，QN)得到对应的(γ，κ)；步骤I：当n取[0,N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：Ti＝Xf‑n×P+i×P＝Xf+(i‑n)×P，n∈[0，N]   (13.4)公式(13.4)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Ti表示第i个peak的TRE均值，对于落在Ti附近的片段，认为该片段具有拷贝数i；对于没有落在Ti附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，计算peak的MAF的期望fb，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{fb}；同时计算各个peak的TRE均值和方差或标准差；步骤J：根据步骤F计算的P和步骤I计算的{fb}构建如公式(19)所示的用“贝叶斯信息准则”校正后的混合高斯分布模型，然后对模型极大似然估计；其中，步骤J包括如下几步：J1：以步骤F计算的P构建如公式(17)所示的高斯分布模型：公式(17)中，L(es；γ,κ)表示基因组片段TRE的似然函数，N表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σi表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，es为第s个window的TRE观测值，Si表示第i个peak的TRE均值即步骤I中的Ti，pi表示第s个window的拷贝数为i的权重，对所有的i，pi均取值为1；J2：以步骤I计算的fb构建如公式(18)所示的高斯分布模型：公式(18)中，L(fs；γ,κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；Fi,j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；fs表示该片段内所有HGSNV的MAF的观测值均值，由步骤E得到；σi,j表示该片段内所有HGSNV的MAF观测值的标准差，由步骤E得到；pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1，pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1；J3：将(17)与(18)相加得到混合高斯模型，然后对混合模型进行BIC(Bayesian Information Criterion)校正得到最终混合模型如公式(19)：BIC(es,fs；γ,κ)＝‑2×log L(fs；γ,κ)‑2×log L(es；γ,κ)+I×log(N)+J×log(M)  (19)公式(19)中，BIC(es,fs；γ,κ)表示混合模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N是基因组中window的数量，M是基因组中HGSNV的个数，对[0,N]范围内的每一个整数值n，通过步骤G得到Qn，或者通过步骤I得到所有peak的MAF期望的集合{fb}，由一对(P，{fb})构建一个公式(19)所示的模型；步骤K：以0.001为分辨率，对[P‑m,P+m]区间的所有P值，重复步骤G～J，得到一系列不同的(P，Qn)与对应的似然函数值，取最大的似然函数值对应的(P，Qn)作为最合适的P和Q值，m是0到0.5之间的一个值；步骤L：查询步骤H的结果，找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。