互信息關聯模型在基因調控網絡構建中的應用*

王銳 張獻志 曾麗華 陳活良 張潔瓊

(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.