CN110555530B - 一种基于分布式的大规模基因调控网络构建方法 - Google Patents

Abstract

本发明提供一种基于分布式的大规模基因调控网络构建方法，涉及医学信息学领域，本发明首先利用基因间的互信息值构建出基因调控的概要网络，根据结点的连结情况将网络分为若干个子网络，在每个子网络利用贝叶斯等方法在分布式环境下计算结点的调控方向，再将子网络进行合并，从而得到最终的基因调控网络。该发明可以提高基因调控网络的构建效率，通过分析基因调控网络让我们更系统地剖析细胞的功能，更深刻地洞见生命的本质。

Description

一种基于分布式的大规模基因调控网络构建方法

技术领域

本发明涉及医学信息学领域，具体涉及一种基于分布式的大规模基因调控网络构建方法。

背景技术

基因调控网络为人们认识复杂的生命过程提供了一个重要且有效的途径。基因调控网络是由若干个基因和这些基因间相互制约的关系构成的。基因间相互制约的关系简单来说就是，一个基因的表达受其它基因的影响，而该基因又会影响其他基因的表达。对基因调控网络的研究将有助于理解生命组织内部基因及其产物的生成过程和调控关系，可以实现对基因功能的整体认识和把握，从而对寻找和识别人类致病因子的研究有着重要的作用。基因调控网络是由一组基因、蛋白质、小分子以及它们之间的相互调控作用所构成的一种生化网络.它是一类基本且重要的生物网络。作为系统生物学的核心领域，构建和分析基因调控网络将有利于我们更系统地剖析细胞的功能，更深刻地洞见生命的本质。

随着基因芯片技术的不断发展，产生了大量的基因表达数据，这为构建基因调控网络奠定了基础。目前，对于基因调控网络的构建主要有两个方向，即通过生物学实验对网络进行绘制和通过重建算法对网络进行重建。而这两种方法都面临一些挑战，生物学实验所建立的基因调控网络往往不够完整准确，利用重建算法构建网络由于其复杂度较高，也存在着重建规模有限的问题。由于基因调控网络十分复杂，而所需构建的网络规模也较为庞大，利用重建算法对基因调控网络进行重建时会受限于单台计算机的计算能力，面对大量的基因表达数据以及复杂的基因调控网络就会显得力不从心。如何弥补单台计算机的计算缺陷、怎样充分利用好海量的基因表达数据构建大型的基因调控网络成为一种挑战。

随着大数据计算平台的不断发展，像Hadoop，Spark等分布式计算平台凭借着超强的计算能力在科学计算方面得到充分的运用，我们可以运用分布式计算的高效计算能力，解决构建大规模的基因调控网络效率低的问题。

发明内容

针对现有技术的不足，本发明的目的是为了解决在面对大规模的基因数量时，构建网络效率低的困境。首先利用基因间的互信息值构建出基因调控的概要网络，根据结点的连接情况将网络分为若干个子网络，在每个子网络利用贝叶斯等方法在分布式环境下计算结点的调控方向，再将子网络进行合并，从而得到最终的基因调控网络。该发明可以提高基因调控网络的构建效率，通过分析基因调控网络让我们更系统地剖析细胞的功能，更深刻地洞见生命的本质。

一种基于分布式的大规模基因调控网络构建方法，具体步骤如下：

1、一种基于分布式的大规模基因调控网络构建方法，其特征在于：包括以下步骤：

步骤1：利用基因之间的互信息值建立概要网络；

步骤2：将概要网络中的每个基因利用KNN算法进行网络分解，增加target gene的亲本集，对局部网络进行构建；

步骤3：将得到的局部网络放到分布式平台中，将每一个子网络在集群中不同的计算节点进行计算，通过分布式计算框架MapReduce对局部精确网络进行构建，即进行贝叶斯网络的构建；

步骤4：将由各个节点构建的局部精确网络进行合并，形成精确网络，利用条件互信息CMI去除假阳边，同时计算一阶CMI(i,j|k)，二阶CMI(i,j|k,l)，如果CMI(i,j|k)或CMI(i,j|k,l)小于设定阈值β，则将i和j相连的冗余边删除；

步骤5：将删除冗余边后的基因调控网络afterNet与删除前基因调控网络beforeNet对比，如果不一致，重复步骤2、步骤3、步骤4，直到对比一致，网络稳定不在变化，输出最终的基因调控网络。

所述步骤1的具体步骤为：

步骤1.1：使用标准离散化方法对基因数据进行预处理，将预处理后的基因表达数据处理成矩阵A形式，在基因的表达数据矩阵A中，其中每一列表示一个基因，每一行表示一个样本；

步骤1.2：将数据矩阵A中的列向量A_i与列向量A_j依次计算互信息值I_ij，其中A_i代表基因表达数据矩阵A的第i列，是第i个基因的表达数据；A_j代表基因表达数据矩阵A的第j列，是第j个基因的表达数据，I_ij代表基因i与基因j的互信息值；

步骤1.3：人为对阈值α进行设定，若互信息值I_ij大于等于α，则说明基因i与基因j存在调控关系，否则便不存在调控关系，由此生成概要网络G。

所述步骤2的具体步骤为：

步骤2.1：将概要网络进行无向化处理；

步骤2.2：依次遍历所述步骤1中的概要网络G，基因i相关联的基因j，基因j关联的基因，形成以target gene基因i为中心，与基因i相关联及间接相关的局部网络。

所述步骤3的具体步骤为：

步骤3.1：在分布式计算框架的第一阶段Mapper阶段中，以局部网络的targetgene为key值对局部网络进行划分，这样局部网络就会根据key值将局部网络分到各个节点上；

步骤3.2：在分布式计算框架的第二阶段Reducer中，利用globalMIT计算出targetgene的最优父结点，得到具有调控方向的局部精确网络。

本发明的有益效果：

本发明能够快速的构建出基因调控网络，是分析基因网络的基础，这样能够充分利用好基因数据，使构建大规模基因调控网络称为一种可能。基因网络的构建有助于我们更好的理解生命组织内部基因以及其产物的生成过程和调控关系，可以实现对基因功能的整体认识和把握，这对寻找和识别人类致病因子的研究有着重要的应用。

附图说明

图1为本发明实施例的基于分布式的大规模基因调控网络的构建方法流程图；

图2为本发明实施例的概要网络进行KNN网络分解；

图3为本发明实施例的概要网络分解网络流程图；

图4为本发明实施例ACC与阈值α之间的变化趋势图；

图5为本发明实施例ACC与阈值β之间的变化趋势图；

图6为本发明实施例的基因个数与时间趋势图；

图7为本发明实施例的分布式节点个数与时间趋势图；

具体实施方式

为了使本发明的目的、技术方案及优势更加清晰，下面结合附图和具体实施例对本发明做进一步详细说明。此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

一种基于分布式的大规模基因调控网络构建方法，其流程如图1所示，包括以下步骤：

步骤1：利用基因之间的互信息值来建立概要网络，概要网络只揭示基因间是否存在关系，而基因间的调控方向关系不得而知，建立概要网络的目的就是为了减少不相关或相关系数较小基因之间的计算，同时也为下面的基因调控网络的分解打下基础；

本实施例利用生成器生成大肠杆菌的模拟数据，本次使用的基因个数分别为150，200，250，300，Sample的个数分别为42，84，126，168，210，本次实验中使用基因个数为200，sample个数为126的数据，利用互信息值构建出概要网络；

步骤1.1：使用标准离散化方法对基因数据进行预处理，消除由于实验技术所导致的基因芯片的数据缺失或数据异常，使各个样本和平行实验的数据处于相同的水平，将预处理后的基因表达数据处理成矩阵A形式，将矩阵A的基因表达数据进行离散化，将离散化的基因表达数据进行数据拼接，最终形成sample为120的基因表达数据；

步骤1.2：将数据矩阵A中的列向量A_i与列向量A_j依次计算互信息值I_ij，其中A_i代表基因表达数据矩阵A的i列，也就是第i个基因的表达数据；A_j代表基因表达数据矩阵A的j列，也就是第j个基因的表达数据，I_ij为基因i与基因j的互信息值；

概要网络的构建采用互信息的方法来衡量基因之间的关系，互信息通常用于测量两个基因变量X和Y之间的可靠性，向量x和y代表基因X，Y的表达值向量。因此两个基因之间的相关性可以通过对基因表达数据的计算得出。两个离散变量X和Y之间的互信息可以通过熵的形式表现如下：

其中，p(x,y)为向量x和y的联合概率分布函数，p(x)和p(y)分别为向量x和y的边缘概率分布；H(X)和H(Y)分别为变量X和变量Y的熵，H(X,Y)为变量X和Y的联合熵。

由于基因表达数据普遍采用高斯分布假设，熵可以通过下面的高斯分布的概率密度函数来估计：

其中C是变量X的协方差矩阵，|C|是矩阵的行列式，N是样本数量，n是在C中的基因个数，一般来说，如果样本的数量几乎等于基因个数，则协方差矩阵通常用于估计基因表达谱中分布的协方差矩阵，可以认为是矩阵的很好的一种近似。然而，当样本数量小于基因个数时，正则化协方差矩阵用于估计基因表达分布的协方差矩阵。重复样本的数量将影响该方法的性能，并且增加的重复样本可以增强GRN推理算法的功效。因此，变量X的熵可以表示为：

通过以上公式分析，变量X与变量Y的MI(互信息值)可以由下面的公式得出：

其中，|C(X)|是变量X的协方差矩阵的行列式，|C(Y)|是变量Y的协方差矩阵的行列式，|C(X,Y)|是变量X和Y的协方差矩阵的行列式。将互信息计算公式转化成协方差的形式可以很大程度上提高互信息计算的效率。

应用到本实施例中具体公式为：

其中，MI(X_i,Y_j)为向量X_i和Y_j的互信息，|C(X_i)|是向量X_i的协方差矩阵的行列式，|C(Y_j)|是向量Y_j的协方差矩阵的行列式，|C(X_i,Y_j)|是向量X_i和Y_j的协方差矩阵的行列式。

得出高MI值表明变量(基因)X和Y之间可能存在密切关系，而低MI值意味着它们的独立性。如果变量(基因)X和Y彼此独立，则表明MI(X,Y)＝0。

步骤1.3：人工选择阈值α，若互信息值I_ij大于等于α，则说明基因i与基因j存在调控关系，否则就不存在调控关系。由此生成概要网络G。我们选用小规模10个基因进行仿真。我们得到ACC与阈值α之间的变化趋势图，如图4所示，根据ACC的变化趋势，我们将α范围设定为0.025≤α≤0.03，从而减少参数对构建网络的影响。

步骤2：将概要网络中的每个基因利用KNN算法进行网络分解，增加target gene的亲本集，对局部网络进行构建；

步骤2.1：将概要网络进行无向化处理。

步骤2.2：依次遍历所述步骤1中的概要网络G，基因i相关联的基因j，基因j关联的基因，形成以target gene基因i为中心，与基因i相关联及间接相关的局部网络；

k-nearest neighbor(KNN)算法KNN是一个基本而简单的分类算法，是一种基本分类与回归方法，在一个图中G(V,E)中，V代表一组结点，E代表结点之间的边。在图中k代表各个结点之间的最短路径距离。也就是说，结点Vi的k近邻(KNN)是由一组结点组成，这些结点到结点Vi的最短路径是k。其基本做法是：给定测试实例，基于某种距离度量找出训练集中与其最靠近的k个实例点，然后基于这k个最近邻的信息来进行预测，K值太小使得基因结点的搜索范围过小影响网络的构建精度。K值太大又会增加计算量影响构建网络的效率，实际上，k＝2的基因或结点的k近邻包含马尔可夫毯(在可信的贝叶斯网络中，一个结点的马尔可夫毯即为该结点的父结点、子结点以及子结点的父结点)结点的一部分，包括k＝1的所有k近邻，以及该结点k＝2的k近邻。所以我们将K设置为2；

如图2所示，与基因结点A直接相连的基因结点B、C、D，与基因结点B、C直接相连的分别为E和F，因为在本实施例中选择k＝2，而k代表各个结点之间的最短路径距离。所以A的亲本集为{B、C、D、E、F}。

步骤3：将得到的局部网络放到分布式平台中，将每一个子网络在集群中不同的计算节点进行计算，通过分布式计算框架MapReduce对局部精确网络进行构建，即进行贝叶斯网络的构建；

步骤3.1：在分布式计算框架的第一阶段Mapper阶段中，以局部网络的targetgene为key值对局部网络进行划分，这样局部网络就会根据key值将局部网络分到各个节点上；

如图3所示，在图中以1，2，3为target gene组成的局部概要网络，在Mapper中我们把target gene为key值，这样以target gene为中心组成的局部网络就会被分散到分布式的各个节点中。

步骤3.2：在分布式计算框架的第二阶段Reducer中，利用globalMIT计算出targetgene的最优父结点，得到具有调控方向的局部精确网络。

在此步骤中，我们在各个结点上构建局部精确网络，贝叶斯网络是一种有向的无环图，我们需要在分布式各个结点上构建出有方向的精确贝叶斯网络，利用的是globalMIT工具。该工具的具体内容如下：

globalMIT是一个从基因表达数据中学习全局最佳DBN结构的工具箱。globalMIT引入的基于信息理论的评分指标，称为互信息测试(MIT)。利用MIT，在多项式时间内有效地实现了学习全局最优DBN的任务。

互信息测验(MIT)是学习BN(贝叶斯网络)最近引入的一种评分标准。为了更好的理解MIT，我们让X＝{X₁,...,X_n}表示具有相应{r₁,...,r_n}离散状态的n个变量的集合，D表示我们的N个观测数据集，G表示DBN(动态贝叶斯网络)，并且

是G中X_i的父亲集合，具有相应的

离散状态，s_i＝|Pa_i|，MIT分数定义为：

其中，I(X_i,Pa_i)是X_i与其父亲集合的互信息值，

是

的值(显着性水平1-α的卡方分布)，l_iσ_i(j)被定义为：

其中σ_i＝{σ_i(1),...,σ_i(s_i)}是Pa_i的索引集{1...s_i}的任何排列，第一个变量的状态数最大，第二个是具有第二大状态数的变量，依此类推。

在各个节点上，每个子图以一个target gene基因为中心。然而，每个结点X_i的父集可以由G中的其他结点组成，即识别最佳子图的计算复杂度相当高，仍然无法计算大规模网络中每个结点的所有可能子图的贝叶斯分数。一般来说，基因X_i的邻居基因最有可能调节它。因此，我们限制每个结点的父结点的大小通过大致计算X_i每个结点的最大贝叶斯分数。我们可以识别每个中心结点的父结点集的所有可能子网结构的得分，以选择最优贝叶斯具有最大联合概率分布得分的子网。

步骤4：由各个结点构建的精确网络进行合并，形成精确网络，利用互信息往往高估了基因之间的调控强度，这是目前研究中普遍存在的问题，不考虑一个目标基因由另外两个或多个基因共同调控，从而导致更多的冗余边。在这一步中我们利用CMI(条件互信息)删除冗余边，通过计算一阶CMI(i,j|k),二阶CMI(i,j|k,l)，如果CMI(i,j|k)(或者CMI(i,j|k,l)小于一个预定的阈值β，就将i和j相连的冗余边删除。用此方法删除冗余的边。

所述冗余边为在正确的基因调控网络中本来不该存在的边，也就是错误的边。

类似地，在基因表达数据的高斯分布假设下，通过使用下面的等效公式可以容易地计算给定变量(基因)Z的两个变量(基因)X和Y的CMI。

显然，当给定Z时X和Y有条件地独立时，CMI(X,Y|Z)＝0。另外，这个等价表达式是给定一个或多个变量Z计算两个变量X和Y之间的CMI的有效方法，例如，如果条件变量Z＝(Z1，Z2)由两个变量Z1和Z2组成，我们可以获得二阶CMI。

人工选择阈值β，若互信息值CMI大于等于β，我们选用小规模10个基因进行仿真。我们得到ACC与阈值β之间的变化趋势图，如图5所示，根据ACC的变化趋势，我们将β设定0.024≤β≤0.03范围，从而减少参数对构建网络的影响。

步骤5：将删除冗余边后的基因调控网络afterNet与删除前基因调控网络beforeNet对比，如果不一致，重复步骤2、步骤3、步骤4。直到对比一致，网络稳定不在变化，输出最终的基因调控网络。

实验结果显示构建实例中基因个数为120，sample为126的时间为7389.984s，构建网络效率很高。我们将sample个数定为126，基因个数100、150、200、150、300。得到了构建网络时间随基因个数增加的趋势图，如图6所示，可以看出随着基因个数的增加，构建网络的时间会不断升高。另外，我们控制分布式平台的结点，将其结点设置为4,3,2,1，得到了趋势图如图7所示，显示出构建网络的速度会随着集群节点的个数的减少不断的降低。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解；其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；因而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (4)

1.一种基于分布式的大规模基因调控网络构建方法，其特征在于：包括以下步骤：

步骤1：利用基因之间的互信息值建立概要网络；

步骤2：将概要网络中的每个基因利用KNN算法进行网络分解，增加target gene的亲本集，对局部网络进行构建；

步骤3：将得到的局部网络放到分布式平台中，将每一个子网络在集群中不同的计算节点进行计算，通过分布式计算框架MapReduce对局部精确网络进行构建，即进行贝叶斯网络的构建；

步骤4：将由各个节点构建的局部精确网络进行合并，形成精确网络，利用条件互信息CMI去除假阳边，同时计算一阶CMI(i,j|k)，二阶CMI(i,j|k,l)，如果CMI(i,j|k)或CMI(i,j|k,l)小于设定阈值β，则将i和j相连的冗余边删除；

步骤5：将删除冗余边后的基因调控网络afterNet与删除前基因调控网络beforeNet对比，如果不一致，重复步骤2、步骤3、步骤4，直到对比一致，网络稳定不在变化，输出最终的基因调控网络。

2.根据权利要求1所述的一种基于分布式的大规模基因调控网络构建方法，其特征在于：

所述步骤1的具体步骤为：

步骤1.1：使用标准离散化方法对基因数据进行预处理，将预处理后的基因表达数据处理成矩阵A形式，在基因的表达数据矩阵A中，其中每一列表示一个基因，每一行表示一个样本；

步骤1.2：将数据矩阵A中的列向量A_i与列向量A_j依次计算互信息值I_ij，其中A_i代表基因表达数据矩阵A的第i列，是第i个基因的表达数据；A_j代表基因表达数据矩阵A的第j列，是第j个基因的表达数据，I_ij代表基因i与基因j的互信息值；

步骤1.3：人为对阈值α进行设定，若互信息值I_ij大于等于α，则说明基因i与基因j存在调控关系，否则便不存在调控关系，由此生成概要网络G。

3.根据权利要求1所述的一种基于分布式的大规模基因调控网络构建方法，其特征在于：所述步骤2的具体步骤为：

步骤2.1：将概要网络进行无向化处理；

步骤2.2：依次遍历所述步骤1中的概要网络G，基因i相关联的基因j，基因j关联的基因，形成以target gene基因i为中心，与基因i相关联及间接相关的局部网络。

4.根据权利要求1所述的一种基于分布式的大规模基因调控网络构建方法，其特征在于：所述步骤3的具体步骤为：

步骤3.1：在分布式计算框架的第一阶段Mapper阶段中，以局部网络的target gene为key值对局部网络进行划分，这样局部网络就会根据key值将局部网络分到各个节点上；

步骤3.2：在分布式计算框架的第二阶段Reducer中，利用globalMIT计算出targetgene的最优父结点，得到具有调控方向的局部精确网络。

Images