CN109243523A - 一种基于乳腺癌疾病的调控网络构建及分析方法 - Google Patents

Abstract

本发明提出一种基于乳腺癌疾病的基因调控网络构建及分析方法，流程包括：乳腺癌相关基因初步筛选；基因调控网络的构建；根据基因调控网络结果，进行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性；取网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据作为为筛选出的相关基因；在人类基因中有效地筛选出了乳腺癌相关基因，从而建立与乳腺癌疾病相关的基因调控网络，并通过节点中心性分析得到重要基因，可以促进从基因学的角度对乳腺癌疾病的研究，为找到干预乳腺癌疾病发生的有效途径奠定基础。

Description

一种基于乳腺癌疾病的调控网络构建及分析方法

技术领域

本发明属于医学信息学领域，涉及一种基于乳腺癌疾病的基因调控网络构建及分析方法。

背景技术

遗传表达确定了细胞和组织的表型和发育状态。人类常见的疾病，如恶性肿瘤和神经退行性疾病都可以追溯到异常基因表达。通常，基因不是单独表达的，一个基因表达可以影响其他基因的表达，同时也受其他基因表达的影响。基因之间的相互作用和相互关联的调节关系形成了基因调控网络。基因调控网络的建模和分析为乳腺癌疾病提供了遗传视角的治疗和诊断。因此，构建与疾病相关的基因调控网络对探索人体机制和疾病治疗具有重要意义。

现在的网络建模方法有效地解决了基因调控网络构建的问题，但也存在一些缺点。例如，贝叶斯网络模型可以准确地建立监管网络，但仅适用于构建小规模网络。互信息可以支持大规模网络建模，但不能描述基因之间的调节关系。多模型组合可以改善上述两个问题，但仍不能有效满足基因组规模调节网络建模的需求。

发明内容

针对现有技术的不足，本发明的目的是提出基于乳腺癌疾病的基因调控网络构建及分析方法，利用基因对之间的互信息值进行逐轮筛选，在全基因组数据中筛选出与乳腺癌疾病相关的200个基因进行基因调控网络构建并进行网络分析，为乳腺癌疾病的研究人员提供有价值的参考。

本发明提出一种基于乳腺癌疾病的基因调控网络构建及分析方法，包括以下步骤：

步骤1：乳腺癌相关基因初步筛选：将乳腺癌基因分为基准基因、相关基因与筛选出的相关基因，基准基因指现有技术中已经认定与乳腺癌疾病有关系的基因，相关基因指除现有技术中以外的所有基因，筛选出的相关基因指用本申请的方法筛选出与乳腺癌疾病关系紧密的基因，包括步骤1.1～步骤1.4：

步骤1.1：将基因表达数据分为两个矩阵，X矩阵为乳腺癌基准基因的表达数据，Y矩阵是待筛选的其余乳腺癌相关基因的表达数据，在X矩阵与Y矩阵中，每一行为一个基因，每一列为一个样本；

步骤1.2：矩阵Y中的向量Y_j依次与X中的向量X_i计算每个基因之间的互信息值，Y_j和X_i代表第j个和第i个基因的数据；

步骤1.3：人工选择阈值I_sn，若互信息值大于等于I_sn，则向量Y_j对应的基因添加到X矩阵中，并在矩阵Y中删除，阈值I_sn为动态选择的，根据每一轮的互信息结果，对每轮的互信息的值大小进行排序，I_sn即为当前轮中第sn个阈值，筛选出来互信息值最大的前sn个基因；

步骤1.4：重复步骤1.2至1.3，直到矩阵X中的基因数量达到预设定的值M个；

步骤2：基因调控网络的构建：利用步骤1中获得基因的表达数据矩阵X，通过BNFinder2工具进行调控网络构建，设置调控节点集合的数量，评分函数为BDE评分函数，得到一个概率网络，网络节点代表基因，有向边代表基因之间的调控关系；具体包括步骤2.1～步骤2.3：

步骤2.1：数据处理：根据X矩阵，确定基因个数及样本数量，并为每个基因设定潜在调控基因；

步骤2.2：设定基因调控网络结构的结构学习及参数学习，采用基于评分搜索的结构学习方法，遍历所有可能的结构，然后用BDE评分函数衡量各个结构，进而找出最好的结构；具体包括步骤2.2.1～步骤2.2.4：

步骤2.2.1：对于某个基因，首先获取该基因的潜在调控基因及其表达数据；

步骤2.2.2：为基因的每个潜在调控节点赋予权重；

步骤2.2.3：获取调控基因为空集的基因的BDE分数；

步骤2.2.4：遍历基因的潜在调控基因的集合，每次遍历结果即为该基因的一种调控结构，计算该结构的BDE分数，存入到一个堆栈中，遍历结束后，堆栈中留下的为评分最高的对应的调控基因的集合，堆栈堆的大小即为设定的调控节点集合的数量；

步骤2.3：根据网络的学习结果整合调控关系，得到基因调控网络结果。；

步骤3：根据基因调控网络结果，进行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性，将BNFinder2得到的基因调控网络，通过Python的扩展包networkx计算各个基因的度中心性，接近中心性，中介中心性和特征向量中心性；

步骤4：取网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据作为为筛选出的相关基因。

所述网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据，N≤10，且矩阵X中的基因数量M>N。

所述步骤1.3中人工选择阈值I_sn，其中sn等于需要筛选出来的基因数量除以需要的轮数。

所述基准基因包括基因名称如下：AKT1、AKT2、APC、ARID1A、ARID1B、ARID2、ASXL1、BAP1、BRCA1、BRCA2、CASP8、CDH1、CDKN1B、CDKN2A、CCND1、MDM2、99Mb、ERBB2、ZNF217、FGFR1/ZNF703、MYC、GATA3、KRAS、MAP2K4、MAP3K1、MAP3K13、MIL2、MIL3、NCOR1、NF1、PIK3CA、PTEN、RB1、SETD2、SF3B1、SMAD4、SMARCD1、STK11、TBX3和TP53。

所述互信息的计算方法，如下公式：

其中，I(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的协方差矩阵的行列式。

所述节点度中心性(Degree Centrality)：在一个网络图G＝(V,E)中，节点v的度d_v指的是与v相连的E中边的数量，公式为：

其中，C_B(v)为节点度中心性，m为该节点的度，n为基因的数量；

所述接近中心性(Closeness Centrality)：考量每个结点到其它结点的最短路的平均长度。距离其它结点越近，那么它的中心度越高，公式为：

其中，C_cl(v)为接近中心性，n为基因的数量，d(v,u)为节点v和u之间的最短路径。

所述中介中心性(Betweenness Centrality)：以经过某个节点的最短路径数目来刻画节点重要性的指标。一个结点充当“中介”的次数越高，它的中介中心度就越大，公式为：

其中，C_B(v)为中介中心性，σ(s,t|v)表示经过节点v的s到t的最短路径条数，σ(s,t)表示s到t的最短路径条数。

特征向量中心性(Eigenvector Centrality)：取决于其相邻节点的数量，即该节点的度，也取决于其相邻节点的重要性。特征向量中心性基于其相邻节点的中心性来计算节点的中心性。节点i的特征向量中心性是由以下等式定义的向量y的第i个元素：

Ay＝λy

其中，A是网络的邻接矩阵，特征值为λ，y为解上述方程得出的特征向量。

有益技术效果：

本发明是一种基于乳腺癌疾病的基因调控网络构建及分析方法，在人类基因中有效地筛选出了乳腺癌相关基因，从而建立与乳腺癌疾病相关的基因调控网络，并通过节点中心性分析得到重要基因，可以促进从基因学的角度对乳腺癌疾病的研究，为找到干预乳腺癌疾病发生的有效途径奠定基础。

附图说明

图1为本发明实施例的基于乳腺癌疾病的基因调控网络的构建及分析方法流程图；

图2为本发明实施例的乳腺癌相关基因逐轮筛选策略单轮筛选过程流程图；

图3为本发明实施例的乳腺癌相关基因的基因调控网络构建过程流程图；

图4为本发明实施例的乳腺癌相关基因的基因调控网络节点分析过程流程图；

图5为本发明实施例的调控基因的网络；

图6为本发明实施例的互信息计算结果截图

图7为本发明实施例的乳腺癌相关基因的基因节点中心性分析计算结果截图。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明：一种基于乳腺癌疾病的基因调控网络的构建及分析方法，其流程如图1所示，包括乳腺癌相关基因初步筛选过程、基因调控网络建模过程和基因调控网络分析过程，具体包括：

步骤1：乳腺癌相关基因初步筛选：将乳腺癌基因分为基准基因、相关基因与筛选出的相关基因，基准基因指现有技术中已经认定与乳腺癌疾病有关系的基因，一共40个基因，相关基因指除现有技术以外的所有基因，本实施例取24236个基因，筛选出的相关基因指用本申请的方法筛选出与乳腺癌疾病关系紧密的基因，包括步骤1.1～步骤1.4：

步骤1.1：将基因表达数据分为两个矩阵，X矩阵为乳腺癌基准基因的表达数据，Y矩阵是待筛选的其余乳腺癌相关基因的表达数据，在X矩阵与Y矩阵中，每一行为一个基因，每一列为一个样本；

步骤1.2：矩阵Y中的向量Y_j依次与X中的向量X_i计算每个基因之间的互信息值，Y_j和X_i代表第j个和第i个基因的数据；

基因筛选过程中采用的互信息的方法来衡量基因之间的关系。互信息通常用于测量两个变量X和Y之间的可靠性，因此两个基因之间的相关性可以通过对基因表达数据的计算得出。两个离散变量X和Y之间的互信息可以通过熵的形式表现如下：

I(X,Y)＝H(X)+H(Y)-H(X,Y)

其中，H(X)和H(Y)分别为变量X和变量Y的熵，H(X,Y)为变量X和Y的联合熵，离散变量X的熵H(X)可以表示如下：

其中，p(x)是向量X的边缘概率分布函数。向量X和Y的联合熵可以表示为

其中，p(x,y)为X和Y的联合概率分布函数。熵通过高斯核概率密度估计可以表示如下：

其中，C是变量X的协方差矩阵，n为C中的变量(基因)数量，由此，向量X和Y的互信息可以表示为协方差的形式：

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

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

其中，I(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的协方差矩阵的行列式。

如图6所示，为本发明实施例的互信息计算结果截图，探针编号是原始的基因表达数据中基因名称对应的探针编号，一个基因可能对应多个探针，最后一列数据代表基准基因来自第几轮的筛选出来的基因，0则代表是来自最开始的40个基因中的。

步骤1.3：人工选择阈值I_sn，若互信息值大于等于I_sn，则向量Y_j对应的基因添加到X矩阵中，并在矩阵Y中删除，阈值I_sn为动态选择的，根据每一轮的互信息结果，对每轮的互信息的值进行排序，I_sn即为当前轮中第sn个阈值，筛选出来互信息值最大的前sn个基因，其中sn等于需要初步筛选出来的基因数量除以需要的轮数，需要初步筛选出来的基因数为200，本实施例轮数设置为10；

步骤1.4：重复步骤1.2至1.3，直到矩阵X中的基因数量达到预设定的值M个，M＝240；

图2为本发明具体实施方式中乳腺癌相关基因逐轮筛选策略单轮筛选过程流程图，每轮过程一样，直到筛选到需要初步筛选出来的基因数量，本实施例设置为200个；

所述基准基因包括基因名称如下：AKT1、AKT2、APC、ARID1A、ARID1B、ARID2、ASXL1、BAP1、BRCA1、BRCA2、CASP8、CDH1、CDKN1B、CDKN2A、CCND1、MDM2、99Mb、ERBB2、ZNF217、FGFR1/ZNF703、MYC、GATA3、KRAS、MAP2K4、MAP3K1、MAP3K13、MIL2、MIL3、NCOR1、NF1、PIK3CA、PTEN、RB1、SETD2、SF3B1、SMAD4、SMARCD1、STK11、TBX3和TP53。

步骤2：基因调控网络的构建：利用步骤1中获得基因的表达数据矩阵X，包括240个基因数据，通过BNFinder2工具进行调控网络构建，设置调控节点集合的数量，评分函数为BDE评分函数，得到一个概率网络，网络节点代表基因，有向边代表基因之间的调控关系；具体包括步骤2.1～步骤2.3：

步骤2.1：数据处理：根据X矩阵，确定基因个数及样本数量，并为每个基因设定潜在调控基因，如图5所示，为本发明实施例的调控基因的网络，箭头表示调控方向，虚线表示负调控，实线是正调控；

步骤2.2：设定基因调控网络结构的结构学习及参数学习，采用基于评分搜索的结构学习方法，遍历所有可能的结构，然后用BDE评分函数衡量各个结构，进而找出最好的结构；具体包括步骤2.2.1～步骤2.2.4：

步骤2.2.1：对于某个基因，首先获取该基因的潜在调控基因及其表达数据；

步骤2.2.2：为基因的每个潜在调控节点赋予权重；

步骤2.2.3：获取调控基因为空集的基因的BDE分数；

步骤2.2.4：遍历基因的潜在调控基因的集合，每次遍历结果即为该基因的一种调控结构，计算该结构的BDE分数，存入到一个堆栈中，遍历结束后，堆栈中留下的为评分最高的对应的调控基因的集合，堆栈堆的大小即为设定的调控节点集合的数量；

步骤2.3：根据网络的学习结果整合调控关系，得到基因调控网络结果；

图3为本发明具体实施方式中乳腺癌相关基因的基因调控网络构建过程流程图。利用上一步骤获得的240个基因的表达数据通过BNFinder2工具进行调控网络构建，调控节点集合的数量设置为5，评分函数为BDE评分函数，可以得到一个概率网络，网络节点代表基因，有向边表示基因之间的调控关系。在此过程中，依次基于每个基因学习它的调控网络，待全部基因的子网络都构建完后合并成一个全局网络，即最终的基因调控网络。

本实施方式中，基因调控网络的构建采用的BNFinder2工具，该工具的具体内容如下：

BNFinder2是一个构建贝叶斯网络准确并有效的工具，静态贝叶斯和动态贝叶斯网络(DBN)都适用。一个贝叶斯网络定义包括一个有向无环图(DAG)和一个条件概率表集合。DAG中每一个节点表示一个随机变量{X₁，X₂，……，X_n}，可以是可直接观测变量或隐藏变量，而有向边表示随机变量间的条件依赖；条件概率表中的每一个元素对应DAG中唯一的节点，存储此节点对于其所有直接前驱节点的联合条件概率节点X的联合概率可以表示成：

其中，pa(i)表示节点i之“因”，或称pa(i)是i的parents(父母)。

BNFinder2采用基于评分搜索的结构学习方法，遍历所有可能的结构并采用合适的评分函数衡量结构好坏，进而找出最好的结构。采用的评分函数为BDe评分函数(BayesianDirichlet)，利用先验知识和数据寻求概率最大的网络结构，函数如下：

其中，Γ为Gamma函数，m_i'_jk是数据D中节点V_i状态为k且父节点状态组合为j的样本数，m′_ij＝∑_km′_ijk，α_ijk为服从Dirichlet分布的函数的参数，

步骤3：根据基因调控网络结果，进行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性，将BNFinder2得到的基因调控网络，通过Python的扩展包networkx计算各个基因的度中心性，接近中心性，中介中心性和特征向量中心性；

所述节点度中心性(Degree Centrality)：在一个网络图G＝(V,E)中，节点v的度d_v指的是与v相连的E中边的数量，公式为：

其中，C_B(v)为节点度中心性，m为该节点的度，n为基因的数量；

所述接近中心性(Closeness Centrality)：考量每个结点到其它结点的最短路的平均长度。距离其它结点越近，那么它的中心度越高，公式为：

其中，C_cl(v)为接近中心性，n为基因的数量，d(v,u)为节点v和u之间的最短路径。

所述中介中心性(Betweenness Centrality)：以经过某个节点的最短路径数目来刻画节点重要性的指标。一个结点充当“中介”的次数越高，它的中介中心度就越大，公式为：

其中，C_B(v)为中介中心性，σ(s,t|v)表示经过节点v的s到t的最短路径条数，σ(s,t)表示s到t的最短路径条数。

特征向量中心性(Eigenvector Centrality)：取决于其相邻节点的数量(即该节点的度)，也取决于其相邻节点的重要性。特征向量中心性基于其相邻节点的中心性来计算节点的中心性。节点i的特征向量中心性是由以下等式定义的向量x的第i个元素：

Ay＝λy

其中，A是网络的邻接矩阵，特征值为λ，y为解上述方程得出的特征向量。

如图7所示截图为根据因调控网络结果，行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性计算结果。

步骤4：取网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据作为筛选出的相关基因,，N＝10。

图4为本发明具体实施方式中乳腺癌相关基因的基因调控网络节点分析过程流程图。对生成的基因调控网络进行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性。将BNFinder得到的基因调控网络通过Python的扩展包networkx计算各个基因的度中心性，接近中心性，中介中心性和特征向量中心性。依据节点中心性的分析结果，挑选出四种度量指标中值最高的10个节点，作为筛选出的相关基因的结果，如表1所示。

表1筛选出的相关基因的结果

Claims (6)

1.一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，包括以下步骤：

步骤1：乳腺癌相关基因初步筛选：将乳腺癌基因分为基准基因、相关基因与筛选出的相关基因，基准基因指现有技术中已经认定与乳腺癌疾病有关系的基因，相关基因指除基准基因以外所有基因，筛选出的相关基因指用本申请的方法筛选出与乳腺癌疾病关系紧密的基因，包括步骤1.1～步骤1.4：

步骤1.1：将基因表达数据分为两个矩阵，X矩阵为乳腺癌基准基因的表达数据，Y矩阵是待筛选的其余乳腺癌相关基因的表达数据，在X矩阵与Y矩阵中，每一行为一个基因，每一列为一个样本；

步骤1.2：矩阵Y中的向量Y_j依次与X中的向量X_i计算每个基因之间的互信息值，Y_j和X_i代表第j个和第i个基因的数据；

步骤1.3：人工选择阈值I_sn，若互信息值大于等于I_sn，则向量Y_j对应的基因添加到X矩阵中，并在矩阵Y中删除，阈值I_sn为动态选择的数值，根据每一轮的互信息结果，对每轮的互信息的值大小进行排序，I_sn即为当前轮中第sn个阈值，筛选出来互信息值最大的前sn个基因；

步骤1.4：重复步骤1.2至1.3，直到矩阵X中的基因数量达到预设定的值M个；

步骤2：基因调控网络的构建：利用步骤1中获得基因的表达数据矩阵X，通过BNFinder2工具进行调控网络构建，设置调控节点集合的数量，评分函数为BDE评分函数，得到一个概率网络，网络节点代表基因，有向边代表基因之间的调控关系；具体包括步骤2.1～步骤2.3：

步骤2.1：数据处理：根据X矩阵，确定基因个数及样本数量，并为每个基因设定潜在调控基因；

步骤2.2：设定基因调控网络结构的结构学习及参数学习，采用基于评分搜索的结构学习方法，遍历所有可能的结构，然后用BDE评分函数衡量各个结构，进而找出最好的结构；具体包括步骤2.2.1～步骤2.2.4：

步骤2.2.1：对于某个基因，首先获取该基因的潜在调控基因及其表达数据；

步骤2.2.2：为基因的每个潜在调控节点赋予权重；

步骤2.2.3：获取调控基因为空集的基因的BDE分数；

步骤2.2.4：遍历基因的潜在调控基因的集合，每次遍历结果即为该基因的一种调控结构，计算该结构的BDE分数，存入到一个堆栈中，遍历结束后，堆栈中留下的为评分最高的对应的调控基因的集合，堆栈堆的大小即为设定的调控节点集合的数量；

步骤2.3：根据网络的学习结果整合调控关系，得到基因调控网络结果；

步骤3：根据基因调控网络结果，进行节点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性，将BNFinder2得到的基因调控网络，通过Python的扩展包networkx计算各个基因的度中心性，接近中心性，中介中心性和特征向量中心性；

步骤4：取网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据作为为筛选出的相关基因。

2.根据权利要求1所述一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，所述网络节点的度中心性，接近中心性，中介中心性和特征向量中心性的前N个数据，N≤10，且矩阵X中的基因数量M>N。

3.根据权利要求1所述一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，所述步骤1.3中人工选择阈值I_sn，其中sn等于需要筛选出来的基因数量除以需要的轮数。

4.根据权利要求1所述一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，所述基准基因包括基因名称如下：AKT1、AKT2、APC、ARID1A、ARID1B、ARID2、ASXL1、BAP1、BRCA1、BRCA2、CASP8、CDH1、CDKN1B、CDKN2A、CCND1、MDM2、99Mb、ERBB2、ZNF217、FGFR1/ZNF703、MYC、GATA3、KRAS、MAP2K4、MAP3K1、MAP3K13、MIL2、MIL3、NCOR1、NF1、PIK3CA、PTEN、RB1、SETD2、SF3B1、SMAD4、SMARCD1、STK11、TBX3和TP53。

5.根据权利要求1所述一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，所述互信息的计算方法，如下公式：

其中，I(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的协方差矩阵的行列式。

6.根据权利要求1所述一种基于乳腺癌疾病的基因调控网络构建及分析方法，其特征在于，所述点中心性分析，包括节点的度中心性，接近中心性，中介中心性和特征向量中心性，用如下公式进行计算：

节点度中心性：在一个网络图G＝(V,E)中，节点v的度d_v指的是与v相连的E中边的数量，公式为：

其中，C_B(v)为节点度中心性，m为该节点的度，n为基因的数量；

所述接近中心性：

其中，C_cl(v)为接近中心性，n为基因的数量，d(v,u)为节点v和u之间的最短路径；

所述中介中心性：

其中，C_B(v)为中介中心性，σ(s,t|v)表示经过节点v的s到t的最短路径条数，σ(s,t)表示s到t的最短路径条数；

特征向量中心性：节点i的特征向量中心性是由以下等式定义的向量y的第i个元素：

Ay＝λy

其中，A是网络的邻接矩阵，特征值为λ，y为解上述方程得出的特征向量。