基于改進SEM算法的基因調控網絡構建方法

摘 要:動態貝葉斯網絡(DBN)是基因調控網絡的一種有力建模工具。貝葉斯結構期望最大算法(SEM)能較好地處理構建基因調控網絡中數據缺失的情況，但SEM算法學習的結果對初始參數設置依賴性強。針對此問題，提出一種改進的SEM算法，通過隨機生成一些候選初始值，在經過一次迭代后得到的參數中選擇一個最好的初始值作為模型的初始參數值，然后執行基本的SEM算法。利用啤酒酵母細胞周期微陣列表達數據，構建其基因調控網絡并與現有文獻比較，結果顯示該算法進一步提高了調控網絡構建的精度。

關鍵詞:基因調控網絡; 動態貝葉斯網絡; 貝葉斯結構期望最大化算法

中圖分類號:TP181

文獻標志碼:A

文章編號:1001-3695(2010)02-0450-03

doi:10.3969/j.issn.1001-3695.2010.02.012

Method for modeling gene regulation network based onimproved structure expectation maximization algorithm

GE Ling-ling， WANG Hao， YAO Hong-liang

(College of Computer Science Technology， Hefei University of Technology， Hefei 230009， China)

Abstract:Dynamic Bayesian network(DBN) is a powerful moedling tool for gene rugulation network.Missing data in building gene regulation network is better dealt with SEM(Bayesion structure expectation maximization) algorithm，however，the result of learning by SEM algorithm has strong dependence on the initial parameters. This paper proposed an improver SEM algorithm， which randomly generated a number of candidate initial parameters and selected the best parameter as whole model’s initial parameter to execute basic SEM algorithm after a iterative process.Comparing gene regulation network constructed with yeast cycle gene expression data by improved SEM algorithm with existing literature，the result indicates further improve the accuracy of constructing regulation network.

Key words:gene rugulation network; dynamic Bayesian network(DBN); Bayesian structure expectation maximization algorithm

0 引言

基因調控網絡的研究目的是期望從系統的角度全面揭示基因組的功能和行為。目前構建基因調控網絡主要有布爾網絡[1]、微分方程[2]和貝葉斯網絡[3]等方法，這些方法都在不同層次上對真實的調控網絡進行了抽象化。布爾網絡是定性研究基因調控網絡，而微分方程又是通過精細的數學分析來量化描述生物過程，但缺乏抗噪聲能力，計算量大，魯棒性能不佳。貝葉斯網絡模型是這兩個極端的折中。利用貝葉斯網絡構建基因調控網絡是目前生物信息學研究的熱點[3]。

根據處理的基因表達數據類型的區別，貝葉斯網絡分為靜態和動態貝葉斯網絡兩種方法。前者適用于處理無時序信息的芯片表達數據，后者適用于處理有時序信息的數據。動態貝葉斯網絡在考慮時間因素后，通過劃分時間點，可以對生物過程中反饋等循環調控進行描述，克服了靜態貝葉斯不能解決有向有環圖的缺陷。

在構建基因調控網絡的實驗中所采用的微陣列數據噪聲大且有很多缺失值，這種缺失值經常會對結果的精確性有影響。然而目前用于基于貝葉斯網絡構建基因調控網絡的貪心搜索GS、 MWST和K2等網絡學習算法，常常被用來處理完備數據下的結構學習[4]。Friedman等人[5]將EM思想引入到丟失數據情況下的結構學習，借鑒參數學習的EM算法，提出貝葉斯結構期望最大(Bayesian structure expectation maximization，SEM)算法。該算法在一定程度上提高了學習效率，但是存在學習精度低、對初始參數值依賴的缺點。

本文對SEM算法的初始值θ0的選取進行了改進，在執行基本SEM算法之前先隨機生成一些候選初始參數，對這些初始參數分別執行EM算法的一次迭代，在迭代得到的值中用最大評分函數進行最優選擇，把使當前網絡評分最大的參數設置為整個算法的初始值θ0，提高了算法的準確性。 將此方法應用于啤酒酵母細胞周期的基因表達譜數據中，構建基因調控網絡。與現有文獻比較，結果表明改進后的算法進一步提高了構建調控網絡的精度。

1 動態貝葉斯網絡

動態貝葉斯網絡(DBN)是建立在時間序列數據集上的一種貝葉斯網絡，把不同時間點上的隨機變量區分開來，當做不同的隨機變量來構建貝葉斯網絡，從而避免了出現自環的情況。

從靜態貝葉斯轉變到動態貝葉斯，需要作一些假設和簡化處理。假設條件如下:

a)在一個有限的時間內條件概率變化過程對于所有t是一致平穩的。

b)動態概率過程是馬氏的，即未來時刻的概率只與當前時刻有關而與過去時刻無關。

基于以上假設，建立在隨機過程時間軌跡上的聯合概率分布的動態貝葉斯由兩部分組成:

a)一個先驗網絡B0。定義在初始狀態x0上的聯合概率分布。

b)一個轉移網絡B→。定義在變量xt與xt+1上的轉移概P(xt|xt+1)(對所有t都成立)。設X={X1，…，Xn}是動態貝葉斯網絡的隨機變量集，xit表示變量xi在t時對應的隨機變量，則動態貝葉斯網絡在X={X1，…，Xn}上的聯合概率分布可以表示為

PB(X1，…，Xn)=PB0(X0)∏T-1t=0PB→(Xt|Xt-1)(1)

靜態貝葉斯網絡無法描述圖1(a)中如X1→X2→X5→X1的環狀反饋結構，但是在生物過程中包括很多像反饋這樣的循環調控過程。動態貝葉斯網絡考慮時間因素后，通過劃分時間點，可以將上述反饋調控作如圖1(b)中X1(t)→X2(t+2)→X5(t+3)→X1(t+4)形式的描述。

2 改進的SEM算法

2.1 基本SEM算法

SEM算法是通過調用貝葉斯網絡的推理算法填充數據集D，將不完備數據下的結構學習問題轉換為較容易解決的完備數據下的結構學習問題。SEM算法主要分為結構搜索和參數學習兩步。進行結構搜索時， SEM算法使用期望充分統計因子代替不存在的充分統計因子，這樣可使打分函數具有可分解形式，再進行局部搜索，以試圖發現得分更高的網絡結構;然后在選定的網絡結構上尋找使得分最大的參數。該方法能夠在一定程度上提高學習效率，并為具有缺省數據的貝葉斯結構學習提供一個框架。

其基本過程是:在每次迭代中，首先利用參數EM算法來實現參數最大化，然后利用BIC評分來進行模型選擇。其中，給定初始模型M0，模型Mi(i>0)是在EM算法第i次迭代中生成，且模型序列M0，…，Mn中的模型個數與迭代的次數相同。

算法可描述如下:

隨機選定模型結構和參數M0及θ0

loopn=0，1，…，直至算法收斂

{loop l=0，1，…，直至算法收斂或l=lMAX

{ 使得θn，l+1=argmaxθQ(θ:Mn，θn，l);}

θn+1=θn，l+1;

將Score(M:θn，D)和Score(M:θn+1，D)中得分最大的模型賦給Mn+1;將Mn+1模型的參數賦給θn+1，0}

其中:Score(M:θn，D)和Score(M:θn+1，D)為BIC評分函數;argmaxθ Q(θ:Mn，θn，l)為求似然函數期望最大化。

由上述算法描述可以看出，SEM算法的執行結果對初始參數有很強的依賴性，不好的初始值會導致學習過程的循環次數增加，降低算法的時間性能和結果的學習精度。因此，怎樣選取到一個好的初始參數對最后的網絡結果有較大的影響。針對這個問題，本文提出了改進的SEM算法。

2.2 改進的SEM算法

改進的SEM算法是針對初始值θ0的選取進行的改進，其主要思想是在執行基本的SEM算法之前，通過數據處理來選擇出一個最佳的初始參數值，作為基本SEM算法的輸入。獲得最佳初始參數值主要包括以下幾個步驟:

a)設數據集為D，給定初始模型M0，在執行SEM算法之前，隨機生成k個初始參數值，分別記為θ01，θ02，…，θ0k，計算當前θ0i(i=1，…，k)的似然函數期望:

L(θ|θ0i)=∑L∑Xi ln P(DL，XL|θ)P(XL|DL，θ0i)(2)

其中:DL和XL分別表示當前數據集和所有的變量。由式(2)得到θ0i(i=1，2，…，k)，分別對應的似然函數期望，記為L(θ|θ0i)(i=1，2，…，k)。

b)通過最大化當前期望似然函數值，選擇下一個估計:

θ1i=arg max E[P(D|θ)|D，θ0i，M0](3)

得到k個結果記為θ1i(i=1，2，…，k)。其中arg max()表示尋找具有最大評分的參量。設Q0i=E[P(D|θ)|D，θ0i，M0]，則式(3)可以寫成

θ1i=arg max Q0i(θ，θ0i)(4)

得到k個結果記為θ1i(i=1，2，…，k)。

c)在這些結果中選取一個最佳值記為θ0:

θ0=arg max Q(θ1i)(i=1，2，…，k)(5)

d)把由式(5)得到的θ0代入SEM算法的初始參數θ0，執行基本的SEM算法。

從θ0i選擇下一個估計θ1i的計算過程可以反復迭代很多次，在最不影響基本SEM算法時間性能的前提下，選定每個θ0i只進行一次迭代。經過這種初始化處理，可以從候選初始參數中選擇出使當前網絡評分得到最大值的參數作為基本算法的初始參數值，有效地減少了基本SEM算法執行過程中循環迭代的次數，提高了調控網絡構建的精度。

3 實驗與分析

3.1 實驗環境和數據

本文的實驗環境MATLAB 7.0，運行環境為操作系統Windows XP，CPU P43.0 GHz，內存512 MB，硬盤80 GB。在網絡構建實驗中采用了Leary等人[4]編寫的基于MATLAB的BNT(Bayesian network toolbox)Structure Learning Package。該軟件包是對Murphy等人開發的BNT工具箱在靜態網絡結構學習方面的擴充。

DNA微陣列使得生物學家能夠在基因組層次上研究任何種類細胞在任何時間、條件下的基因表達模式。本文使用的是由 Spellman等人[6]于 1998年提出的啤酒酵母 (saccharomyces cerevisiae)細胞周期微陣列表達譜數據集。這個數據集是采用三種不同的同步化方法，最終從酵母細胞中提取出由 800 個左右表達水平符合周期變化的基因組成的，包括了這些基因在不同實驗條件下或不同時間點上的具體表達水平。這些數據包含兩個短的時間片(兩個時間點對應實驗條件為CLN3、CLB2)和四個中等時間片(時間點分別為18、24、17和14個，對應實驗條件分別為ALPHA、CDC15、CDC28和ELU)。本文在基因網絡的學習中，采用了四個時間片。

目標網絡如圖2(a)所示，來源于KEGG[7]，是一個以CDC28 為中心的細胞周期通路。圖2(c)表示的是Kim等人[8]用動態貝葉斯模型在相同的數據下構建出來的基因調控網絡圖。

本文在實驗中采用改進的SEM算法，分別選取了k=10，15，20構建基因調控網絡，隨著k值的增大，學習的精度增加，但是同時數據處理步驟的計算量也增大。結果表明當k=15時，算法執行的時間性能基本不變，但是學習精度明顯提高，綜合效率最好。構建出來的基因調控網絡如圖2(b)所示。

在給出比較結果之前，首先介紹一下對實驗結果的表示方法和評價標準。網絡中的有向邊表示兩個基因之間存在調控關系，并且標志出調控與被調控的雙方。筆者使用三種符號來表示預測出的調控關系是否與現有生物學資料相符。有圓圈符號的，表示已經有文獻證實的調控關系;有三角符號的有向邊，表示調控關系存在，但是方向與預測結果相反，或者預測的調控關系至多只跳過了一個基因，即如果有X1→X2→X3 的調控關系存在，而算法給出的預測結果是X1→X3， 那么這條邊就會被標記為三角;如果有向邊的旁邊標注了叉號，則表示現有的文獻并沒有記載這樣的調控關系。值得注意的是，這樣的預測未必是錯誤的，而只能說現有的實驗并不能夠支持這樣的調控關系預測而已。

3.2 實驗結果與分析

為了定量地評價重構網絡的真實性，本文引入了文獻[8]中提到的 sensitivity(敏感度)和 specificity(特異度)兩個指標:sensitivity=正確估計的邊數 / 目標網絡中所有邊數， 是指網絡中實際存在的調控關系中，被算法正確發現的調控關系所占的比例;specificity=正確估計的邊數/所有估計的邊數，是指算法預測的調控關系總數中，正確預測所占的比例。換句話說，sensitivity 表示的是算法對于調控關系的敏感程度，specificity 衡量了算法預測的特異程度。在實驗中計算specificity和sensitivity時，目標網絡中總共的路徑數目是19。顯然，通過本文中的方法構建出來的網絡敏感度和特異度都比文獻[8]中的結果好，結果表明本文提出的改進的SEM算法構建出來的基因調控網絡更佳。本文的實驗結果與文獻[8]中結果的比較如表1所示。

表1 本文的實驗結果與文獻[8]中結果的比較

指標文獻[8]改進的SEM算法

正確估計

錯誤估計

方向相反或者跳過一個基因數

specificity

sensitivity43826.7%21.1%64540.0%31.6%

4 結束語

基因調控的理論和應用已經取得了許多重大的成果，但是隨著研究的不斷深入，所面臨的問題和挑戰也越來越多。SEM算法由于其在處理缺失值或部分可觀察數據上的優勢很適合從微陣列基因數據上構建基因調控網絡。本文提出了一種改進的SEM算法，通過數據的初始化處理，從候選初始參數中選取相對最好的初始值，然后執行SEM算法，從而更好地構建出基因調控網絡。通過實驗證實了本文方法更有效，在與文獻[8]的數據比較后可以看出結果網絡更接近于最佳網絡。但是，目前微陣列數據也存在一定的問題，如數據本身的噪聲以及在離散化過程中產生出來的噪聲，這些都可能會對網絡的結果造成影響;同時，本文的實驗中僅僅選取了節點很少的網絡，而對于具體不同的網絡選取什么樣的k值沒有找到一個很好的評判標準，在應用到不同規模的基因調控網絡構建方面有一定的局限性，需要將來進一步的完善。

參考文獻:

[1]AKUTSU T， KUHARA S， MIYANO S. Algorithms for identifying Boolean networks and related biological networks based on matrix multiplication and fingerprint function[J]. Journal of Computational Biology，2000，7(3-4):331-343.

[2]WAHDE M， HERTZ J.Coarse-grained reverse engineering of genetic regulatory networks[J]. Biosystems， 2000，55(1-3):129-136.

[3]FRIEDMAN N， LINIAL M， NACHMAN I， et al.Using Bayesian network to analyze expression data[J]. Journal of Computational Bio-logy， 2000，7(3-4):601-620.