方躍勝， 姚宏亮

(1.安徽水利水電職業(yè)技術(shù)學(xué)院，安徽合肥 231603;2.合肥工業(yè)大學(xué)計(jì)算機(jī)與信息學(xué)院，安徽合肥 230009)

0 引言

基因調(diào)控網(wǎng)絡(luò)(GRN，Gene Regulation Network)的研究目的是期望從系統(tǒng)的角度全面揭示基因組的功能和行為.基因調(diào)控網(wǎng)絡(luò)是由細(xì)胞中相互作用的基因片段以及其它起調(diào)控作用的物質(zhì)共同構(gòu)成的調(diào)控網(wǎng)絡(luò)，它表征了細(xì)胞中基因之間的調(diào)控關(guān)系，指導(dǎo)著基因到mRNA的表達(dá)過程[1].在基因調(diào)控網(wǎng)絡(luò)的建模方面主要的數(shù)學(xué)模型和方法有:線性模型、布爾網(wǎng)絡(luò)模型、聚類分析和貝葉斯網(wǎng)絡(luò)等[2].以上是系統(tǒng)生物學(xué)中常用的建模方法.但由于這些模型對數(shù)據(jù)和建模對象的要求很高，目前在實(shí)際應(yīng)用中還有許多問題有待解決.生物信息學(xué)研究的一個(gè)熱點(diǎn)是利用貝葉斯網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)，貝葉斯網(wǎng)絡(luò)被認(rèn)為是構(gòu)建基因調(diào)控網(wǎng)絡(luò)的最有前途的方法[3].

貝葉斯網(wǎng)絡(luò)亦稱信念網(wǎng)絡(luò)(Belief Network)，它是一種不確定性處理模型，用來模擬人類推理過程中因果關(guān)系的[4].根據(jù)處理的基因表達(dá)數(shù)據(jù)無序與有序的區(qū)別，貝葉斯網(wǎng)絡(luò)分為靜態(tài)和動態(tài)兩個(gè)分支.學(xué)習(xí)貝葉斯網(wǎng)就是要尋找一種能按某種測度最好地與給定實(shí)例數(shù)據(jù)集擬合的網(wǎng)絡(luò).貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)的關(guān)鍵在于尋找到一個(gè)好的搜索算法和一個(gè)對網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行有效評估的方法[5].

目前用于貝葉斯構(gòu)建基因調(diào)控網(wǎng)絡(luò)的貪心搜索GS、MWST和K2等網(wǎng)絡(luò)學(xué)習(xí)算法，常被用來處理完備數(shù)據(jù)下的結(jié)構(gòu)學(xué)習(xí).EM(Expectation Maximization)算法通過對數(shù)據(jù)進(jìn)行充分統(tǒng)計(jì)來得到后驗(yàn)概率分布的平均估計(jì)，然后通過最大似然得到和數(shù)據(jù)擬合的最好的網(wǎng)絡(luò).另外，EM算法對用于學(xué)習(xí)的數(shù)據(jù)依賴小，數(shù)據(jù)可以不完備.標(biāo)準(zhǔn)的EM算法是用于丟失數(shù)據(jù)的參數(shù)估計(jì)的，它是對固定網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行參數(shù)學(xué)習(xí).Friedman[6]等人提出了一種擴(kuò)展的EM算法可以進(jìn)行結(jié)構(gòu)學(xué)習(xí)，稱為MS－EM算法(Model Selection EM algorithm)[7].

圖1 改進(jìn)的MS－EM算法實(shí)現(xiàn)結(jié)構(gòu)

1 改進(jìn)的MS－EM算法

本文采用改進(jìn)的MS－EM算法進(jìn)行結(jié)構(gòu)學(xué)習(xí).因?yàn)閺挠?jì)算上說，最大化一個(gè)固定模型的參數(shù)比搜索一個(gè)更好的模型代價(jià)小的多，這個(gè)參數(shù)的最大化可以使用參數(shù)EM算法得到，參數(shù)EM算法運(yùn)行指定的步數(shù)，或者運(yùn)行到收斂為止，然后才考慮改變模型的選擇來獲取評分的改進(jìn).這種算法在某些情況下還可以避免局部極大化.改進(jìn)的算法描述如下:

圖2 系統(tǒng)開發(fā)流程圖

圖3 實(shí)驗(yàn)平臺的進(jìn)入界面

隨機(jī)的選擇 M0和 θ0，loop n=0，1，….直至算法收斂，loop l=0，1，….直至算法收斂或者 l=lmax，使得，尋找一個(gè)模型 Mn+1，使得最大化讓參數(shù)).改進(jìn)的算法實(shí)現(xiàn)結(jié)構(gòu)如圖1所示.

搜索引擎模塊負(fù)責(zé)選擇等待評價(jià)的模型，對于每個(gè)待選模型，搜索引擎調(diào)用評分函數(shù)模塊，它執(zhí)行某種評分函數(shù)(例如:MDL或者Bayesian評分函數(shù)，本次設(shè)計(jì)使用MDL評分函數(shù))，為了評價(jià)得分，評分模塊需要調(diào)用相應(yīng)的統(tǒng)計(jì)值，統(tǒng)計(jì)模塊通過計(jì)算樣本數(shù)據(jù)的期望來“回答”評分函數(shù).這個(gè)算法根據(jù)統(tǒng)計(jì)模塊和推理算法的不同，會有細(xì)微的變化.

2 算法的設(shè)計(jì)與實(shí)現(xiàn)

2.1 基因數(shù)據(jù)來源及預(yù)處理

目前研究基因調(diào)控網(wǎng)絡(luò)的主要手段之一，是使用概率圖模型中的貝葉斯網(wǎng)絡(luò)[8].貝葉斯網(wǎng)絡(luò)具有豐富的隨機(jī)特性，所以特別適合于處理微陣列數(shù)據(jù)這樣由于實(shí)驗(yàn)條件限制而含有大量噪聲的數(shù)據(jù).本文中使用Spellman[9]等人于1998年提出的啤酒酵母(Saccharomyces Cerevisiae)細(xì)胞周期微陣列表達(dá)譜數(shù)據(jù)集.這些基因表達(dá)數(shù)據(jù)集描述了酵母菌在細(xì)胞周期過程中的基因表達(dá)特性.算法中要求數(shù)據(jù)是離散的，所以在編寫算法中將連續(xù)的數(shù)據(jù)進(jìn)行了離散化處理.

2.2 系統(tǒng)的總體要求

開發(fā)環(huán)境:用VC作界面，調(diào)用MATLAB下的M腳本文件實(shí)現(xiàn)功能.

功能要求:(1)采用EM算法從樣本學(xué)習(xí)貝葉斯網(wǎng)絡(luò)的結(jié)構(gòu)，輸出學(xué)習(xí)所用時(shí)間和得分;(2)實(shí)現(xiàn)可視化界面，進(jìn)行可視化操作，輸出最終得到的基因調(diào)控網(wǎng)絡(luò)圖.

2.3 系統(tǒng)的實(shí)現(xiàn)流程

本系統(tǒng)是在Windows環(huán)境下，采用VC6.0作界面，以調(diào)用 MATLAB引擎的形式執(zhí)行 MATLAB7.0下編寫的M腳本文件和函數(shù)、參數(shù)傳入和傳出，從而把VC的可視化與MATLAB強(qiáng)大的科學(xué)計(jì)算這兩個(gè)優(yōu)點(diǎn)相結(jié)合.系統(tǒng)開發(fā)流程如圖2所示.(至于腳本文件與函數(shù)的編寫，開發(fā)環(huán)境設(shè)置及VC與MATLAB的參數(shù)傳遞，這里不再贅述.)

根據(jù)M腳本文件，設(shè)計(jì)VC界面，實(shí)現(xiàn)MATLAB下的抽象算法的可視化.實(shí)驗(yàn)中選取了兩組數(shù)據(jù)分別進(jìn)行試驗(yàn)，并分別與參考文獻(xiàn)中的正確網(wǎng)絡(luò)圖進(jìn)行比較，以察看算法的正確率.該系統(tǒng)是靜態(tài)概率模型構(gòu)建基因調(diào)控網(wǎng)絡(luò).實(shí)驗(yàn)平臺的VC界面設(shè)計(jì)及運(yùn)行效果如圖3所示.

點(diǎn)擊“進(jìn)入”:進(jìn)入實(shí)驗(yàn)平臺的MS－EM算法學(xué)習(xí)界面，如圖4、5所示.

圖4 算法學(xué)習(xí)界面

圖5 打開文件獲取路徑名

圖6 輸入最大迭代次數(shù)執(zhí)行EM算法

圖7 結(jié)果分析頁面

(1)第一組數(shù)據(jù)

打開數(shù)據(jù)組1，構(gòu)建酵母菌的細(xì)胞周期調(diào)控網(wǎng)絡(luò)，使用的數(shù)據(jù)選擇Futcher[10]所列的主要轉(zhuǎn)錄因子:MBP1，SWI4，SWI6，F(xiàn)KS1，ACE2，CDC28，CLN3，CLB2，SIC，CLN2，HHT1，MCM1，F(xiàn)KH1，NDD1和SWI5(FKH2除外，因?yàn)樵谒玫臄?shù)據(jù)庫中，F(xiàn)KH2基因的缺失值的比例過大).其中MBP1，SWI4和SWI6主要結(jié)合在G1晚期基因的啟動子，MCM1，F(xiàn)KH2和NDD1結(jié)合于G1/M期基因，SWI5和ACE2結(jié)合于M/G1期基因，另外FKH1在G1、S和G1/M期都有結(jié)合.

下面是EM算法學(xué)習(xí)運(yùn)行界面:可輸入最大迭代次數(shù)，然后開始學(xué)習(xí).點(diǎn)擊“開始學(xué)習(xí)”按紐，則調(diào)用MATLAB引擎，執(zhí)行em腳本文件.如圖6所示.

算法運(yùn)行耗時(shí)，得分:為了評價(jià)重構(gòu)的網(wǎng)絡(luò)，并且選擇與真實(shí)情況最接近的網(wǎng)絡(luò)結(jié)構(gòu)，需要對網(wǎng)絡(luò)進(jìn)行評分.在同樣的標(biāo)準(zhǔn)下，得分越高的網(wǎng)絡(luò)就越接近真實(shí)情況.如圖7所示.

顯示最終網(wǎng)絡(luò)圖:針對基因調(diào)控網(wǎng)絡(luò)建模問題來說，貝葉斯網(wǎng)絡(luò)的節(jié)點(diǎn)表示的是不同的基因，而節(jié)點(diǎn)之間的邊表示兩個(gè)基因之間是否有基于概率的調(diào)控關(guān)系，邊的方向表示的是調(diào)控方和被調(diào)控方.實(shí)驗(yàn)中共使用了15個(gè)基因作為節(jié)點(diǎn)，基因之間調(diào)控關(guān)系如圖8所示.

圖8 基因調(diào)控網(wǎng)絡(luò)圖

(2)第二組數(shù)據(jù)

數(shù)據(jù)來源:我們使用的仍然是從先前的對酵母細(xì)胞周期的基因調(diào)控(Spellman et al.，1998)的研究得到的公開的數(shù)據(jù).共18個(gè)基因，包括組蛋白基因 HHT1，HHF1，HHF2，HTA1，HTA2，HTB1 和HTB2，還包括在細(xì)胞周期中比組蛋白基因執(zhí)行不同功能(如DNA初始復(fù)制，DNA損傷修復(fù)等)的基因.這些酵母菌數(shù)據(jù)從http://genome－www.stanford.edu/cellcycle/data/rawdata/individual.html 下載，并根據(jù) http://www.yeastgenome.org/SGD 數(shù)據(jù)庫將系統(tǒng)名稱轉(zhuǎn)為標(biāo)準(zhǔn)名稱.

為了進(jìn)一步評估算法的性能，我們選擇第二組數(shù)據(jù)開始學(xué)習(xí)，同樣迭代30次.如圖9、10所示.

算法耗時(shí)，得分:

圖9 算法執(zhí)行

圖10 算法結(jié)果

最終基因調(diào)控網(wǎng)絡(luò)圖:基因之間調(diào)控關(guān)系如圖11所示.

圖11 基因調(diào)控網(wǎng)絡(luò)圖

3 算法的實(shí)驗(yàn)與分析

3.1 結(jié)果比較

(1)實(shí)驗(yàn)一

參考 Simon I，Barnett J，Hannett N，et a1.Serial regulation of transcriptional regulators in the yeast cell cycle.Cell，2001;106(6):697－708 中正確的調(diào)控關(guān)系，繪制基因調(diào)控網(wǎng)絡(luò)圖.由實(shí)驗(yàn)結(jié)果繪制的基因調(diào)控網(wǎng)絡(luò)圖如圖12所示.

圖12 由數(shù)據(jù)組一得到的基因調(diào)控網(wǎng)絡(luò)

圖13 基因調(diào)控關(guān)系參考圖

圖14 由數(shù)據(jù)組二得到的基因調(diào)控網(wǎng)絡(luò)

實(shí)線表示已經(jīng)被驗(yàn)證存在的邊結(jié)果中有一條，占總邊數(shù)的6%;點(diǎn)劃線表示錯誤的邊(與正確邊方向相反的邊)，由于偶然或數(shù)據(jù)噪聲造成，結(jié)果中有兩條，占總邊數(shù)的12%;長劃線表示的邊為未被現(xiàn)有文獻(xiàn)證明的，其中有些可能是正確的.

(2)實(shí)驗(yàn)二

正確的基因調(diào)控關(guān)系參考如圖13所示，基因之間的調(diào)控關(guān)系是通過PathwayAssist軟件獲得的，作為本次實(shí)驗(yàn)的參考.

通過本次實(shí)驗(yàn)繪制的基因調(diào)控網(wǎng)絡(luò)圖如圖14所示:

圖中實(shí)線表示可由PathwayAssist得出的邊，即已被證明是正確的邊，共五條，占總邊數(shù)的19%;點(diǎn)劃線表示錯誤的邊，即與已證明的正確邊反向的邊，占12%;長劃線表示為現(xiàn)有文獻(xiàn)還不能證明的邊.

3.2 結(jié)果分析

學(xué)習(xí)結(jié)果中文獻(xiàn)尚未證明的邊較多，正確的邊相對較少，分析其原因，一是現(xiàn)有技術(shù)正確率都不太高，且尚未被文獻(xiàn)證明的邊有可能是正確的，因?yàn)槟壳瓣P(guān)于基因調(diào)控關(guān)系的已有知識還是比較匱乏的.貝葉斯網(wǎng)絡(luò)方法不僅可以構(gòu)建已經(jīng)存在的基因關(guān)系，而且可以發(fā)現(xiàn)一些基因間的新關(guān)系.二是所使用的基因數(shù)據(jù)有限.若要使得學(xué)習(xí)結(jié)果與真實(shí)網(wǎng)絡(luò)更加擬和，可提高所使用樣本的數(shù)量，樣本數(shù)量越多，學(xué)習(xí)結(jié)果就與真實(shí)越趨近.

4 結(jié)束語

本文借助于VC6.0和MATLAB7.0通過分析、設(shè)計(jì)、到編程完整實(shí)現(xiàn)了基因調(diào)控網(wǎng)絡(luò)EM算法的學(xué)習(xí)，并設(shè)計(jì)了友好的用戶使用界面，用于顯示每一步的結(jié)果并與正確結(jié)果進(jìn)行比較.實(shí)驗(yàn)結(jié)果表明改進(jìn)后的算法進(jìn)一步降低了時(shí)間性能，提高了構(gòu)建調(diào)控網(wǎng)絡(luò)的精度.但仍存在不足之處，需進(jìn)一步改進(jìn)EM算法，這也是今后應(yīng)該解決的問題.另外，考慮到EM算法耗時(shí)太長，本文只是在一個(gè)很小規(guī)模的模型(15個(gè)節(jié)點(diǎn))上進(jìn)行了試驗(yàn)，今后應(yīng)該在該改進(jìn)MS－EM算法的時(shí)間性能的基礎(chǔ)上，對大規(guī)模的網(wǎng)絡(luò)進(jìn)行學(xué)習(xí)，從而把EM算法應(yīng)用于實(shí)際生活當(dāng)中.

[1]DE Jong H.Modeling and Simulation of Genetic Regulation System:a Literature Review[J].Journal of Computational Biology，2002，9(1):67－103.

[2]Friedman N，Linial M，Nachman I，et a1.Using Bayesian Network to Analyze Expression Data[J].Journal of Computational Biology，2000，7(1):6O1－62O.

[3]Heckerman D，Geiger D，Chickering D M.Learning Bayesian Networks:the Combination of Knowledge and Statistical Data[J].Machine Learning，1995，20:197－ 243.

[4]D.Heckerman，Atutorialon Learning with Bayesian Networks[J].Microsoft Research Tech Report，MSR－ TR－ 95－ 06.1996，300－302.

[5]強(qiáng)波，王正志.基于動態(tài)貝葉斯網(wǎng)構(gòu)建基因調(diào)控網(wǎng)絡(luò)[J].生物醫(yī)學(xué)工程研究，2008，27(3):145－149.

[6]Nir Friedman，“Learing Belief Networks in the Presence of Missing Values Hidden Variables”[J].Computer Science Division U-niversity of alifornia，Berkeley CA94720，1998:139－147.

[7]Liu.C and Rubin，D.B(1994).The ECME Algorithm:A Simple Extension of MS－EM and ECM with Faster Monotone Convergence[J].Biometrika 81.

[8]Lauritzen S L.The EM Algorithm for Graphical Association Models with Missing Data[J].Comput.Stat Data Anal，1995，(19):191－201.

[9]Spellman P T，Sherlock G，Zhang M Q，et a1.Comprehensive Dentification of Cell Cycle－Regulated Genes of the Yeast Saccharo Myce Scerevisiae by Microarray Hybridization[J].Molecular Biology of the Cell，1998，9:3273－3297.

[10]Futcher B.Transcripti0nal Regulation Networks and Yeast Cell Cycle[J].Current Opinion in Cell Biology，2002，14:676－683.