本发明属于计算生物学领域，涉及一种基于线性模型的基因共调控网络中关键调控子识别方法。
背景技术：
：在后基因组时代，理解基因、非编码RNA、蛋白质和其他相关生物分子的功能，提示生物过程的实现机理成为当前计算系统生物学和生物信息学最重要的研究目标之一。其中，基因调控的研究是一个非常重要的课题。了解基因表达的调控机制，对于我们认识生物学过程和疾病的发生机制都起到了重要的作用。在真核生物中，有两类重要的调控因子：转录因子(TF)和microRNA(miRNA)，分别在转录水平及转录后水平上调控靶基因的表达水平。转录因子是一类具有特定功能的蛋白质，它通过结合到基因的启动子区域来开启基因的转录过程。miRNA是近年来研究发现的一种新的基因调控元件，是在真核生物中发现的一类内源性的具有调控功能的非编码RNA，其大小长约20-25个核苷酸。转录因子、miRNA在基因表达调控中起到了重要的作用，这种调控作用遍及各种生物活动以及疾病发生过程。在此基础上，研究发现转录因子和miRNA存在着广泛的相互作用和合作调控，它们组成了一个复杂的共调控网络。共调控网络中包含转录因子调控miRNA、转录因子调控靶基因、miRNA调控转录因子和靶基因的调控作用，这些调控作用体现了细胞分子生命过程和功能执行的各个阶段，所以共调控网络包含比单个网络更丰富的生物信息。因此，有效地识别共调控网络上的关键调控子对疾病的临床治疗和药物设计都有重要作用，这将可能会给人类疾病的治疗提供一种新的手段。随着高通量技术的迅速发展，产生了大量的基因组学、转录组学和蛋白质组学等组学数据，为生物分子功能研究提供了新的契机。以往对于关键点的识别算法，主要是集中在蛋白质相互作用网络上对于关键蛋白的识别。与蛋白质相互作用网络相比，转录调控网络的进化研究则更困难。首先，可信的转录调控网络数据仍然不容易获取；其次，就已有的转录调控网络来看，由于网络本身的功能特性，其展现出的拓扑特性与蛋白质相互作用网络有较大不同，且加上调控作用的有向性，使得调控网络展现出的拓扑特性更加复杂。因此对于调控网络上关键调控子的识别也比关键蛋白的识别更为复杂。近年来，对于调控网络的研究越来越多，已经有了多种基于计算的方法来识别调控网络上的关键调控子，主要是有下面这几类方法：基于信息流模型(RWR)、排名算法(PageRanking)、构建分类器(SVM、Regularizedleast-squaresclassification)、贝叶斯网络、基于回归模型等等。然而，已有的方法或多或少存在一些问题：比如不能处理大数据、时间复杂度太高、精度有待提高等等。2015年，Alexandra等提出了MIPRIP方法，利用线性模型来识别调控网络上的关键调控子，实验结果表明，基于线性模型的方法能有效的识别具有重要生物意义的调控子。然而，该方法只是单纯的考虑了转录因子与基因之间的关系，并没有考虑到共调控网络中调控子之间的相互作用和合作调控关系，同时识别精度也有待提高。因此，有必要设计一种基于线性模型的基因共调控网络中关键调控子识别方法。技术实现要素：本发明所要解决的技术问题是提供一种基于线性模型的基因共调控网络中关键调控子识别方法。该基于线性模型的基因共调控网络中关键调控子识别方法只需根据基因表达谱数据和基因调控关系就能比较准确地识别基因共调控网络中具有生物意义的关键调控子。发明的技术解决方案如下：一种基于线性模型的基因共调控网络中关键调控子识别方法，包括以下步骤：步骤1)构建基因共调控网络：输入基因表达谱数据、基因调控关系和蛋白质相互作用数据(Protein-ProteinInteraction，PPI)，过滤掉其中存在无表达谱数据节点的作用关系对，建立基因共调控网络GCN(geneco-regulatorynetworks)，基因共调控网络GCN中共包含三种节点：调控子miRNA(microRNA)、调控子TF和基因gene，节点之间存在作用边：miRNA-gene、TF-gene和gene-gene；若基因共调控网络GCN中任意两点，存在作用关系则边权为1，否则为0；步骤2)对已知疾病基因，分别计算调控子miRNA、调控子TF和邻接基因的活动值；活动值即miRNA、TF和邻接基因对已知疾病基因的影响值；步骤3)在已构建的基因共调控网络GCN中，利用基因表达谱数据以及步骤2)中得到的调控子和邻接基因的活动值构建线性模型，预测已知疾病基因的表达，获取已知疾病基因的预测表达值；步骤4)根据已知疾病基因的预测表达值和真实表达值之间的差值最小化将步骤3)构建的线性模型转化为最优化问题，基于混合整数线性规划思想对最优化问题进行求解，最终识别基因共调控网络中的关键调控子。进一步地，所述用于预测已知疾病基因的表达构建的线性模型表达式如下：其中，i表示已知的疾病基因，m、t、g分别表示调控子miRNA、调控子TF、已知的疾病基因i的邻接基因；g′i,s表示已知的疾病基因i在样本s中预测的表达值，β0是指线性模型的附加权重(additiveoffset)，M、T、G分别表示miRNA集、TF集、gene集；βm、βt、βg分别表示m、t、g的优化参数，在步骤4)中的最优化问题处理时利用优化器可以直接计算得出；esm,i、tst,i、gsg,i分别表示m、t、g与i的作用边权值，取值为0或1；actm,s、actt,s、actg,s分别表示m、t、g在样本s中的活动值；所述样本s是指已知疾病的某个观测个体的数据。进一步地，所述的根据基因预测表达值和真实表达值之间的差值最小化将线性模型转化为最优化问题，表示为：其中，gi,s、g′i,s分别表示疾病基因i在样本s中真实的表达值和预测的表达值，O与S分别表示的是已知的疾病基因集和该疾病的总样本集；采用Gurobi优化器对所述优化问题进行求解，记录每个调控子在求解优化问题过程中被优化器选择的次数，根据选择次数对所有调控子进行排名，取排名前50的调控子作为最终的候选调控子。当安装好Gurobi优化器后，只需在R语言中导入gurobi包，便可直接调用gurobi函数进行优化问题处理，该gurobi函数存在三个输入参数：优化模型、timeLimit和OutputFlag，timeLimit一般取值600，OutputFlag取默认值0，所述优化模型即为构建的线性模型利用已知疾病基因的预测表达值和真实表达值之间的差值最小化将构建的线性模型转化为最优化问题而得到的。为了得到一系列典型的不同大小的模型，通过约束基因的调控子个数来构建线性模型。对于每个已知疾病基因，分别设定调控子个数为1至k来构建线性模型。进一步地，所述的调控子miRNA、调控子TF和邻接基因的活动值分别由下面两种方法进行计算：1)计算调控子miRNA和调控子TF的活动值：第一步，首先计算调控子r的所有目标基因的基准表达值：其中，r表示调控子，为调控子miRNA或者调控子TF；表示调控子r的目标基因gt的基准表达值，取值为基因gt在调控子r表达水平趋于0的所有样本中表达值的平均值；e(r)->0表示调控子r表达水平趋于0；目标基因的基准表达值是指无调控作用影响时的目标基因的表达值；第二步，计算目标基因基准表达值与调控子影响后的真实表达值之间的差值，即目标基因的表达水平变化值有：其中，ygt,s表示调控子r的目标基因gt在样本s中的真实表达值，表示调控子r的目标基因gt的表达水平变化值；第三步，根据目标基因的表达水平变化值构建简单线性模型，求解调控子的活动值actr,s：其中，G’表示调控子r的目标基因集，分别表示调控子r的目标基因集的表达水平变化值总和与基准表达值总和；3)计算邻接基因的活动值，采用基于邻接基因对其所有作用基因的表达影响累积效应来求解，即：其中，N表示样本s中的基因总数，gsg,i表示基因g与样本s中的基因i的作用边权，gi,s表示样本s中的基因i在样本s中的表达值，所述样本s是指已知疾病的某个观测个体的数据。进一步地，对所述步骤2)获得的调控子和邻接基因的活动值进行归一化处理后，再用于步骤3)中的线性模型的构建。有益效果本发明提供了一种基于线性模型的基因共调控网络中关键调控子识别方法(co-BOTLM)，利用基因表达谱数据和基因调控关系，通过构建线性模型预测已知疾病基因的表达来完成基因共调控网络中关键调控子的识别。与已有基于线性模型识别关键调控子的方法相比，本发明所述的co-BOTLM方法具有以下优势：1)应用于共调控网络，共调控网络包含比单个网络更丰富的生物信息，因此所识别的调控子可能具有更重要的生物意义；2)加入蛋白质相互作用数据(PPI信息)，考虑基因的表达可能受到邻接基因的影响；3)引用新的方法计算调控子和邻接基因的活动值，有效提高了癌症基因表达预测的精度。本发明实现简单，只需根据基因表达谱数据和基因调控关系就能比较准确地识别基因共调控网络中的关键调控子。通过实验证明，本发明所述的方法co-BOTLM能有效的识别基因共调控网络中的关键调控子，且识别的关键调控子都具有很重要的生物意义。同时，通过对比其他方法，精确度也有所提高。具体实验结果图对比和分析详见实施例。附图说明图1是本发明co-BOTLM的流程图。具体实施方式以下将结合附图和具体实施例对本发明做进一步详细说明：实施例1：一、基于线性模型的基因共调控网络中关键调控子识别模型本发明将基因共调控网络中的关键调控子定义为：利用基因表达谱数据和基因调控关系，通过构建线性模型预测已知疾病基因的表达，从而识别出的在共调控网络中严重影响疾病基因表达的调控子。为了清晰描述基于线性模型的基因共调控网络中关键调控子识别模型，发明人将该模型的相关定义如下：提出的构建线性模型预测已知疾病基因表达，其表达形式如下：基于线性模型的基因共调控网络中关键调控子识别模型的目标是识别在共调控网络中严重影响疾病基因表达的调控子。通过利用基因表达谱数据和基因调控关系构建线性模型预测已知疾病基因的表达，来完成基因共调控网络中关键节点的识别。基于线性模型的基因共调控网络中关键调控子识别方法的整个流程如图1所示。首先输入基因表达谱数据、基因调控关系和PPI数据。方法co-BOTLM可以划分为4个子过程：1)构建基因共调控网络；2)考虑到基因的表达可能受到调控子和邻接基因的影响，因此针对已知疾病基因，分别计算miRNA、TF和邻接基因的活动值(即miRNA、TF和邻接基因对已知疾病基因的影响值)；3)在已得到的基因共调控网络中利用基因的表达谱数据构建线性模型，预测已知疾病基因的表达；4)根据基因预测表达值和真实表达值之间的差值最小化将线性模型转化为最优化问题，并基于混合整数线性规划思想(MILP)进行求解，最终识别基因共调控网络中的关键调控子，整个识别过程结束；采用Gurobi优化器对所述优化问题进行求解，记录每个调控子在求解优化问题过程中被优化器选择的次数，根据选择次数对所有调控子进行排名，取排名前50的调控子作为最终的候选调控子。当安装好Gurobi优化器后，只需在R语言中导入gurobi包，便可直接调用gurobi函数进行优化问题处理，该gurobi函数存在三个输入参数：优化模型、timeLimit和OutputFlag，timeLimit一般取值600，OutputFlag取默认值0，所述优化模型即为利用已知疾病基因的预测表达值和真实表达值之间的差值最小化将构建的线性模型转化为最优化问题而得到的。为了得到一系列典型的不同大小的模型，通过约束基因的调控子个数来构建线性模型。对于每个已知疾病基因，分别设定调控子个数为1至k来构建线性模型。在本实例中，k值为5(经多次实验，当k取值为5时，实验效果达到最佳)。二、基于线性模型的基因共调控网络中关键调控子识别方法有效性验证为了验证方法co-BOTLM的有效性，将co-BOTLM方法应用于一组卵巢癌数据集上。实验数据集包括：卵巢癌样本数据、基因调控关系、PPI数据、已知的卵巢癌有关的疾病基因。卵巢癌样本数据下载自TCGA数据库，共385个样本，通过过滤掉表达值绝对值过小或在各样本中无明显差异表达的基因，最后得到一个共385个样本，其中包含559个miRNA和12456个基因的卵巢癌表达谱数据集。作用关系数据包括miRNA-gene、TF-gene和PPI数据，分别下载自MicroCosmwebsite、ENCODE数据库和BioGrid数据库。通过将卵巢癌表达谱数据集和作用关系相互映射，最终构建了一个miRNA-TF基因共调控网络，网络中共包含三种类型的节点：12381个基因、559个miRNA和75个TF，节点之间存在的作用关系：59660对gene-gene、241722对miRNA-gene和9877对TF-gene。对于已知的卵巢癌有关疾病基因，从DDOC数据库下载共得到379个，过滤掉无表达谱数据或无调控作用关系的疾病基因，最终剩余123个。在本实例中进行了三折交叉验证实验，将方法co-BOTLM与Alexandra等人提出的MIPRIP方法在预测精度方面进行了对比，引用皮尔逊相关系数PCC来计算co-BOTLM方法预测的疾病基因表达数据和真实表达数据之间的相似性，PCC值越大，则相似性越高，进而表明co-BOTLM方法构建的线性模型准确度越高，因此实验结果的精度也更高。实例中PCC值使用R语言的cor函数进行计算。同时，在本实例中，还对co-BOTLM方法所识别出的调控子进行了特性及功能富集分析。1.实验结果分析，验证算法有效性表1：miRNA-TF基因共调控网络中排名前20的调控子No.识别的关键调控子目标基因个数优化器选择次数1hsa-mir-106a*377502hsa-mir-586508433hsa-mir-423-5p496384hsa-mir-515-3p512345hsa-mir-181a-2*496346hsa-mir-768-3p530327hsa-mir-663480328hsa-mir-539382319hsa-mir-2064773010hsa-mir-509-3p5523011hsa-mir-362-3p5122512hsa-mir-378*5192413hsa-mir-520c-3p5662414hsa-mir-33a5232415hsa-mir-29a*4952316hsa-mir-193a-3p4962317hsa-mir-6014842318FOXA21692319hsa-mir-26b4662220hsa-mir-30b54122在本实例中，三折交叉验证实验后，最终取得平均PPC值为0.535，表明本发明中的线性模型预测的基因表达值和真实表达值有比较高的相似性，因此证明了co-BOTLM方法构建的线性模型准确度比较高，能够有效识别网络中的关键调控子。实验运行完成后，根据优化器对所有调控子的选择次数，对其进行排名，取前50个作为本实例中的候选关键调控子。在上面表1中，列出了排名前20的调控子，可以看出，除FOXA2以外的任意调控子所调控的基因均不少于300个，并且其中很多基因已被证实和卵巢癌有关。由于TF实验数据过少，FOXA2的目标基因偏少。由此表明被识别的调控子与卵巢癌基因共调控网络中大量的基因存在作用关系，可能与大量基因(包括已知的卵巢癌疾病基因)的表达有关，因此在此共调控网络中具有至关重要的作用。2.方法co-BOTLM与MIPRIP方法实验对比，验证算法精确性表2：方法MIPRIP实验结果的PCC值No.1234510.33299070.43121500.44364490.47317760.489345820.31952370.42214950.45000000.46878500.485140230.32140190.43411210.45710280.47682240.4916822注：1-3：表示三折交叉验证实验，1-5：表示构建线性模型的调控子个数k值表3：方法co-BOTLM实验结果的PCC值No.1234510.50187500.57098210.59401790.61125000.622767920.48580360.55758930.58696430.60258930.616428630.49562500.55187500.56919640.59187500.6059821MIPRIP方法与本发明的co-BOTLM方法均是基于线性模型来识别特定疾病的关键调控子，然而存在三个不同之处：1)MIPRIP方法应用于调控网络，co-BOTLM方法应用于共调控网络，转录因子和miRNA存在着广泛的相互作用和合作调控，因此共调控网络包含比单个网络更丰富的生物信息；2)对于影响疾病基因表达的因素，除转录因子和miRNA之外，co-BOTLM方法还考虑了邻接基因对其可能产生的影响；3)MIPRIP方法与co-BOTLM方法的转录因子和miRNA活动值计算方式不同。由于MIPRIP方法应用于调控网络，不考虑网络中的共调控关系，因此本实例在进行对比实验室时，将转录因子看成是普通基因。表2、表3分别是方法MIPRIP与方法co-BOTLM实验结果取得的PCC值，从表中可以很明显看出，co-BOTLM方法取得了更高的PCC值，平均PCC值为0.571，而MIPRIP方法的平均PCC值为0.433。很明显，方法co-BOTLM预测的基因表达值与真实表达值有更高的相似性，因此，实验间接表明方法co-BOTLM精度性更高，所识别的关键调控子可靠性较高。3.实验结果功能富集分析，验证结果的有效性表4：排名前10的调控子GO富集分析Ncellularcomponentassemblycellularcomponentassemblyo.：调控子排名，富集的GO术语：按P-value(越小越好)排名前3的GO术语，GO个数：P-value<0.05的GO术语个数，P-value：<0.05表明富集程度高。表5：排名前10的调控子KEGG通路富集分析No.：调控子排名，富集的KEGG通路：按P-value(越小越好)排名前3的KEGG通路，KEGG个数：P-value<0.05的KEGG个数，P-value：<0.05表明富集程度高。为了验证本发明中的co-BOTLM方法所识别的关键调控子是具有生物意义的，在本次实例中，分别使用R语言的GOstats对所识别的关键调控子进行了GO富集分析和KEGG通路富集分析。表4和表5分别显示的是排名前10的调控子的GO和KEGG通路富集分析结果。很明显，从表4可以看出，本发明中的co-BOTLM方法所识别的排名前10的调控子大部分均富集了300个以上的GO术语，其中比较频繁富集的GO术语有：cellularcomponentorganization、cellularprocess、celldeath、negativeregulationofdendriticcelldifferentiation等等，表明了所识别的调控子大量参与了细胞有关的生命活动过程。hsa-mir-515-3p与hsa-mir-768-3p所富集的GO术语个数少于100，原因可能是由于此两个miRNA的目标基因与GOstats库中匹配基因较少，同时，Jiang等人在2016年已证明，由于hsa-mir-768-3p下调与黑色素瘤细胞的蛋白质合成中的MEK/ERK-mediated加强有关，因此hsa-mir-768-3p有可能在卵巢癌中有着潜在的预后功能。同理，从表5中很明显可以看出，排名前10的调控子大部分均富集了至少5条以上的KEGG通路，其中比较频繁富集的生物过程有：Prostatecancer、pathwaysincancer、signalingpathway、ErbBsignalingpathway等等，表明了所识别的调控子参与了大量的癌症及信号通路，与癌症有着密切的关系。综上所述，充分证明了实验所识别的调控子参与了大量的生物过程，尤其是与细胞活动及癌症有关的生物过程，因此具有很重要的生物意义。当前第1页1 2 3