阅读说明：本技术 互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法 (Method for identifying heterogeneous cancer driver genes of mutually exclusive constraint graph Laplace ) 是由 习佳宁 黄庆华 于 2020-06-23 设计创作，主要内容包括：本发明提供了一种互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法。首先,获取癌症基因组变异数据与基因互作关系网络；然后,采用矩阵化模型对癌症的异质性进行描述,并通过互斥性约束矩阵分解对异质性癌症的样本参数进行差异化估计；接着,构建联合关联互作网络正则化的互斥性约束矩阵分解优化函数,并通过迭代求解对局部样本中的受互作影响驱动基因参数进行修正；最后,采用离群值检验方法识别驱动基因。本发明能够解决癌症样本的参数差异化估计与受互作影响局部样本驱动基因的有效识别问题,实现从异质性癌症样本的基因变异数据中识别出在局部样本发生变异的驱动基因。(The invention provides a heterogeneous cancer driver gene identification method based on mutual exclusion constraint graph Laplace. Firstly, acquiring a cancer genome variation data and gene interaction relation network; then, describing heterogeneity of the cancer by adopting a matrix model, and performing differential estimation on sample parameters of the heterogeneous cancer through mutual exclusion constraint matrix decomposition; then, constructing a mutual exclusion constraint matrix decomposition optimization function for joint correlation interaction network regularization, and correcting the parameters of the interaction-influenced driving genes in the local samples through iterative solution; finally, the driver genes were identified using outlier testing. The method can solve the problems of parameter differentiation estimation of the cancer sample and effective identification of the local sample driver gene influenced by interaction, and realizes identification of the driver gene with variation in the local sample from gene variation data of the heterogeneous cancer sample.)

互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法

技术领域

本发明属生物信息学、基因组数据挖掘技术领域，具体涉及一种互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法。

背景技术

癌症是一种高发恶性疾病，主要由驱动基因的变异所导致。但在癌症基因组中，存在大量与癌变无关的伴随性变异，对判别驱动基因造成严重混淆。由于驱动基因变异比伴随性变异更倾向于在多例样本中同时发生，现有研究主要通过癌症样本的基因变异数据，将驱动基因视为多例样本中的高频变异基因，并通过考察基因变异率的统计显著性，寻找多样本中显著高频变异的驱动基因。如Lawrence等人在文献"Lawrence S,Stojanov P,Polak P,et al.Mutational heterogeneity in cancer and the search for newcancer-associated genes[J].Nature,2013,499(7457):214"中提出了一种基于各基因背景变异率修正的变异频率统计检验方法，以识别癌症样本中显著高频变异的基因。Kumar等人则在文献"Kumar D,Swamidass J,Bose R.Unsupervised detection of cancer drivermutations with parsimony-guided learning[J].Nature genetics,2016,48(10):1288"中对驱动基因预测进一步采用精简化约束，以减少识别中的假阳性结果。然而，由于癌症还具有肿瘤异质性，即不同样本的变异驱动基因存在较大差异，加剧了驱动基因在局部样本分布的复杂性。对于具有肿瘤异质性的癌症，由于不同样本中的变异驱动基因差别较大，导致仅在局部样本中的变异驱动基因相对于全体样本呈现的变异率较低。尽管现有研究在异质性癌症各样本的亚类已知时，可对各类样本内的高变异率基因分别进行识别，但在样本类型缺失时则无法区分具有差异性的各局部样本，因而无法识别异质性癌症的局部样本驱动基因。

由于驱动基因还可受其他变异基因的互作影响而发生异常，相关研究还将变异率作为各基因的功能异常影响力，通过基因互作关系传播进行建模，以筛选结果中受影响程度较高的驱动基因。如Raphael等人在文献"Leiserson D,Vandin F,Wu H,et al.Pan-cancer network analysis identifies combinations of rare somatic mutationsacross pathways and protein complexes[J].Nature genetics,2015,47(2):106"中将基因变异频率作为影响力，通过基因相互作用关系对影响力进行传播，并将传播后得分作为各基因受相互作用的影响程度，以识别癌症样本共有的关联性驱动基因。由于上述传播会对大量未变异基因造成影响，进而导致预测结果的假阳性问题，对此Cho等人在文献"ChoA,Shim E,Kim E,et al.MUFFINN:cancer gene discovery via network analysis ofsomatic mutation data[J].Genome biology,2016,17(1):129"中仅允许高频变异基因对直接相互作用基因产生影响，通过限制传播过程以避免多级传播所造成的识别误判。为在统计显著性层面过滤无关基因，Horn等人在文献"Horn H,Lawrence S,Chouinard R,etal.NetSig:network-based discovery from cancer genomes[J].Nature methods,2018,15(1):61"中通过相互作用关系对高频变异基因影响力进行显著性描述，进而提升共有关联性驱动基因的预测率。然而，上述基于互作关系传播的研究仅考虑了基因变异率在全体样本层面的影响，仍然无法考察基因互作关系对局部样本造成的影响。对于异质性癌症，由于在建模中基因互作关系对局部样本的影响缺失，导致受基因互作影响的驱动基因在局部样本漏检。

综上，目前研究存在以下问题：1)异质性癌症的样本类型缺失情况时，发生变异的驱动基因在不同样本差异较大，导致其在全体样本间的变异率较低而难以有效识别；2)现有方法均以基因在全体样本的变异频率作为标准，但在异质性癌症中，受基因互作影响的驱动基因会在局部样本发生漏检。因此，由于缺少对样本类型缺失时的癌症样本差异化方法思路，对受互作影响的局部样本驱动基因漏检问题认识不足，制约了异质性癌症的致病机理与临床诊疗研究发展。

发明内容

为了克服现有技术的不足，本发明提供一种互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法。首先，获取癌症基因组变异数据与基因互作关系网络；然后，采用矩阵化模型对癌症的异质性进行描述，并通过互斥性约束矩阵分解对异质性癌症的样本参数进行差异化估计；接着，构建联合关联互作网络正则化的互斥性约束矩阵分解优化函数，并通过迭代求解对局部样本中的受互作影响驱动基因参数进行修正；最后，采用离群值检验方法识别驱动基因。本发明能够解决癌症样本的参数差异化估计与受互作影响局部样本驱动基因的有效识别问题，实现从异质性癌症样本的基因变异数据中识别出在局部样本发生变异的驱动基因。

一种互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法，其特征在于步骤如下：

步骤1：获取癌症基因组变异数据与基因互作关系网络，通过基因注释数据库对基因名进行ID统一化；

步骤2：根据所输入的癌症样本中各个基因是否发生变异的情况，将各样本构建为0/1向量，向量的元素为各基因是否变异的0/1取值，并拼接所有0/1向量构成变异矩阵X＝[x_ij]_m×n，其中，m为样本个数，n为基因个数，i＝1,2,…,m，j＝1,2,…,n，x_ij为矩阵X中第i行第j列元素，若第i个样本的第j个基因存在变异则x_ij值为1，否则值为0；

步骤3：对所输入的m个样本和n个基因，设置空矩阵U＝[u_ik]_m×r和V＝[v_jk]_n×r，其中，r为给定的参数维数，满足r＜＜m,n，矩阵U记为样本参数矩阵，其每一个元素记为样本参数，u_ik表示第i个样本在r元空间的第k个参数，矩阵V记为基因参数矩阵，其每一个元素记为基因参数，v_jk表示第j个样本在r元空间的第k个参数，k为r元空间的各维度参数的索引，k＝1,2,…,r；

步骤4：求解如下互斥性约束的矩阵分解模型，得到矩阵U和V中的样本参数和基因参数：

其中，l表示当前正在对r元空间的第l个参数进行考察，l＝1,2,…,r，ε为可调阈值，取值范围为(0,0.1]；

步骤5：按下式计算得到基因参数的图拉普拉斯正则化项Reg_Y(V)：

其中，Y表示基因之间相互作用关系所组成的集合，Y＝{(i,j)|i与j存在互作关系}，i和j分别代表第i个基因和第j个基因，(i,j)表示两个存在相互关系的基因所组成的基因对；v_i为矩阵V中的第i行，v_j为矩阵V中的第j行；I为示性函数，I_Y(i,j)为示性函数中的元素，如果基因互作网络中的基因i与基因j存在邻边，则I_Y(i,j)取值为1，否则I_Y(i,j)取值为0；

迭代求解以下联合图拉普拉斯正则化项的矩阵分解优化函数，得到融合基因互作网络的样本参数矩阵U′和基因参数矩阵V′：

其中，λ表示正则化项的调谐参数，取值为大于零的实数；

步骤6：对于步骤5所得的矩阵U′，将其第k列元素最大值的索引所对应的样本作为第k组亚群的局部样本，k＝1,2,…,r，得到所有r个组亚群的局部样本；

步骤7：对r组亚群的具备样本，分别采用离群值检验方法进行驱动基因检测，得到每组局部样本所对应的原假设分布，具体为：

首先，对于r组亚群中的第k组亚群，k＝1,2,…,r，选取第k组亚群的局部样本在变异矩阵X中所对应的行，构成第k组亚群所对应的第k个子矩阵，并采用重启动的随机游走算法，利用基因互作关系集合Y对第k个子矩阵进行随机游走处理；然后，对随机游走得到的游走后矩阵进行随机重排，再将重排后的矩阵中的所有1×n维的行向量进行相加，加和后的1×n维向量为n个基因的分布水平采样,向量中n个维度的取值代表n个基因在本次采样的分布水平，将前述随机重排采样重复10000次，得到10000次n个基因的采样结果，将n个基因的10000次采样结果取值，对n个基因分布构建n个取值频次分布图，作为第k组亚群的局部样本的n个基因所分别对应的n个原假设分布；

步骤8：在第k组局部样本的n个原假设分布中，第j个基因所对应的原价设分布为n个原假设分布中的第j个原假设分布，将基因参数矩阵V′中的元素v_j′_k值，与第j个原假设分布的横坐标进行比对，将分布函数比对位置的右侧面积作为检验p值，通过Benjamini-Hochberg错误发现率矫正算法对检验p值进行错误发现率矫正，得到矫正后的p值，如果矫正后的p值小于0.05，则认为第j个基因是第k组局部样本的驱动基因，令j的取值从1取到n，k的取值从1取到r，按前述过程进行处理，则可以得到每一个基因是否为所有局部样本的驱动基因的识别结果。

本发明的有益效果是：1)针对异质性癌症的样本在基因组变异层面存在的较大差异，在样本类型缺失情况下，通过对癌症样本参数进行差异化估计，解决异质性癌症的样本难以区分问题。2)针对受基因互作影响驱动基因在局部样本的漏检问题，通过将基因互作关系对驱动基因参数进行正则化，实现基因参数在互作影响层面的关联性桥接，提升对局部样本驱动基因识别性能。

附图说明

图1是本发明的互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法流程图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

如图1所示，本发明提供了一种互斥性约束图拉普拉斯的异质性癌症驱动基因识别方法，其具体实现过程如下：

步骤1：通过The Cancer Genome Atlas(TCGA)和International Cancer GenomeConsortium(ICGC)等数据库获取癌症基因变异数据。从STRING相互作用数据库收集基因互作关系网络的数据。通过基因注释数据库Database for Annotation,Visualization andIntegrated Discovery(DAVID)对基因名进行ID统一化，以消除不同来源数据中基因的同物异名现象。

步骤2：采用矩阵化模型对癌症的异质性进行描述，矩阵的行与列分别描述差异性局部样本及其相应驱动基因。通过多元参数对不同局部样本进行描述，参数取值的差异代表各局部样本间存在的区别，以有效量化样本的差异性。具体为：根据所输入的癌症样本中各个基因是否发生变异的情况，将各样本构建为0/1向量，，向量的元素为各基因是否变异的0/1取值，并拼接所有0/1向量构成变异矩阵X＝[x_ij]_m×n，其中，m为样本个数，n为基因个数，i＝1,2,…,m，j＝1,2,…,n，x_ij为矩阵X中第i行j列元素，若第i个样本的第j个基因存在变异则x_ij值为1，否则值为0。

步骤3：采用多元参数描述基因在局部样本的异常程度，以反映异质性癌症各局部样本的驱动基因，记为样本参数和基因参数。即对所输入的m个样本和n个基因，设置空矩阵U＝[u_ik]_m×r和V＝[v_jk]_n×r，其中，r为用户自定的参数维数，在给定时需满足r＜＜m,n，矩阵U记为样本参数矩阵，其每一个元素记为样本参数，u_ik表示第i个样本在r元空间的第k个参数，矩阵V记为基因参数矩阵，其每一个元素记为基因参数，v_jk表示第j个样本在r元空间的第k个参数,k＝1,2,…,r。

步骤4：采用互斥性约束矩阵分解对变异数据的低维参数进行估计，以保障不同局部样本的差异性：

其中，l表示当前正在对r元空间的第l个参数进行考察，k为r元空间的各维度参数的索引，ε为一个用户可调的阈值，可取值范围为(0,0.1]。

上式通过将协方差限制在较小取值的方式实现互斥性约束，从而突显基因组变异区别较大的不同局部样本，以对异质性癌症的样本参数进行差异化估计。在估计结果中，样本参数与基因参数通过矩阵乘法表示可近似重建基因变异的原始数据矩阵，其中各行描述基因变异层面的样本间区别，各列则表示各基因在局部样本的异常程度。

步骤5：将基因互作网络中各基因是否存在邻边的关系，通过示性函数I进行表示，即基因i与基因j如果存在邻边则I_Y(i,j)取值为1，否则I_Y(i,j)取值为0。根据基因在互作网络中的关联关系，采用图拉普拉斯方法对各基因参数的距离进行正则化：

其中，Reg_Y(V)即为基因参数的图拉普拉斯正则化项。Y表示基因之间相互作用关系所组成的集合，即Y＝{(i,j)|i与j存在互作关系}，此处i和j分别代表第i个和第j个基因，(i,j)表示两个存在相互关系的基因所组成的基因对。v_i为矩阵V中的第i行，v_j为矩阵V中的第j行。

上式可反映在各局部样本的基因参数关于互作网络的关联程度。通过将关联性桥接的正则化项，与样本参数和基因参数进行联合估计，构建以下联合关联互作网络正则化的互斥性约束矩阵分解优化函数：

其中，λ表示正则化项的调谐参数，由用户自定，取值为大于零的实数。

通过对上式进行迭代求解，存在互作关系的局部样本基因参数距离逐渐邻近，从而对局部样本中的受互作影响驱动基因参数进行修正。

步骤6：上述最终所得的用于反映局部样本指示关系的样本参数矩阵U′，将其第k列元素最大值的索引所对应的样本作为第k组亚群的局部样本，k＝1,2,…,r，得到所有r个组亚群的局部样本；

步骤7：对r组亚群的具备样本，分别采用离群值检验方法进行驱动基因检测，得到每组局部样本所对应的原假设分布，具体为：

首先，对于r组亚群中的第k组亚群，选取第k组亚群的局部样本在变异矩阵X中所对应的行，重新构成第k组亚群所对应的第k个子矩阵。采用重启动的随机游走算法，利用基因互作关系集合Y对第k个子矩阵进行随机游走处理，然后，对随机游走得到的游走后矩阵进行随机重排，再将重排后的矩阵中的所有1×n维的行向量进行相加，加和后的1×n维向量则为n个基因的分布水平采样,向量中n个维度的取值代表这n个基因再本次采样的分布水平。将上述随机重排采样重复10000次，得到10000次全体n个基因的采样结果，将n个基因的10000次采样结果取值，对n个基因分布构建n个取值频次分布图，作为第k组亚群的局部样本的n个基因所分别对应的n个原假设分布

步骤8：在第k组局部样本的n个原假设分布中，第j个基因所对应的原价设分布为n个原假设分布中的第j个原假设分布。将基因参数矩阵V′中的元素v_j′_k值，与第j个原假设分布的横坐标进行比对，将分布函数比对位置的右侧面积作为检验p值。通过Benjamini-Hochberg错误发现率矫正算法对检验p值进行错误发现率矫正，得到矫正后的p值，如果矫正后的p值小于0.05，则认为第j个基因是第k组局部样本的驱动基因，令j的取值从1取到n，k的取值从1取到r，按上述过程进行处理，则可以得到每一个基因是否为所有局部样本的驱动基因的识别结果。