沈威 ，鄭明，劉桂霞，邢翀，吳佳楠，周春光，周柚

(1. 吉林大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院，吉林 長春，130012；2. 北華大學(xué) 計算機(jī)科學(xué)技術(shù)學(xué)院，吉林 長春，132021)

隨著互補(bǔ)基因技術(shù)(cDNA)和微陣列技術(shù)的進(jìn)步，人們已經(jīng)能夠在整個基因組范圍內(nèi)進(jìn)行基因調(diào)控網(wǎng)絡(luò)的構(gòu)建。這些通過實(shí)驗獲得的數(shù)據(jù)提供了在一些特定時間點(diǎn)各種不同的生化環(huán)境下基因表達(dá)的、很有價值的描述。并且這些數(shù)據(jù)讓研究者能從基因的時序表達(dá)值下推斷基因調(diào)控網(wǎng)絡(luò)。盡管如此，需要的基因調(diào)控網(wǎng)絡(luò)是存在于基因組中的成千上萬的基因所構(gòu)成的網(wǎng)絡(luò)，極其復(fù)雜。所以，在系統(tǒng)生物學(xué)中一項極具挑戰(zhàn)的工作就是從各種數(shù)據(jù)尤其是時序數(shù)據(jù)中推斷出基因調(diào)控網(wǎng)絡(luò)。在可觀測系統(tǒng)和從算法中獲得的網(wǎng)絡(luò)之間是存在差異的，這種差異導(dǎo)致了算法精確度的差異，目前很難提供一個高度精確的算法來實(shí)現(xiàn)對任意實(shí)驗數(shù)據(jù)的網(wǎng)絡(luò)推斷。盡管如此，研究更精確的網(wǎng)絡(luò)推斷模型仍然是生物信息學(xué)的研究熱點(diǎn)。許多方法被提出并應(yīng)用在基因調(diào)控網(wǎng)絡(luò)推斷上，這些方法包括布爾網(wǎng)絡(luò)[1]、貝葉斯網(wǎng)絡(luò)[2-3]和微分方程模型[4]等。布爾網(wǎng)絡(luò)提供了一種網(wǎng)絡(luò)結(jié)構(gòu)關(guān)系的概念層方法，網(wǎng)絡(luò)一般規(guī)模龐大，常常包含成百上千的基因，但是，網(wǎng)絡(luò)中節(jié)點(diǎn)的狀態(tài)只有開和關(guān)2種狀態(tài)。貝葉斯網(wǎng)絡(luò)是一種能夠描述網(wǎng)絡(luò)關(guān)系的概率模型，網(wǎng)絡(luò)模型中每條邊都存在1個概率，貝葉斯網(wǎng)絡(luò)的圖形構(gòu)成的是1個有向無環(huán)圖，節(jié)點(diǎn)一般包括幾十個節(jié)點(diǎn)，規(guī)模適中。微分方程模型則是一種高度精確的網(wǎng)絡(luò)構(gòu)建模型，以變化率的形式表示出各個基因的影響關(guān)系，是一個動態(tài)方程組。 這種方式構(gòu)成的網(wǎng)絡(luò)往往只包含幾個基因，一般在20個以內(nèi)，所構(gòu)成的網(wǎng)絡(luò)在數(shù)學(xué)上最精確，即規(guī)模小而精。但是無論采用哪種方法，仍然存在精度不高的問題，因此，需要一個更好的方法來進(jìn)一步改進(jìn)微分方程模型。改進(jìn)策略之一是改變方程的形式，如把方程中某個參數(shù)進(jìn)行分解或改變，更可以把微分問題化為積分問題[5]等，但這種策略在本質(zhì)上并沒有改變精度。另一種改進(jìn)策略是引入奇異值分解(SVD)[6]的方法。在運(yùn)算過程中把方程組中表達(dá)數(shù)據(jù)所構(gòu)成的矩陣或其轉(zhuǎn)置矩陣進(jìn)行分解，之后進(jìn)行矩陣運(yùn)算并得到符合奇異值分解規(guī)則的解即權(quán)值矩陣作為最終網(wǎng)絡(luò)結(jié)果。在多個基因的實(shí)驗數(shù)據(jù)條件下，方程組的解往往并不唯一，而奇異值分解求出的是惟一解，這個解可能在無數(shù)解集中并不是最優(yōu)解。針對這個問題，本文作者在通過奇異值分解求得的特解基礎(chǔ)上求出符合奇異值分解規(guī)則的通解即解的集合，并把這個通解作為一個取值范圍，結(jié)合改進(jìn)形式的微分方程模型策略用優(yōu)化算法求解和優(yōu)化網(wǎng)絡(luò)的權(quán)值矩陣并得出網(wǎng)絡(luò)的最終結(jié)果。

1 算法的數(shù)學(xué)框架

用于推斷網(wǎng)絡(luò)的方法可以分為2步。首先，要使用奇異值分解方法來構(gòu)建一個解的集合。這個集合是由符合方程組規(guī)律的解集所組成的，如果規(guī)模很大，一般提供通解。然后，在這個集合或標(biāo)準(zhǔn)下進(jìn)行模型的優(yōu)化，使用線性回歸分析進(jìn)行微分方程的求解，并最終得到所求網(wǎng)絡(luò)。

為了簡化，僅對穩(wěn)態(tài)系統(tǒng)或者接近穩(wěn)態(tài)的系統(tǒng)進(jìn)行分析。構(gòu)建網(wǎng)絡(luò)的數(shù)學(xué)模型采用的是微分方程模型。這個模型可以簡化為每個方程均為式(1)所示形式的方程所組成的方程組：

其中，xi(t)為基因i在時刻t的濃度；N為網(wǎng)絡(luò)所包含基因的總數(shù)；wij為基因j對基因i進(jìn)行影響的權(quán)值參數(shù)，這個值是正數(shù)時表示正調(diào)控，負(fù)數(shù)時表示負(fù)調(diào)控，0表示在一定閾值范圍內(nèi)可以視為不發(fā)生相互作用；bi表示在沒有調(diào)控輸入的情況下基因 i的變化，也包括了外界對基因i的刺激所產(chǎn)生的表達(dá)濃度的變化。

在實(shí)驗中，可以提供一些恒定的外界刺激，即bi是恒定的。為了方便，更多的時候設(shè)定bi為0，即不施加任何外界刺激。然后，對N個基因在T個時間段進(jìn)行測量，即得到所需數(shù)據(jù)，數(shù)據(jù)組織形式如下：

其中：x為基因表達(dá)濃度；下標(biāo)表示基因標(biāo)識，上標(biāo)表示測試時間序列；共有N個基因、T個時間測試點(diǎn)。根據(jù)式(2)，式(1)可以寫成式(3)所示形式：

為了求出權(quán)值矩陣W，往往把方程組進(jìn)行變形和改進(jìn)，這可以在一定程度上提高精度，但要從根本上提高精度需要進(jìn)行奇異值分解?？梢詫矩陣進(jìn)行奇異值分解：

其中：矩陣A是一個不必滿足滿秩條件的對角矩陣。為了計算方便，可以在計算過程中進(jìn)行調(diào)整，使得矩陣A中對角線上的值是按照大小進(jìn)行排序的。即非零值第1個值最大，其他值遞減。后面的幾個對角線上的值為0，非對角線上的值全為0。對角線上零值的個數(shù)可以為 0，即此矩陣可以為滿秩矩陣。在式(4)中，奇異值分解得到的結(jié)果U與V都是正交矩陣，且都滿足：

式中：E為單位矩陣。

對式(3)兩邊同時減BN×T，然后，方程兩邊同時右乘XN×T的逆矩陣，并把式(4)代入式(3)，得到：

在式(6)中，并不是每個矩陣都存在逆矩陣，特別是矩陣 U，這里的逆矩陣指的就是廣義逆矩陣[7]。在廣義逆矩陣定義下，任意一個矩陣都存在其逆矩陣。但其廣義逆矩陣往往不是方陣，更不是奇異陣。

由式(6)得到基因調(diào)控網(wǎng)絡(luò)構(gòu)建所需的權(quán)值矩陣W。研究者們一般就把此特解作為最終的網(wǎng)絡(luò)結(jié)果，但該結(jié)果并不精確，因為符合式(6)條件的值很多，數(shù)量級很大，甚至無窮多解，即N＞＞T，這樣所求的值其實(shí)并不一定是最優(yōu)解。在計算過程中，最好能夠盡可能多的得到解集或者取值范圍，并使用優(yōu)化算法進(jìn)行計算。本文算法在特解的基礎(chǔ)上，構(gòu)造下式來求通解：

式中：C為任意常數(shù)。

構(gòu)造了這個通解后，就可以得到所求解的取值范圍，從而可以修改微分方程模型，并計算優(yōu)化，在此范圍內(nèi)求得所求網(wǎng)絡(luò)。

從式(1)出發(fā)對權(quán)值矩陣 W 進(jìn)行變形。首先需要定義此算法的能量函數(shù)，即誤差最小化函數(shù)。此誤差為所求值和觀測值之間的差：

其中：rik為式(1)的殘基；vik為rik的1個權(quán)值，則

從式(9)可以看出：rik符合殘基概念。本文所要做的就是在式(7)所示的網(wǎng)絡(luò)權(quán)值取值范圍內(nèi)使殘基最小化。得到所要進(jìn)行的目標(biāo)之后，需要把方程進(jìn)一步改進(jìn)，以便提高效率，較少誤差。在式(8)和(9)中，所需計算的只有wij，而網(wǎng)絡(luò)中共有N個基因且權(quán)值矩陣W為一方陣，則所需計算的參數(shù)為N個。一旦基因數(shù)過于龐大，則所需計算的參數(shù)個數(shù)將隨指數(shù)級增長，不適合求解，因此，需要把權(quán)值矩陣W進(jìn)一步優(yōu)化并分解，這樣權(quán)值w將變?yōu)椋?/p>

經(jīng)過式(10)的轉(zhuǎn)化之后，權(quán)值矩陣W中的值只與t和x這2個變量相關(guān)，從而使得計算非常簡便。

優(yōu)化算法方面，本文采用粒子群算法[8]和遺傳算法[9]相結(jié)合的方式，通過最小化式(8)，求得最小值，然后驗證是否符合式(7)的方式，求得最終的網(wǎng)絡(luò)。

2 算法的實(shí)現(xiàn)

本文算法的實(shí)現(xiàn)部分用Java語言編寫，所用環(huán)境為MyEclipse 8.0版本。算法的流程步驟如下。

(1) 讀取基因表達(dá)數(shù)據(jù)，其格式符合國際表達(dá)數(shù)據(jù)標(biāo)準(zhǔn)(MIAME[10])。

(2) 對已經(jīng)讀取的基因表達(dá)數(shù)據(jù)所構(gòu)成的矩陣運(yùn)用式(4)進(jìn)行奇異值分解并得到式(4)中的U與V。

(3) 運(yùn)用式(6)求得所需計算的權(quán)值矩陣W的特解W特，并根據(jù)式(7)求得W的通解W通。

(4) 隨機(jī)生成大量權(quán)值構(gòu)成粒子群算法和遺傳算法的初始群體。并根據(jù)步驟(3)得到的取值范圍淘汰不符合要求的隨機(jī)權(quán)值。淘汰的權(quán)值重新生成，直到所有的隨機(jī)權(quán)值符合W通的條件為止。

(5) 運(yùn)用式(10)分解權(quán)值，將權(quán)值降維為只有2個變量的方程組。并運(yùn)用遺傳算法和粒子群算法(PSO)最小化式(8)，得到最小值。

(6) 若步驟(5)計算得到的值符合 W通的條件則算法終止，否則進(jìn)入步驟(5)。

算法流程圖如圖1所示。

圖1 本文算法流程圖Fig.1 Flow chart of proposed algorithm

3 實(shí)驗結(jié)果與分析

3.1 測試方法

為了測試本文算法的效果，進(jìn)行比較實(shí)驗。比較對象為單純使用式(1)的基本微分方程模型(方法 1)和代入式(10)的改進(jìn)型微分方程模型(方法2)及本文模型(方法3)。這3種模型的比較實(shí)驗將分為模擬數(shù)據(jù)實(shí)驗和真實(shí)數(shù)據(jù)實(shí)驗2部分。

采用 Bansal等[11]所提出來的正預(yù)測值(PPV)和敏感度(SE)評價算法優(yōu)劣。

3.2 模擬數(shù)據(jù)實(shí)驗

數(shù)據(jù)來源采用生成無標(biāo)度網(wǎng)絡(luò)[12]后微分方程逆向反推基因表達(dá)數(shù)據(jù)的方法進(jìn)行計算獲得。與隨機(jī)生成的網(wǎng)絡(luò)相比，無標(biāo)度網(wǎng)絡(luò)更符合基因調(diào)控網(wǎng)絡(luò)的生物學(xué)現(xiàn)象。在生成的無標(biāo)度網(wǎng)絡(luò)中，每個節(jié)點(diǎn)的連接數(shù)都接近于冪率分布 p(k)～k。本文用增加重導(dǎo)向自增算法[13]進(jìn)行構(gòu)建。參數(shù)r設(shè)定為3，節(jié)點(diǎn)數(shù)N為20。根據(jù)這種算法和參數(shù)構(gòu)建的網(wǎng)絡(luò)如圖2所示。在圖2中，原點(diǎn)為節(jié)點(diǎn)0，逆時針依次為節(jié)點(diǎn)1，2，…，20。方塊的為有向線段末端，指向被調(diào)控節(jié)點(diǎn)。網(wǎng)絡(luò)生成以后，用式(1)來進(jìn)行數(shù)據(jù)推斷，所有的bi均設(shè)置為0。數(shù)據(jù)的時間點(diǎn)設(shè)定為 50。隨機(jī)噪聲參數(shù)設(shè)定為-0.5～0.5之間。

由圖2所示的圖模擬出基因表達(dá)數(shù)據(jù)后，經(jīng)過運(yùn)算得到結(jié)果如圖3所示。在圖3中，所表示的含義與圖2的一樣，不同的是圓圈代表的是自我調(diào)控關(guān)系。

為了避免結(jié)果巧合，運(yùn)行程序100次，并把N的取值定為 20。計算后取平均值并評估算法的優(yōu)劣。3種方法的比較結(jié)果如表1所示。

從表1可知：本文算法(即方法3)的正預(yù)測值和敏感度明顯高于其他2種算法的結(jié)果，其運(yùn)行時間也比其他2種方法的運(yùn)行時間短。

圖2 參數(shù)r=3，N=20的無標(biāo)度網(wǎng)絡(luò)圖Fig.2 Scale-free network with parameter r=3, N=20

圖3 根據(jù)圖2所生成的數(shù)據(jù)所計算的基因調(diào)控網(wǎng)絡(luò)圖Fig.3 Gene regulatory network calculated with proposed algorithm from data simulated by Fig. 2

表1 3種算法模擬數(shù)據(jù)效果比較Table 1 Compared results of three kinds of algorithm in simulation data

方法2比方法1在運(yùn)行時間上更有優(yōu)勢的原因在于方法2把方法1中的N2個參數(shù)化簡為只有2個參數(shù)，使得運(yùn)算時間大大縮短。本文方法為方法2提供了一個很好的候選解的集合，減少了方法2尋優(yōu)的空間和運(yùn)行時間，但是，由于奇異值分解運(yùn)算本身消耗了一定的時間，因此，與方法2的運(yùn)行時間相比，本文方法的運(yùn)行效率略有提高。

3.3 真實(shí)數(shù)據(jù)實(shí)驗

真實(shí)數(shù)據(jù)的計算采用的是 GEO[14]基因表達(dá)數(shù)據(jù)庫上下載的基因表達(dá)數(shù)據(jù)，這里采用的是經(jīng)典的酵母基因表達(dá)數(shù)據(jù)，即采用GDS38數(shù)據(jù)的子集，其表達(dá)的結(jié)果圖可以從日本KEGG[15]數(shù)據(jù)庫中的pathway里得到。這個數(shù)據(jù)包括3個子集，這3個子集包括3種實(shí)驗條件：Alpha，elu和cdc15，它們分別包括18，14和24個時間點(diǎn)。抽取出來的1個網(wǎng)絡(luò)子集如圖4所示。經(jīng)過本文算法計算，對應(yīng)圖4所生成的基因調(diào)控網(wǎng)絡(luò)如圖5所示。在圖5中，所有的含義與圖2中的一樣，原點(diǎn)作為第1個初始點(diǎn)，圓圈表示自調(diào)控，方塊表示被調(diào)控段。

圖4 KEGG數(shù)據(jù)庫的pathway中酵母基因調(diào)控網(wǎng)絡(luò)的1個子集圖Fig.4 A part of yeast cell cycle pathway available from KEGG

圖5 用本文算法算得的基因調(diào)控網(wǎng)絡(luò)圖(對應(yīng)圖4)Fig.5 Gene regulated network calculated by proposed algorithm (according to Fig. 4)

根據(jù)式(1)和(10)的微分方程模型，運(yùn)行100次后取平均值，最終結(jié)果如表2所示。

從表 2可以看出：本文算法(即方法 3)比其他 2種算法更優(yōu)越。從運(yùn)行時間上來看，本文算法比其他2種算法相比更小。與模擬實(shí)驗結(jié)果相比，優(yōu)勢不是很明顯。

表2 3種算法運(yùn)算真實(shí)數(shù)據(jù)結(jié)果比較Table 2 Compared results of three kinds of algorithm in real data

4 結(jié)論

(1) 提出了一種具有通解的構(gòu)建基因調(diào)控網(wǎng)絡(luò)算法。這種新的融合算法不僅在運(yùn)算過程中提供了一個候選解的集合，使得運(yùn)算過程中免去了非解的不必要計算，縮短了運(yùn)算時間，同時，使得運(yùn)算結(jié)果更精確。通過與其他2種傳統(tǒng)微分方程模型算法的對比實(shí)驗，證實(shí)了本文算法的效率較高，符合預(yù)期結(jié)果。

(2) 本文算法的目標(biāo)是從有限的數(shù)據(jù)中得到最接近于真實(shí)的基因調(diào)控網(wǎng)絡(luò)的結(jié)果，提高實(shí)驗水平。但是，本文的基因調(diào)控網(wǎng)絡(luò)算法整體的精度，還需進(jìn)一步提高，這有待進(jìn)一步研究。

[1] D’haeseleer P, Liang S, Somogyi R. Genetic network inference:from co-expression clustering to reverse engineering[J].Bioinformatics 2000, 16: 707-726.

[2] Dojer N, Gambin A, Mizera A, et al. Applying dynamic Bayesian networks to perturbed gene expression data[J]. BMC Bioinformatics, 2006, 7: 249.

[3] Beal M J, Falciani F, Ghahramani Z, et al. A Bayesian approach to reconstructing genetic regulatory networks with hidden factors[J]. Bioinformatics, 2005, 21: 349-356.

[4] Bansal M, Gatta G D, di Bernardo D. Inference of gene regulatory networks and compound mode of action from time course gene expression profiles[J]. Bioinformatics, 2006, 22:815-822.

[5] Eugene N, Emmanuel B. Regulatory network reconstruction using an integral additive model with flexible kernel functions[J].BMC Systems Biology, 2008, 2: 8.

[6] Nelson P A, Kahana Y. Spherical harmonics, singular-value decomposition and the head-related transfer function[J]. Journal of Sound and Vibration, 2001, 239(4): 607-638.

[7] Liang M L, Dai L F. The left and right inverse eigenvalue problems of generalized reflexive and anti-reflexive matrices[J].Journal of Computational and Applied Mathematics, 2010, 234:743-749.

[8] Brits R, Engelbrecht A P, van den Bergh F. Locating multiple optima using particle swarm optimization[J]. Appl Math Comput,2007, 189(2): 1859-1883.

[9] Chen D Y, Chuang T R, Tsai S C. Jgap: A java-based graph algorithms platform[J]. Software Pract Exper, 2001, 31(7):615-635.

[10] Tim F R, Philippe R S, Paul T S, et al. A simple spreadsheet-based, MIAME-supportive format for microarray data: MAGE-TAB[J]. BMC Bioinformatics, 2006, 7: 489.

[11] Bansal M, Belcastro V, Ambesi-Impiombato A, et al. How to infer gene networks from expression profiles[J]. Mol Sys Biol,2007, 3: 78.

[12] Barabási A L, Albert R. Emergence of scaling in random networks[J]. Science, 1999, 286: 509-512.

[13] Silva A C, da Silva J K L, Mendes J F F. Scale-free network with Boolean dynamics as a function of connectivity[J]. Phys Rev E,2004, 70(6): 66140-66147.

[14] Wilhite S E, Barrett T. Strategies to explore functional genomics data sets in IVCBI’s GEO database[J]. Methods Mol Biol, 2012,802: 41-53.

[15] KEGG Cell cycle-yeast-Saccharomyces cerevisiae[EB/OL].[2010-04-13].http://www.genome.jp/kegg/pathway/sce/sce0411 1.html.