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

权利要求书

1.一种用于计算癌症样本中癌症细胞纯度和染色体倍性的方法，所述方法包括以下步骤：

步骤A：

获取配对的癌症组织样本和正常组织样本的全基因组测序数据，并将测序数据比对到参考基因组；

步骤B：

从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，read是指高通量测序平台产生的测序序列；HGSNV是指杂合生殖系细胞单碱基变异；MAF表示主要等位基因分数，指两种等位基因分数中的较大值，等位基因分数的计算方法为覆盖某一等位基因的read数量除以覆盖该位点总read数量的比值；MAF计算公式如(1.1)所示：

公式(1.1)中，n^r为包含与参考基因组相同等位基因的read数量，n^a为包含另一种等位基因的read的数量，n^t表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；

步骤C：

根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量，其中，window是指按照一定长度划分的基因组片段，该长度代表window大小；

步骤D：

使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC-seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段，其中，TRE表示癌症片段读长富集程度e_s，指癌症样本中某一片段s内read数量与相应正常样本中对应片段read数量的比值，定义公式如下：

公式(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)所示的类自回归模型，计算两个相邻拷贝数片段对应的peak的TRE的差值即P，其中，peak是指基因组所有window的TRE分布中，聚集在一起的TRE簇；以及其中，遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：

公式(12)和(13)中，X_t表示0到M_t之间的TRE值；t表示扩大了1000倍的TRE值；M_t表示TRE的最大值；变量P表示两个相邻拷贝数片段对应的peak的TRE的差值；C(X_t)表示在TRE为X_t的位点，对应的window数量；C(X_t+1000×P)表示在TRE为X_t+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；

步骤G：

根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际观测peak之前最多可能存在理论peak的数量N，最后当第一个实际观测peak之前存在的理论peak数量为n时，计算Q的值，以Q_n表示，其中步骤G包括：

G1：

根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的X_f作为第一个实际观测peak的TRE均值：

公式(13.1)中，i表示第i个peak，C(X_f+P×i)表示在TRE为X_f+P×i的位点，对应的window数量，n_m表示M_t以内peak的最大数量，M_t表示TRE的最大值；

G2：

使用公式(13.2)，根据步骤F计算的P和步骤G1计算的X_f，计算在X_f之前最多可能存在的peak数量N：

公式(13.2)中，X_f表示第一个实际观测peak的TRE均值，P表示两个相邻拷贝数片段对应的peak的TRE的差值，floor表示向下取整数；

G3：

利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Q_n的值：

Q_n＝X_f-n×P+2×P＝X_f+(2-n)×P，n∈[0，N] (13.3)

公式(13.3)中，n表示第一个实际观测peak之前存在的理论peak数量，取值范围是0到N之间的整数，P表示两个相邻拷贝数片段对应的peak的TRE的差值，X_f表示第一个实际观测peak的TRE均值，Q_n表示在X_f之前理论上存在n个peak时的Q值；

步骤H：

使用步骤F计算的P与步骤G计算的Q_n，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：

公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，由此对(P，Q_N)得到对应的(γ，κ)；

步骤I：

当n取[0，N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：

T_i＝X_f-n×P+i×P＝X_f+(i-n)×P，n∈[0，N] (13.4)

公式(13.4)中，n表示第一个实际观测peak之前存在的理论peak数量，取值范围是0到N之间的整数，P表示两个相邻拷贝数片段对应的peak的TRE的差值，X_f表示第一个实际观测peak的TRE均值，T_i表示第i个peak的TRE均值，

对于落在T_i附近的片段，认为该片段具有拷贝数i；对于没有落在T_i附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，计算peak的MAF的期望f_b，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{f_b}；同时计算各个peak的TRE均值和标准差；

步骤J：

根据步骤F计算的P和步骤I计算的{f_b}构建如公式(19)所示的用“贝叶斯信息准则”校正后的最终的混合高斯模型，然后对最终的混合高斯模型极大似然估计；其中，步骤J包括如下几步：

J1：

以步骤F计算的P构建如公式(17)所示的高斯分布模型：

公式(17)中，L(e_s；γ，κ)表示基因组片段TRE的似然函数，N_w表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σ_i表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，e_s为第s个window的TRE观测值，Sⁱ表示第i个peak的TRE均值即步骤I中的T_i，p_i表示第s个window的拷贝数为i的权重，对所有的i，p_i均取值为1；

J2：

以步骤I计算的f_b构建如公式(18)所示的第一混合高斯模型：

公式(18)中，L(f_s；γ，κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；F^i，j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；f_s表示该片段内所有HGSNV的MAF的均值，由步骤E得到；σ_i，j表示该片段内所有HGSNV的MAF的标准差，由步骤E得到；p_i，j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，p_i，j取值均为1，p_i表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，p_i取值均为1；

J3：

将(17)与(18)相加得到第二混合高斯模型，然后对第二混合高斯模型进行BIC(Bayesian Information Criterion)校正得到最终的混合高斯模型如公式(19)：

BIC(e_s，f_s；γ，κ)＝-2×logL(f_s；γ，κ)-2×logL(e_s；γ，κ)+I×log(N_w)+J×log(M) (19)

公式(19)中，BIC(e_s，f_s；γ，κ)表示最终的混合高斯模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N_w是基因组中window的数量，M是基因组中HGSNV的个数，

对[0，N]范围内的每一个整数值n，通过步骤G得到Q_n，或者通过步骤I得到所有peak的MAF期望的集合{f_b}，由一对(P，{f_b})构建一个公式(19)所示的模型；

步骤K：

以0.001为分辨率，对[P-m，P+m]区间的所有P值，重复步骤G～J，得到一系列不同的(P，Q_n)与对应的似然函数值，取最大的似然函数值对应的(P，Q_n)作为最合适的P和Q值，m是0到0.5之间的一个值；

步骤L：

查询步骤H的结果，找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。