阅读说明：本技术 一种基于结构预测的基因调控网络构建方法 (Gene regulation and control network construction method based on structure prediction ) 是由 王之琼 郭上慧 曲路渲 信俊昌 钱唯 于 2019-09-17 设计创作，主要内容包括：本发明提供一种基于结构预测的基因调控网络构建方法。包括以下内容：首先计算系数矩阵,通过计算基因之间的Pearson系数、互信息及最大互信息来确定基因之间的相关性,作为筛选潜在父节点集的依据；然后进行结构预测,将获得的基因之间的系数矩阵作为判定基因潜在父节点集的依据,为每个基因选取潜在父节点集；最后进行基因调控网络的结构学习和参数学习。本发明通过基于Person系数、互信息及最大互信息相结合的方法预测基因的潜在父节点集,缩小结构学习的搜索范围,在一定程度上减少了基因调控网络的构建时间,提升计算性能,可以更加快速、准确地构建大规模基因调控网络,进一步了解生物的基因调控机制。(The invention provides a gene regulation and control network construction method based on structure prediction. The method comprises the following steps: firstly, calculating a coefficient matrix, and determining the correlation among genes by calculating Pearson coefficients, mutual information and maximum mutual information among the genes to serve as a basis for screening a potential father node set; then, carrying out structure prediction, taking the obtained coefficient matrix between the genes as a basis for judging a potential father node set of the genes, and selecting the potential father node set for each gene; and finally, performing structure learning and parameter learning of the gene regulation network. The method predicts the potential father node set of the gene by a method based on the combination of the Person coefficient, the mutual information and the maximum mutual information, reduces the search range of structure learning, reduces the construction time of the gene regulation and control network to a certain extent, improves the calculation performance, can construct the large-scale gene regulation and control network more quickly and accurately, and further knows the gene regulation and control mechanism of the organism.)

一种基于结构预测的基因调控网络构建方法

技术领域

本发明涉及医学信息学领域，尤其涉及一种基于结构预测的基因调控网络构建方法。

背景技术

基因表达调控即一个基因的表达被其它基因影响的过程，主要包括转录水平和翻译水平上的调控。由于基因表达的最终形式是蛋白质，需要转录形成mRNA作为模板最终产生蛋白质，因此转录水平上的基因调控是关键。一个基因的转录受到另外基因表达产物的影响，起到激励或抑制作用，其本身产生的蛋白质也可能影响其它基因。这种复杂的调控关系最终构成基因调控网络。了解生物的基因调控机制，可以在遗传视角上理解各种生物过程的发生，揭示生物体的不良过程，对发现生物系统的关键原理及细节至关重要。

目前的一些基因调控网络模型在一定程度上解决了基因调控网络的构建问题，但同时也存在一些问题。如布尔模型在定性的角度上研究基因调控网络，过于粗糙和简化；微分方程模型则通过微分方程来定量、精准地描述基因调控网络，但是同时会由于参数过多造成难以优化的问题，计算量巨大。贝叶斯网络模型通过概率模型来描述基因调控网络，用概率表示调控关系，但是随着网络复杂性的提高，该模型的计算量将大幅上升。

发明内容

针对现有技术的不足，本发明的目的是提出基于结构预测的基因调控网络构建方法。全局网络的构建可以分解为针对子网络的构建问题，在每个子网络的构建过程中，由于需要查找所有剩余基因可能组成的集合来得到最优的子网络结构，计算时间会随着基因数量的增长呈指数上升趋势，因此采用基于Pearson系数、互信息及最大互信息相结合的方法来预测每个基因的可能父节点集，在这个父节点集内搜索最优结构，极大程度上减小了搜索空间，提升计算性能。

为解决上述技术问题，本发明提出了一种基于结构预测的基因调控网络构建方法，包括以下步骤：

步骤1：计算系数矩阵，通过计算基因之间的Pearson系数、互信息及最大互信息来确定基因之间的相关性，作为筛选潜在父节点集的依据；

步骤2：结构预测，利用步骤1中获得的基因之间的系数矩阵，作为判定基因潜在父节点集的依据，为每个基因选取潜在父节点集；

步骤3：结构学习，采用基于评分搜索的结构学习方法，评分函数选择BDe分数，并设置每个基因的调控集合的个数为K，对于基因x_i，以潜在父节点集P_i为搜索空间，遍历所述搜索空间内的基因可能组成的集合，对每一个所述集合计算BDe分数score，根据分数高低评判所述集合作为基因x_i的父节点集的优劣；

步骤4：局部网络合并，每个基因对应一个局部网络G_i，将[G₁，……，G_n]合并成全局网络G[[x₁,G₁],[x₂,G₂],…,[x_n,G_n]]，设基因y∈G_i，则基因x_i和基因y间的存在调控关系表示为y→x_i；

步骤5：参数学习，对于所述全局网络G中的每一条调控关系进行参数学习，所述参数包括调控作用和调控概率，所述调控作用表示为激励或抑制，所述调控概率表示调控基因与靶基因的后验概率。

所述的步骤1计算系数矩阵，通过计算基因之间的Pearson系数、互信息及最大互信息来确定基因之间的相关性，作为筛选潜在父节点集的依据，具体步骤如下：

步骤1.1：将基因表达数据表示为一个矩阵X(m*n)，m表示基因表达数据的样本数，n为基因总数，所述矩阵X(m*n)的每一行为一个基因的表达数据向量，具体表述为第i个基因的表达数据向量表示为X_i，定义3个系数矩阵M_Pearson(n*n)用于存放基因之间的Pearson系数的计算结果，定义M_mi(n*n)用于存放互信息的计算结果，定义M_MIC(n*n)用于存放最大互信息的计算结果；

步骤1.2：由于Pearson系数、互信息和最大互信息三种方法都具有对称性的属性，因此在计算所述3个系数矩阵时只需计算上三角矩阵，利用公式(1)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的皮尔逊系数，将计算得到的皮尔逊系数存放到M_Pearson(n*n)中，具体为M_Pearson[i，i+1]，M_Pearson[i，i+2]，……，M_Pearson[i，n]，同时将计算得到的皮尔逊系数复制到对应的下三角，具体为M_Pearson[i+1，i]，M_Pearson[i+2，i]，……，M_Pearson[n，i]，最终得到M_Pearson的整个矩阵表达式；

式中，r(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的皮尔逊系数，基因y∈{x_i+1,x_i+2,…,x_n}，X_ip表示基因x_i的表达数据向量X_i的第p个样本表达值，Y_jp表示基因y的表达数据向量Y的第p个样本表达值，

表示基因x_i的表达数据向量X_i的样本平均值，

表示基因y的表达数据向量Y的样本平均值，m表示基因表达数据的样本数；

步骤1.3：利用公式(2)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的互信息的值，将计算得到的互信息值存放到M_mi(n*n)中，具体为M_mi[i，i+1]，M_mi[i，i+2]，……，M_mi[i，n]，同时将计算得到的互信息值复制到对应的下三角，具体为M_mi[i+1，i]，M_mi[i+2，i]，……，M_mi[n，i]，最终得到M_mi的整个矩阵表达式，

式中，I(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的互信息值，基因y∈{x_i+1,x_i+2,…,x_n}，|C(X_i)|表示基因x_i的表达数据向量X_i的协方差矩阵的行列式的值，|C(Y)|表示基因y的表达数据向量Y的协方差矩阵的行列式的值，|C(X_i,Y)|表示向量X_i和向量Y的协方差矩阵的行列式；

步骤1.4：利用公式(3)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的最大互信息的值，将计算得到的最大互信息值存放到M_MIC(n*n)中，即M_MIC[i，i+1]，M_MIC[i，i+2]，……，M_MIC[i，n]，同时将计算得到的最大互信息值复制到对应的下三角，具体为M_MIC[i+1，i]，M_MIC[i+2，i]，……，M_MIC[n，i]，最终得到M_MIC的整个矩阵表达式，

式中，mic(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的最大互信息值，基因y∈{x_i+1,x_i+2,…,x_n}，I(X_i,Y)表示所述基因x_i的表达数据向量X_i和基因y的表达数据向量Y的互信息值，a表示坐标轴X轴方向上的网格化的大小，b表示坐标轴Y轴方向上的网格化的大小，参数B取值为B＝m^0.6，其中m表示基因表达数据的样本数。

所述的步骤2结构预测，利用步骤1中获得的基因之间的系数矩阵，作为判定基因潜在父节点集的依据，为每个基因选取潜在父节点集，具体步骤如下：

步骤2.1：定义潜在父节点集比例为λ，则潜在父节点集选择的个数num＝λ*n，n为基因总数，λ的值通过进行实验确定；

步骤2.2：读取M_Pearson中基因x_i与除所述基因x_i之外的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个皮尔逊系数值对应的基因放入P_pearson中；

步骤2.3：读取M_mi中基因x_i与除所述基因x_i之外的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个互信息值对应的基因放入P_mi中；

步骤2.4：读取M_MIC中基因x_i与除所述基因x_i的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个最大互信息值对应的基因放入P_MIC中；

步骤2.5：对P_pearson、P_mi和P_MIC求并集得到集合P_i，并将所述集合P_i作为基因x_i的潜在父节点集。

所述的步骤3结构学习，采用基于评分搜索的结构学习方法，评分函数选择BDe分数，并设置每个基因的调控集合的个数为K，对于靶基因X_i，以潜在父节点集P_i为搜索空间，遍历所述搜索空间内的基因可能组成的集合，对每一个所述集合计算BDe分数，根据分数高低评判所述集合作为基因x_i的父节点集的优劣，具体表述为：

步骤3.1：初始化一个小顶堆S，并定义所述小顶堆S的大小为K，利用公式(4)计算基因x_i的父节点集为空集时的BDe分数score₀，并将

存入所述小顶堆S中；

式中，D＝{D₁,D₂,…,D_n}表示关于基因{x₁,x₂,…,x_n}的样本集，G表示以x₁,x₂,…,x_n为节点的贝叶斯网络，Γ表示Gamma函数，q_i是第i个基因x_i的父节点的表达数据向量中的取值组合的数目，r_i是基因x_i的表达数据向量X_i中的取值数目，m_ijk表示基因x_i的父节点取j值，x_i取k值时的样本数目，α_ijk表示服从Dirichlet分布的函数的参数，

步骤3.2：以潜在父节点集P_i为搜索空间，并对所述搜索空间内的基因进行排列组合，遍历所述搜索空间，根据公式(5)计算得到H个集合，对所述H个集合中的每一个集合利用公式(4)计算BDe分数score_x，并依次与小顶堆堆顶的元素中的BDE分数score进行比较，若score_x>score，则将score_x以及score_x对应的集合存入小顶堆S中，若小顶堆S已存满，则先删除堆顶元素，再将score_x及score_x对应的集合存入小顶堆S中；若score_x≤score，则不进行存入操作，x属于自然数，且满足0<x<H；

式中，g表示所述搜索空间内的基因总数，f表示当前遍历的集合内的基因数量；

步骤3.3：将所述小顶堆S的H个集合存入基因x_i的局部网络G_i中。

所述的步骤5参数学习，对于所述全局网路G中的每一条调控关系进行参数学习，所述参数包括调控作用和调控概率，所述调控作用表示为激励或抑制，所述调控概率表示调控基因与靶基因的后验概率，具体步骤如下：

步骤5.1：对于调控关系y→x_i，基因y作为调控基因，基因x_i作为靶基因，利用公式(1)计算基因x_i和基因y的皮尔逊系数，若所述皮尔逊系数大于0，则定义基因x_i和基因y调控关系为激励关系，用“+”表示；若所述皮尔逊系数小于0，则定义基因x_i和基因y调控关系为抑制关系，用“-”表示；

步骤5.2：对于调控关系y→x_i，利用公式(6)计算基因x_i和基因y的后验概率，

式中，P(Y|X_i)表示通过基因x_i的表达数据向量X_i和基因y的表达数据向量Y计算的后验概率，P(X_i)表示基因x_i的表达数据向量X_i相对于所有调控基因的表达数据向量的联合条件概率，P(Y)表示基因y的表达数据向量Y的概率分布，P(X_i|Y)表示基因y作为基因x_i的父节点时计算的条件概率；

所述联合条件概率的计算公式如下：

式中，pa(x_i)表示基因x_i的父节点集，P(X_i|Y)表示基因y作为基因x_i的父节点时计算的条件概率。

本发明的有益效果是：

本发明是一种基于结构预测的基因调控网络构建方法，通过基于Person系数、互信息及最大互信息相结合的方法预测基因的潜在父节点集，缩小结构学习的搜索范围，在一定程度上减少了基因调控网络的构建时间，提升计算性能，可以更加快速、准确地构建大规模基因调控网络，进一步了解生物的基因调控机制。

附图说明

图1为本发明实施例中的基于结构预测的基因调控网络的构建方法流程图。

图2为本发明实施例中的结构预测阶段流程图。

图3为本发明实施例中的父节点比例实验结果。

图4为本发明实施例中的基因调控网络结构学习过程流程图。

图5为本发明实施例中的基因调控网络参数学习过程流程图。

图6为本发明实施例中的基于结构预测的基因调控网络构建方法时间对比图。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明：一种基于结构预测的基因调控网络的构建方法，其流程如图1所示，包括计算系数矩阵、结构预测过程、基因调控网络的结构学习过程和基因调控网络的参数学习过程，具体包括如下步骤：

步骤1：计算系数矩阵，通过计算基因之间的Pearson系数、互信息及最大互信息来确定基因之间的相关性，作为筛选潜在父节点集的依据，将相关性较强的基因设为当前基因的可能父节点集，本实施实例选取大肠杆菌基因来构建基因调控网络，具体步骤如下：

步骤1.1：将基因表达数据表示为一个矩阵X(m*n)，m表示基因表达数据的样本数，n为基因总数，所述矩阵X(m*n)的每一行为一个基因的表达数据向量，具体表述为第i个基因的表达数据向量表示为X_i，定义3个系数矩阵M_Pearson(n*n)用于存放基因之间的Pearson系数的计算结果，定义M_mi(n*n)用于存放互信息的计算结果，定义M_MIC(n*n)用于存放最大互信息的计算结果；

步骤1.2：由于Pearson系数、互信息和最大互信息三种方法都具有对称性的属性，因此在计算所述3个系数矩阵时只需计算上三角矩阵，利用公式(1)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的皮尔逊系数，将计算得到的皮尔逊系数存放到M_Pearson(n*n)中，具体为M_Pearson[i，i+1]，M_Pearson[i，i+2]，……，M_Pearson[i，n]，同时将计算得到的皮尔逊系数复制到对应的下三角，具体为M_Pearson[i+1，i]，M_Pearson[i+2，i]，……，M_Pearson[n，i]，最终得到M_Pearson的整个矩阵表达式，Pearson系数广泛用于度量两个变量之间的相关程度，其值介于-1与1之间，

式中，r(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的皮尔逊系数，基因y∈{x_i+1,x_i+2,…,x_n}，X_ip表示基因x_i的表达数据向量X_i的第p个样本表达值，Y_jp表示基因y的表达数据向量Y的第p个样本表达值，表示基因x_i的表达数据向量X_i的样本平均值，

表示基因y的表达数据向量Y的样本平均值，m表示基因表达数据的样本数；

步骤1.3：利用公式(2)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的互信息的值，将计算得到的互信息值存放到M_mi(n*n)中，具体为M_mi[i，i+1]，M_mi[i，i+2]，……，M_mi[i，n]，同时将计算得到的互信息值复制到对应的下三角，具体为M_mi[i+1，i]，M_mi[i+2，i]，……，M_mi[n，i]，最终得到M_mi的整个矩阵表达式，互信息通常用于测量两个变量X和Y之间的可靠性，因此两个基因之间的相关性可以通过对基因表达数据的计算得出；

式中，I(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的互信息值，基因y∈{x_i+1,x_i+2,…,x_n}，|C(X_i)|表示基因x_i的表达数据向量X_i的协方差矩阵的行列式的值，|C(Y)|表示基因y的表达数据向量Y的协方差矩阵的行列式的值，|C(X_i,Y)|表示向量X_i和向量Y的协方差矩阵的行列式；

步骤1.4：利用公式(3)计算基因x_i的表达数据向量X_i与基因x_i+1,x_i+2,…,x_n的表达数据向量X_i+1,X_i+2,…X_n的最大互信息的值，将计算得到的最大互信息值存放到M_MIC(n*n)中，即M_MIC[i，i+1]，M_MIC[i，i+2]，……，M_MIC[i，n]，同时将计算得到的最大互信息值复制到对应的下三角，具体为M_MIC[i+1，i]，M_MIC[i+2，i]，……，M_MIC[n，i]，最终得到M_MIC的整个矩阵表达式，最大互信息，属于Maximal Information-based Nonparametric Exploration(MINE)最大的基于信息的非参数性探索，用于衡量两个变量X和Y的线性或非线性的强度，

式中，mic(X_i,Y)表示基因x_i的表达数据向量X_i和基因y的表达数据向量Y的最大互信息值，基因y∈{x_i+1,x_i+2,…,x_n}，I(X_i,Y)表示所述基因x_i的表达数据向量X_i和基因y的表达数据向量Y的互信息值，a、b分别表示坐标轴X轴和Y轴方向上的网格化的大小(本质上就是网格分布)，参数B取值为B＝m^0.6，其中m表示基因表达数据的样本数。

步骤2：结构预测，利用步骤1中获得的基因之间的系数矩阵，作为判定基因潜在父节点集的依据，为每个基因选取潜在父节点集，具体步骤如下：

步骤2.1：定义潜在父节点集比例为λ，则潜在父节点集选择的个数num＝λ*n，n为基因总数，λ的值通过进行实验确定；

如图3所示为本发明实施例的父节点比例实验结果，实验采用大肠杆菌的表达数据，选取的基因数分别为10、20、30、40和50，样本数为1565，设置父节点集比例为分别为0.1、0.2、0.3、0.4，将本方法与基于贝叶斯网络的方法作对比，根据最终的基因调控网络结果进行效果比较，主要的指标包括TPR(真阳率)、FPR(假阳率)、PPV(精确率)、Sensitivity(敏感度)、Specificity(特异度)、RECALL(召回率)、ACC(准确率)和F-score，由结果可以得出，当父节点比例达到0.2时，本方法已与对比方法达到了近似的效果，比例达到0.3和0.4时，效果呈现一种收敛状态，但时间效果降低，因此父节点比例可设为0.2。

步骤2.2：读取M_Pearson中基因x_i与除所述基因x_i之外的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个皮尔逊系数值对应的基因放入P_pearson中；

步骤2.3：读取M_mi中基因x_i与除所述基因x_i之外的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个互信息值对应的基因放入P_mi中；

步骤2.4：读取M_MIC中基因x_i与除所述基因x_i的其他基因之间的皮尔逊系数结果并排序，选择绝对值最高的前num个最大互信息值对应的基因放入P_MIC中；

步骤2.5：对P_pearson、P_mi和P_MIC求并集得到集合P_i，并将所述集合P_i作为基因x_i的潜在父节点集。

图2为本发明具体实施方式中结构预测过程流程图，预测基因x_i的潜在父节点集，其中V₁，V₂，V₃，中的1，……，n分别为与基因x_i进行计算的基因的下标值，value为对应的从系数矩阵中取出的系数值。

步骤3：结构学习，采用基于评分搜索的结构学习方法，评分函数选择BDe(Bayesian Dirichlet)分数，BDe评分函数(Bayesian Dirichlet)，利用先验知识和数据寻求概率最大的网络结构，并设置每个基因的调控集合的个数为K，对于基因x_i，以潜在父节点集P_i为搜索空间，遍历所述搜索空间内的基因可能组成的集合(结构)，对每一个所述集合计算BDe分数score，根据分数高低评判所述集合作为基因x_i的父节点集的优劣，基因调控网络结构学习过程流程图如图4所示，具体步骤如下：

步骤3.1：初始化一个小顶堆S，并定义所述小顶堆S的大小为K，利用公式(4)计算基因x_i的父节点集为空集时的BDe分数score₀，并将

存入所述小顶堆S中；

式中，D＝{D₁,D₂,…,D_n}表示关于基因{x₁,x₂,…,x_n}的样本集，G表示以x₁,x₂,…,x_n为节点的贝叶斯网络，Γ表示Gamma函数，q_i是第i个基因x_i的父节点的表达数据向量中的取值组合的数目，r_i是基因x_i的表达数据向量X_i中的取值数目，m_ijk表示基因x_i的父节点取j值，x_i取k值时的样本数目，α_ijk表示服从Dirichlet分布的函数的参数，

步骤3.2：以潜在父节点集P_i为搜索空间，并对所述搜索空间内的基因进行排列组合，遍历所述搜索空间，根据公式(5)计算得到H个集合(结构)，对所述H个集合中的每一个集合利用公式(4)计算BDe分数score_x，并依次与小顶堆堆顶的元素中的BDE分数score进行比较，若score_x>score，则将score_x以及score_x对应的集合存入小顶堆S中，若小顶堆S已存满，则先删除堆顶元素，再将score_x及score_x对应的集合存入小顶堆S中；若score_x≤score，则不进行存入操作，x属于自然数，且满足0<x<H；

式中，g表示所述搜索空间内的基因总数，f表示当前遍历的集合内的基因数量；

步骤3.3：将所述小顶堆S的H个集合存入基因x_i的局部网络G_i中。

步骤4：局部网络合并，每个基因对应一个局部网络G_i，将[G₁，……，G_n]合并成全局网络G[[x₁,G₁],[x₂,G₂],…,[x_n,G_n]]，设基因y∈G_i，则基因x_i和基因y间的存在调控关系表示为y→x_i。

步骤5：参数学习，对于所述全局网络G中的每一条调控关系进行参数学习，所述参数包括调控作用和调控概率，所述调控作用表示为激励或抑制，所述调控概率表示调控基因与靶基因的后验概率，基因调控网络参数学习过程流程图如图5所示，具体步骤如下：

步骤5.1：对于调控关系y→x_i，基因y作为调控基因，基因x_i作为靶基因，利用公式(1)计算基因x_i和基因y的皮尔逊系数，若所述皮尔逊系数大于0，则定义基因x_i和基因y调控关系为激励关系，用“+”表示；若所述皮尔逊系数小于0，则定义基因x_i和基因y调控关系为抑制关系，用“-”表示；

步骤5.2：对于调控关系y→x_i，利用公式(6)计算基因x_i和基因y的后验概率，

式中，P(Y|X_i)表示通过基因x_i的表达数据向量X_i和基因y的表达数据向量Y计算的后验概率，P(X_i)表示基因x_i的表达数据向量X_i相对于所有调控基因的表达数据向量的联合条件概率，P(Y)表示基因y的表达数据向量Y的概率分布，P(X_i|Y)表示基因y作为基因x_i的父节点时计算的条件概率；

所述联合条件概率的计算公式如下：

式中，pa(x_i)表示基因x_i的父节点集，P(X_i|Y)表示基因y作为基因x_i的父节点时计算的条件概率。

图6为本发明实施例的基于结构预测的基因调控网络构建方法与基于贝叶斯模型的基因调控网络构建方法所用时间对比图，分别采用基因数量为10、20、30、40和50的基因表达数据，记录父节点比例分别为0.1、0.2、0.3、0.4时的时间，结果表明本方法在计算时间方面优于传统方法。