互信息关联模型在基因调控网络构建中的应用*

王锐 张献志 曾丽华 陈活良 张洁琼

(1.第四军医大学卫生勤务学教研室，陕西 西安 710032；2.第四军医大学放射医学教研室，陕西 西安 710032)

学习与记忆能力是人类高级精神活动的特征与基础。学习与记忆机制的研究越来越受到科学领域的广泛关注，这对于揭示大脑的奥秘、并最终认识生命的本质有着重要的意义[1]。相对于其他各种组织器官，人的大脑在功能、结构和基因的表达上表现出更大的复杂性。胎儿海马发育和基因的表达与成人及病变的海马有着极为密切的联系，为了在全基因组水平上系统地研究胎儿阶段中枢神经系统发育情况，香港城市大学科研人员[2]利用基因芯片检测了人脑小脑组织在胎儿的不同发育时期，12060个基因的表达情况，得到了大量的基因表达数据。

发育的最终状态是与发育有关的基因在时间和空间上相互调控的结果[3]。利用时间序列的高通量基因芯片表达数据，通过合适的数学模型，挖掘出基因之间相互调控的网络，这对认识未知基因的功能、整个神经系统的工作原理等重大问题是有积极作用的。目前，已有较多的数学模型应用于调控网络的建模及构建，如灰色理论、Bayesian网络等等[4]，大多数模型都没有考虑基因表达在时-空层面上的调控关系，或者没有考虑基因表达在不同时间点上权重关系，这为深入挖掘基因调控的时-空机制带来了一定的困难。因此，在计算各基因之间的互信息相关系数之前，首先将时间点进行加权，这个方法能够提高调控网络构建的精度[5]。本研究借助互信息关联模型，充分考虑了基因表达时间点的相互关系，利用互信息相关系数来描述发育过程中基因之间的相互关系。

1 材料与方法

1.1 小脑组织基因的选择 我们选择了小脑组织中7个时间点都有表达数据的237个基因。将237个基因投射到Gene Ontology Consortium[6](GO数据库)中，着重选择了生物过程中具有发育功能的43个基因。顺序排列为：DPYSL3, VEGFB, HIF1A, MEIS2, MAP1B, GAS7, NOS3, NNAT, EIF3S2, CNTN2,STMN1, HSF4, MSH6, DOK5, ACTG1, MEF2C, APOD, FEZ1, WASF1, ZIC2, PRKCB1, NEUROD1, APC, H2AFX, E2F5, PLP1, MBD3, BLR1, NRXN1, CRYAB, ADD2, PTN, CLU, NRXN3, DKK3, MET, RPS29, FHL1, HMGB1, NELL2, DACH, RELN, YWHAH。将43个基因投射到20个功能模块中，见表1。

1.2 基因表达各时间点的权重确定 不同发育阶段的基因表达在各时间点上权重是不同的。随着时间的推移，基因表达强度越弱，因此，需要对7个时间点的表达值进行加权。加权处理之前要对基因表达数据进行正则化[7]，采用公式(1)。其正则化的结果，见表2。

(1)

(2)

(3)

表1 43个基因的功能分布Table 1 The distribution of 43 genes’ function

表2 基因表达值正则化结果Table 2 Regularization result of gene expression values

2 结果

2.1 相关系数是用以反映变量之间相关关系密切程度的统计指标。那么，定性变量间的相关度分析可以用信息量来描述定性变量间相关程度的强弱[8]。

式中：p(xi)为发生事件xi的概念，n为可能发生的事件(状态)总数。

2.3 互信息相关系数的计算 互信息相关系数可通过信息熵信息量来确定，也就是信息熵相关系数[11]。它用于描述基因向量genei=(gi1,gi2,…,gin)和向量genej=(gj1,gj2,…,gjn)有多少共同的信息，它的性质类似于非线性相关系数，但更具有不确定性,即当用不同的方法离散连续变量时就会得到不同的结果[12]。

根据各时间点的属性权重，计算出加权后的表达值，采用局部加权线性回归标准化法，对基因表达值进行标准化处理[13]，见公式(5)。对标准化后的数据采用10等分区间[14]。

(4)

(5)

log2是以2为底的对数，p(x)为基因表达值位于第x分位数上的概率值。以基因DPYSL3为例，对各时间点加权后的表达值进行标准化，结果分别为0.0275、0.9519、0.7938、0.6014、1、0.6289和0，将0～1按 10等分进行区间划分，区间之间大小为0.1。结果当中有两个值在第1个区间[0,0.1]，三个值在2区间，一个值在第10个区间，通过公式(6)计算，可信息熵H1=1.94。同理，可得到基因VEGFB的信息熵H2=2.25。上述的两个基因的互信息MI(g1,g2)可通过公式(7)，得到MI(g1,g2)为1.58。通过公式(8)可得到上述两基因相关系数r为0.70。

MI(g1,g2)=H(g1)+H(g2)-H(g1,g2)

(6)

(7)

建立信息熵相关系数矩阵R，相关系数rij表示两个基因之间的调控关系，R可表示为：

(8)

2.4 应用MATLAB 2015a软件可得到43个基因之间调控关系的强弱程度[15]。本研究中对基因之间的作用直接通过相关系数值的形式体现，数值越大强度越强。为了便于相互的比较，我们采用灰度图表示调控矩阵，见图1。

图1 调控矩阵的图像显示Figure 1 Visualization of controlling matrix

注：以灰度表示其调控的强度。白色为调控最大值，黑色为调控最小值

3 讨论

互信息作为相关分析的度量，其最大优势在于能有效刻画变量之间的非线性关系[16]。在大数据相关分析中，最具影响力的研究成果是Reshef等人于2011年发表在《 Science》上的论文“Detecting novel associations in large data sets”[17]，该研究通过互信息定义了两个变量之间的互信息相关系数，用来衡量两个变量之间的相关性。同时，通过对多属性群决策方法的研究，信息熵可用来度量信息量的大小，也是信息有用程度的一种表现形式，采用了以信息熵确定属性权重的方法[18]，客观性较强。 多物种生态系统动态发展理论认为[19]：每种生物都生活在一定的生态系统中，任何一种生物的数量的消长，都会对其他生物产生影响，会影响它所在的生态系统，即包括互惠、竞争和共存三种关系：①互惠关系：即基因之间表达的互相促进。②竞争关系：基因的表达会抑制另一些基因的表达。③共存关系：基因之间的表达没有直接关系，仅共存于同一个调控网络中。在模型建立过程中，我们采用数学函数当中的样条函数进行插值[20]，描述其表达趋势。在网络可视化方面，我们采用了带有色彩梯度变化的灰度图像，该图像可以将基因表达方法很好地表现出来[21]。该方法避免了阈值的选择问题，可以根据实际情况灵活选择。

利用人类胎儿小脑组织发育的基因表达数据，基于非参数统计理论提出了互信息和熵相关系数相结合分析方法，得到了与发育有关基因的相互调控情况。例如31号基因(ADD2)与6号基因(GAS7)、31号基因(ADD2)与38号基因(FHL1)信息熵相关系数分别为0.88和0.98，相关程度很强，通过相关生物学资料可知，ADD2对GAS7和FHL1的表达为抑制作用。ADD2基因与分裂组织发育、形态发生、器官发育、基部发育和促发发育等功能有关，而GAS7、FHL1等几种基因也与形态发生和器官发育有关。

4 结论

本研究从神经科学的研究热点和前沿切入，具有一定的生物与方法学意义。基于互信息关联网络描述基因表达调控的复杂关系。该方法有较广泛的适用性，可为进一步探索生物信息学方法提供理论支持。

[1]郭波涛，李辉智，王文昌，等.酵母基因调控网络的微分方程模型研究[J].中国卫生统计,2006,23(2):129-133.

[2]Ashburner M, Ball C, Blake J,etal. Gene Ontology: tool for the unification of biology [J]. NatGenet，2000,25:25-29.

[3]陈云翔，董骁雄，项华春，等.基于信息熵的群组聚类组合赋权法[J]. 中国管理科学, 2015, 23(6):142-146.

[4]Guo Hongyu.Research on term weighting algorithm based on information entropy theory[J]. Computer Engineering and Applications, 2013, 49(10): 140-146.

[5]易东，杨梦苏，李辉智，等.相关分析在建立基因调控网络中的应用[J].中国卫生统计,2003,20(3):144-146.

[6]J Narula，CJ Williams，A Tiwari，etal.Mathematical Model of a Gene Regulatory Network Reconciles Effects of Genetic Perturbations on Hematopoietic Stem Cell Emergence[J]. Developmental Biology, 2013, 379(2):258-269.

[7]Zhang X, Zhao XM, He K,etal.Inferring gene regulatory networks from gene expression data by path consistency algorithm based on conditional mutual information[J]. Bioinformatics,2012,28(1):98-104.

[8]Qian J, L in J, Luscombe NM,etal. Prediction of regulatory networks:genomeOwide identification of transcrip tion factor tar2 gets from gene expression data [J]. Bioinformatics, 2003, 19(15 ): 1917-1926.

[9]Y Wang，T Joshi，XS Zhang,etal. Inferring gene regulatory networks from multiple microarray datasets[J]. Bioinformatics. 2006,22(19):2413-2420.

[10] 杨英杰，李红燕，谢建平，等. MATLAB 7.X生物信息工具箱的应用——基因芯片分析(三) [J].现代生物医学进展, 2008, 8(4):704-708.

[11] Zhou X，Zhang FM, Hui XB,etal. Group decision-making method based on entropy and experts cluster analysis[J]. Journal of Systems Electronics,2011, 22(3):468-472.

[12] S Feizi,D Marbach,M Médard,etal. Network deconvolution as a general method to distinguish direct dependencies in networks[J]. Nature Biotechnology, 2013, 31(8):726.

[13] 黎妍，张晓飞，易鸣，等. 基因调控网络的边预测[J]. 数学物理学报，2015，35(5):1018-1024.

[14] 张宇镭，党琰，贺平安．利用Pearson相关系数定量分析生物亲缘关系[J].计算机工程与应用, 2005, 41(33):79-82.

[15] 丁晶，赵永龙.以互信息为基础的广义相关系数[J].四川大学学报(工程科学版)， 2002, 34(3):1-5.

[16] B Barzel.Network link prediction by global silencing of indirect correlations [J]. Nature Biotechnology, 2013, 31(8):720-725.

[17] AK Smilde,HAL Kiers,S Bijlsma,etal.Matrix correlations for high-dimensional data [J]. Bioinformatics,2009,25(3):401-405.

[18] TT Cai, CH Zhang, HH Zhou. Optimal rates of covariance matrix estimation [J]. Annals of Statistics,2010,38(4):2118-2144.

[19] Hood L, Heath J R, Phelp sHoodM E,etal.Systems biology and new technologies enable predictive and preventative medicine[J]. Science, 2004, 306(5696) : 640-643.

[20] Wyrick J J， Young R A. Deciphering gene expression regulatory networks[J]. Current Opinion in Genetics & Development, 2002, 12(2):130-136.

[21] Steuer R, Kurths J， Daub C O,etal. The mutual information: detecting and evaluating dependencies between variables[J].Bioinformatics, 2002, 18: 16-25.