緱葵香，宮秀軍，湯 莉

(1. 天津大學(xué)理學(xué)院，天津 300072；2. 天津大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，天津 300072)

時(shí)間序列(time series)微陣列數(shù)據(jù)，可以反映一組基因在生命活動(dòng)周期的時(shí)間序列條件下的表達(dá)水平的變化，其表達(dá)水平變化的時(shí)間延遲關(guān)系可反映基因調(diào)控關(guān)系，因而近年來(lái)被廣泛地用于基因之間轉(zhuǎn)錄調(diào)控關(guān)系尋找和基因調(diào)控網(wǎng)絡(luò)的構(gòu)建．構(gòu)建和分析基因調(diào)控網(wǎng)絡(luò)，從分子水平認(rèn)識(shí)細(xì)胞內(nèi)的生理活動(dòng)和功能，了解通路中的相互作用，以及如何使生物體產(chǎn)生變化．基于微陣列數(shù)據(jù)構(gòu)建基因調(diào)控網(wǎng)絡(luò)，逐漸成為生物信息學(xué)的研究熱點(diǎn)．

很多數(shù)學(xué)模型如布爾網(wǎng)絡(luò)模型、微分方程模型、靜態(tài)貝葉斯網(wǎng)絡(luò)模型和動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)模型[1]等被提出用于構(gòu)建基因調(diào)控網(wǎng)絡(luò)．相對(duì)來(lái)講，布爾網(wǎng)絡(luò)模型是一種簡(jiǎn)單的離散模型，對(duì)系統(tǒng)的模擬是定性的、較為粗糙的．微分方程所給出的模型是確定性的、連續(xù)的，然而由于細(xì)胞中的分子受到熱力學(xué)波動(dòng)和噪聲過(guò)程的支配，基因表達(dá)是一個(gè)隨機(jī)過(guò)程，微分方程模型不能很好地描述基因表達(dá)的隨機(jī)性．靜態(tài)貝葉斯網(wǎng)絡(luò)是一個(gè)有向無(wú)環(huán)的概率圖形，它用概率來(lái)表示調(diào)控關(guān)系，更符合生物學(xué)實(shí)際．靜態(tài)貝葉斯網(wǎng)絡(luò)模型可以成功地從非時(shí)序數(shù)據(jù)中推斷出隨機(jī)變量之間的因果關(guān)系，但是不能對(duì)動(dòng)態(tài)時(shí)間數(shù)據(jù)進(jìn)行建模，因而不能成功地從時(shí)間序列的微陣列數(shù)據(jù)中構(gòu)建基因調(diào)控網(wǎng)絡(luò)．動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)(dynamic Bayesian network，DBN)是在靜態(tài)貝葉斯網(wǎng)絡(luò)的基礎(chǔ)上，引入了時(shí)間維而形成的動(dòng)態(tài)網(wǎng)絡(luò)，動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)能夠?qū)W習(xí)隨機(jī)變量間的概率依存關(guān)系及其隨時(shí)間變化的規(guī)律，并以圖的方式直觀地反映這種關(guān)系，由于微陣列數(shù)據(jù)的時(shí)間特性，動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)可以更精確地描述基因調(diào)控網(wǎng)絡(luò)，所以基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)從時(shí)間序列微陣列數(shù)據(jù)中構(gòu)建基因調(diào)控網(wǎng)絡(luò)成為當(dāng)前的一個(gè)研究熱點(diǎn)．目前，已有許多基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)的相關(guān)算法[2-5]，但這些算法都需要將連續(xù)的基因表達(dá)數(shù)據(jù)離散化，沒(méi)有充分利用數(shù)據(jù)所包含的信息．

筆者提出了一個(gè)基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)的算法，該算法基于時(shí)序互信息來(lái)學(xué)習(xí)動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)，在計(jì)算基因間的時(shí)序互信息時(shí)，該算法考慮了時(shí)間序列微陣列數(shù)據(jù)的時(shí)間特性，并利用協(xié)方差矩陣計(jì)算互信息，沒(méi)有將基因表達(dá)數(shù)據(jù)離散化，與基因表達(dá)數(shù)據(jù)的連續(xù)性相符合，減少了信息的損失．

1 動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)及相關(guān)研究

1.1 動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)

動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)是貝葉斯網(wǎng)絡(luò)的擴(kuò)展，它包含隨時(shí)間變化的隱含的狀態(tài)，使隱含的隨機(jī)進(jìn)程的熵的比率最大化，它是有向有環(huán)圖．設(shè)一個(gè)時(shí)間序列X ∈ Rs×n，s為時(shí)間點(diǎn)的數(shù)目，n為基因的數(shù)目．(t)表示基因i關(guān)于時(shí)間t的表達(dá)水平的變量，在時(shí)刻t，n個(gè)基因表達(dá)水平的向量為 X (t) = (X1(t) , … ,Xn(t ))T，基因 i在 s個(gè)時(shí)間點(diǎn)上表達(dá)水平的向量為 Xi=(Xi(1),… , Xi(s ))T．

動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)假定：①網(wǎng)絡(luò)中的有向弧隨著時(shí)間的變化而展開(kāi)，即假定動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)是一階馬爾科夫的，基因 i在 t時(shí)刻的表達(dá)水平只依賴于它的父節(jié)點(diǎn)在t-1時(shí)刻的表達(dá)水平，而與 t-1時(shí)刻前的表達(dá)狀態(tài)完全獨(dú)立；②基因之間的動(dòng)態(tài)因果關(guān)系在所有的時(shí)間點(diǎn)上都沒(méi)有變化，即隨機(jī)變量的集合以及相應(yīng)的概率轉(zhuǎn)移函數(shù)，對(duì)每個(gè)時(shí)間點(diǎn)來(lái)說(shuō)都是相同的．根據(jù)以上假設(shè)，網(wǎng)絡(luò)的條件概率分布可表示為

式中 P ai(t ? 1 )為基因 Xi(t)的父節(jié)點(diǎn)．一個(gè)典型的動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)如圖1所示．

圖1 一個(gè)簡(jiǎn)單的動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)Fig.1 A simple dynamic Bayesian network

從數(shù)據(jù)中學(xué)習(xí)動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)分為結(jié)構(gòu)學(xué)習(xí)和參數(shù)學(xué)習(xí)兩個(gè)過(guò)程．給定一個(gè)時(shí)間序列的數(shù)據(jù)集D = ( X ( 1),… ,X (n))，從 D中學(xué)習(xí)動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)就是找到一個(gè)與 D最佳匹配的結(jié)構(gòu) G = ( V,E)，其中V = { X1,… ,Xn}為頂點(diǎn)的集合，Xi表示基因 i，E = { (Xi(t ? 1 ),Xj(t) ) |Xi(t ? 1 ) → Xj(t)}為 一 組 有 向 弧的集合，Xi(t ? 1 ) → Xj(t )表示 t-1時(shí)刻的基因 i對(duì) t時(shí)刻的基因 j有調(diào)控關(guān)系，即基因 i在 t-1時(shí)刻的表達(dá)水平影響了基因j在t時(shí)刻的表達(dá)水平．研究人員對(duì)用動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)從時(shí)間序列基因數(shù)據(jù)中構(gòu)建基因調(diào)控網(wǎng)絡(luò)的模型做了大量的研究工作．

1.2 相關(guān)研究

Wu等[2]討論了利用 BDE(Bayesian dirichlet equivalence)打分函數(shù)來(lái)構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的貪婪的爬山搜索算法．給定一個(gè)時(shí)間序列的數(shù)據(jù)集D = ( X ( 1),… ,X (n))，從 D中構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)就是發(fā)現(xiàn)一個(gè)與D最匹配的模型 M =(G,Θ)，G是網(wǎng)絡(luò)結(jié)構(gòu)，Θ是條件概的集合．BDE打分函數(shù)的方法對(duì)父節(jié)點(diǎn)的數(shù)目比較敏感，容易陷入局部最優(yōu)，計(jì)算復(fù)雜性指數(shù)級(jí) O (μ2N)，N為網(wǎng)絡(luò)中節(jié)點(diǎn)的數(shù)目，μ為指定的父節(jié)點(diǎn)的上界．計(jì)算復(fù)雜度高，而且這種方法需要將基因表達(dá)數(shù)據(jù)離散化．

Wu等[2]還討論了利用 MCMC(Markov chain Monte Carlo)方法構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的 MH(metropolis hastings)算法．如果一個(gè)過(guò)程的“將來(lái)”僅依賴“現(xiàn)在”而不依賴“過(guò)去”，則此過(guò)程具有馬爾可夫性，或稱此過(guò)程為馬爾可夫過(guò)程，X(t)=f(X(t-1))．時(shí)間和狀態(tài)都離散的馬爾科夫過(guò)程稱為馬爾科夫鏈．設(shè) S是一個(gè)可數(shù)集，{ X (t),t = 1 ,2,… ,s } 是一組隨機(jī)變量的集合，如果滿足 P (X (t)| X (t ? 1),… ,X (1))=P(X(t) |X (t ? 1 ))，則稱{ X (t),t = 1 ,2,… ,s }為一個(gè)馬爾科夫鏈．對(duì)于一個(gè)馬爾科夫鏈，如果 P (X(t)|X(t ?1))與 t的取值無(wú)關(guān)，則稱這個(gè)馬爾科夫鏈為齊氏馬爾科夫鏈．動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的轉(zhuǎn)移網(wǎng)絡(luò)就是一個(gè)齊氏的馬爾科夫鏈．

這種方法也需要將基因表達(dá)數(shù)據(jù)離散化，運(yùn)算時(shí)間復(fù)雜度高，即使在數(shù)據(jù)充分的情況下，學(xué)習(xí)的結(jié)果也不是很理想．

Zhu等[3]根據(jù)生物學(xué)的細(xì)胞活化信號(hào)路徑等先驗(yàn)知識(shí)，在基于貝葉斯網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)時(shí)，縮小了搜索空間，降低了時(shí)間復(fù)雜度，擴(kuò)大了構(gòu)建基因調(diào)控網(wǎng)絡(luò)的規(guī)模

Wang等[4]針對(duì)趨勢(shì)相關(guān)(兩基因在其表達(dá)水平隨時(shí)間上升與下降的變化趨勢(shì)上相關(guān))關(guān)系在重建基因調(diào)控網(wǎng)絡(luò)中十分重要卻尚未被挖掘利用的問(wèn)題，提出了幾何模式動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)(Gp-DBN)方法．Gp-DBN將每個(gè)基因的表達(dá)數(shù)據(jù)轉(zhuǎn)換為一個(gè)幾何模式，依據(jù)幾何模式確定潛在的調(diào)控子和調(diào)控時(shí)滯，并通過(guò)推理這些幾何模式之間的相關(guān)關(guān)系來(lái)發(fā)現(xiàn)基因間的調(diào)控關(guān)系．該方法解決了挖掘具有趨勢(shì)相關(guān)的基因調(diào)控關(guān)系的問(wèn)題，能夠很大程度地提高重建的基因調(diào)控網(wǎng)絡(luò)的性能．

Liu等[5]指出在基因調(diào)控網(wǎng)絡(luò)中不僅單個(gè)蛋白質(zhì)對(duì)基因具有調(diào)控作用，而且一些蛋白質(zhì)的組合對(duì)基因也具有調(diào)控作用，Liu把這些蛋白質(zhì)的組合看作基因調(diào)控網(wǎng)絡(luò)中的隱藏變量．Liu提出了一個(gè)Semi-fixed模式將基因調(diào)控網(wǎng)絡(luò)表示為具有隱藏變量的貝葉斯網(wǎng)絡(luò)，并用EM算法構(gòu)造Semi-fixed模式的基因調(diào)控網(wǎng)絡(luò)．

Dojer等[6]提出一種從基因擾動(dòng)實(shí)驗(yàn)的數(shù)據(jù)中構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的算法．在構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)前，作者先用 Husmeier離散化方法對(duì)基因表達(dá)數(shù)據(jù)做了預(yù)處理，在構(gòu)建貝葉斯網(wǎng)絡(luò)時(shí)，作者使用 BDE打分函數(shù)，并用從擾動(dòng)實(shí)驗(yàn)的數(shù)據(jù)集中獨(dú)立計(jì)算每個(gè)節(jié)點(diǎn)(基因)的父節(jié)點(diǎn)集的方法代替啟發(fā)式搜索的方法，減少了學(xué)習(xí)時(shí)間，減少了啟發(fā)式搜索的方法帶來(lái)的時(shí)間復(fù)雜性．

以上的基于動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)的算法均需對(duì)基因表達(dá)數(shù)據(jù)離散化，而時(shí)序的微陣列數(shù)據(jù)是連續(xù)的，從而不能充分利用數(shù)據(jù)信息．這里筆者提出了TSMI-GRN(time series mutual informationgene regulation network)算法，該算法基于時(shí)序互信息構(gòu)建基因調(diào)控網(wǎng)絡(luò)，所構(gòu)建的網(wǎng)絡(luò)為動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)，TSMI-GRN算法不需要對(duì)基因表達(dá)數(shù)據(jù)離散化，減少了數(shù)據(jù)信息的損失．

2 TSMI-GRN算法

2.1 時(shí)序互信息

設(shè)D為一個(gè)時(shí)間序列基因表達(dá)數(shù)據(jù)集，包含n個(gè)基因{1,2,…,n}，s個(gè)時(shí)間點(diǎn)，Xi(t)表示基因 i關(guān)于時(shí)間 t的表達(dá)水平的變量，這里要在 D上用 TSMIGRN算法構(gòu)建這 n個(gè)基因的調(diào)控網(wǎng)絡(luò)．設(shè)為一多維平穩(wěn)時(shí)間序列 ．

要在數(shù)據(jù) D上構(gòu)建動(dòng)態(tài)貝葉斯網(wǎng)絡(luò) G = ( V,E)，V = { 1,2,… ,n }，為 n個(gè)基因的集合，E= {(Xi(t ?1),Xj(t) ) |Xi(t ? 1) → Xj(t),i = 1,… ,n, j = 1,… ,n}．在 構(gòu) 建網(wǎng)絡(luò)時(shí)，需要研究 Xi(t ? 1 )與 Xj(t)之間是否存在關(guān)系，即考慮t-1時(shí)刻基因i的表達(dá)水平對(duì)t時(shí)刻基因j的表達(dá)水平的影響，根據(jù)動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)的假設(shè)，t時(shí)刻的狀態(tài)只與 t-1時(shí)刻的狀態(tài)有關(guān)，這里考慮在 t-1時(shí)刻除去i與j兩個(gè)基因外所有基因的表達(dá)水平的條件下，來(lái)研究 Xj(t)和 Xi(t ? 1 )之間的條件獨(dú)立性，這里給出t時(shí)刻基因j與t-1時(shí)刻基因i的時(shí)序互信息的定義．

定義1 設(shè) Yij( t ? 1)= { Xk(t ? 1 ),k ≠ i ,j} ，Xj(t)和Xi(t ? 1 )在 Yij(t ? 1 )條件下的條件互信息為

對(duì)于給定的閾值ε＞0，如果I(Xj( t),Xi(t?1)|Yij(t ?1))＜ε，稱基因j與i相互獨(dú)立，否則稱j與i是相關(guān)的．

在式(2)中，令

則式(2)可簡(jiǎn)化為

這里假定隨機(jī)變量 X1,… ,Xn的聯(lián)合分布是n維高斯分布，線性熵 Hl( X1, X2,… ,Xn)和熵H(X1,X2,… ,Xn)在理論上是等價(jià)的．線性熵

式中：Σ為 n維隨機(jī)變量 ( X1, X2,… ,Xn)的協(xié)方差矩陣．

用 Σ、Σi、 Σj和 Σij分 別 表 示 前 面 所 構(gòu) 造 變 量Z(t) 、 Zi(t) 、 Zj(t) 和 Zij(t)的協(xié)方差矩陣，根據(jù)式(3)和式(4)，給定 Yij(t ? 1 )的條件下，Xj(t)和 Xi(t ? 1 )之間的條件互信息為

2.2 TSMI-GRN算法

定義 2定義了基于時(shí)間序列互信息所構(gòu)建的網(wǎng)絡(luò)結(jié)構(gòu)．

定義2 設(shè)X(t) = { X (t) , X (t) , … ,X (t) }T(t ∈Z)

1 2n為一多維平穩(wěn)時(shí)間序列，圖 G = ( V,E)，對(duì)應(yīng)的頂點(diǎn)集V ={1,2,… ,n}，對(duì) 于 給 定 ε＞0，如 果 Il(Xj(t),Xi(t ? 1 )|Yij(t ? 1 ))＞ ε ，則(Xi(t ? 1 ) → Xj(t) ) ∈ E ，如 果Il(Xj( t),Xi( t ? 1 )|Yij(t ? 1 ))＜ ε ，則(Xi(t ? 1 ) → Xj(t ))?E．

根據(jù)時(shí)序互信息和圖結(jié)構(gòu)的定義，提出了基于TSML-GRN算法從單個(gè)基因表達(dá)數(shù)據(jù)集中構(gòu)建基因調(diào)控網(wǎng)絡(luò)．TSML-GRN算法描述如下．

步驟 1 對(duì)于每個(gè)基因 j，計(jì)算 Il(Xj( t),Xj(t ?1)|Yj(t ? 1 ), j = 1 ,… ,n ).如果Il(Xj( t),Xj( t ? 1 )|Yj(t ? 1))＞ε，則(Xj(t ? 1 ) → Xj(t) ) ∈ E ，表明t-1時(shí)刻的基因j的表達(dá)水平對(duì)t時(shí)刻的基因j的表達(dá)水平具有調(diào)控作用．

步驟2 對(duì)于給定的基因i和j，i≠j，計(jì)算Il(Xj( t),Xi( t ? 1 )|Yij(t ? 1 ))，其中 j = 1 ,2,…, n ,i = 1,2,… ,n.如 果 Il(Xj( t),Xi( t? 1 )|Yij(t? 1))＞ ε ,則(Xi(t? 1 )→ Xj(t ))∈E，表明 t-1時(shí)刻基因 i的表達(dá)水平對(duì) t時(shí)刻基因 j的表達(dá)水平具有調(diào)控作用．計(jì)算 Il(Xi( t),Xj(t ?1)|Yij(t ? 1 ))，如 果 Il(Xi( t),Xj(t ? 1 )|Yij(t ?1))＞ε ，則(Xj(t ? 1 ) → Xi(t) ) ∈ E ，表明t-1時(shí)刻基因j的表達(dá)水平對(duì)t時(shí)刻基因i的表達(dá)水平具有調(diào)控作用．

步驟3 根據(jù)步驟1與步驟2，構(gòu)建了一個(gè)初始的動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu) G = ( V,E)，由于基因調(diào)控網(wǎng)絡(luò)具有稀疏性[7]，再利用一個(gè)打分函數(shù)對(duì)網(wǎng)絡(luò)結(jié)構(gòu)G進(jìn)行優(yōu)化．設(shè)定每個(gè)節(jié)點(diǎn)的父節(jié)點(diǎn)的上界μ=3，對(duì)于G中父節(jié)點(diǎn)數(shù)目大于3的節(jié)點(diǎn)，選取所有小于等于3的父節(jié)點(diǎn)的非空子集，選取使打分函數(shù)I(X , subsetpa(X ) /max{ H (X ) ,H(subsetpa(X ) )}[8]為最大分?jǐn)?shù)的子集作為該節(jié)點(diǎn)父節(jié)點(diǎn)集，得到最終的網(wǎng)絡(luò)結(jié)構(gòu)．

在上述的TSML-GRN算法中，步驟1計(jì)算的是t-1時(shí)刻基因j與t時(shí)刻基因j的網(wǎng)絡(luò)關(guān)系，步驟2計(jì)算的是t-1時(shí)刻基因i與t時(shí)刻基因j，以及t-1時(shí)刻基因 j與 t時(shí)刻基因 i的網(wǎng)絡(luò)關(guān)系．經(jīng)過(guò)步驟1和步驟2的計(jì)算，基于時(shí)序互信息構(gòu)建了一個(gè)初始的基因調(diào)控網(wǎng)絡(luò)，但這個(gè)網(wǎng)絡(luò)中每個(gè)節(jié)點(diǎn)的父節(jié)點(diǎn)的數(shù)目較大，而基因調(diào)控網(wǎng)絡(luò)具有稀疏性，所以步驟3引入了一個(gè)打分函數(shù)來(lái)優(yōu)化步驟1與步驟 2構(gòu)建的初始網(wǎng)絡(luò)結(jié)構(gòu)，得到最終的網(wǎng)絡(luò)結(jié)構(gòu)．下面在真實(shí)的基因表達(dá)數(shù)據(jù)集上測(cè)試 TSMI-GRN算法，根據(jù)實(shí)驗(yàn)結(jié)果，對(duì)算法的有效性進(jìn)行分析．

3 實(shí)驗(yàn)分析

在酵母菌周期細(xì)胞微陣列數(shù)據(jù)集 Cdc15上測(cè)試TSML-GRN算法，該數(shù)據(jù)集從斯坦福微陣列數(shù)據(jù)庫(kù)(http：//smd/stanford.edu/)中獲得．這里選擇了數(shù)據(jù)集Cdc15分別包含 14個(gè)基因和 25個(gè)基因的兩個(gè)子集進(jìn)行測(cè)試．測(cè)試結(jié)果與基因組數(shù)據(jù)庫(kù) KEGG中描述的通過(guò)實(shí)驗(yàn)得到的酵母菌周期細(xì)胞中一些基因的調(diào)控網(wǎng)絡(luò)(http：//www.genome.jp/dbget-bin/show_ pathway?dce04111)相比較，來(lái)分析實(shí)驗(yàn)結(jié)果．在 TSMIGRN算法中，假定了基因表達(dá)數(shù)據(jù)滿足正態(tài)分布．利用 SAS軟件中的 Kolmogorov-Smirnov方法對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了是否為高斯分布的檢驗(yàn)，檢驗(yàn)結(jié)果表明基因表達(dá)水平數(shù)據(jù)大致服從對(duì)數(shù)正態(tài)分布，證明這一假設(shè)是合理的．在實(shí)驗(yàn)中，分別選取了不同的閾值ε為0.2、0.3、0.5、0.6、0.8 和 1.0 對(duì) TSMI-GRN 進(jìn)行了測(cè)試．將測(cè)試結(jié)果與 KEGG數(shù)據(jù)庫(kù)中已有的結(jié)果進(jìn)行比較，發(fā)現(xiàn)0.5ε=時(shí)的實(shí)驗(yàn)結(jié)果與KEGG數(shù)據(jù)庫(kù)中的結(jié)果最接近．圖 2比較了選取不同閾值時(shí)網(wǎng)絡(luò)的準(zhǔn)確度．

圖2 不同閾值下構(gòu)建網(wǎng)絡(luò)的準(zhǔn)確度Fig.2 Accuracy of the constructed network with different thresholds

在基因調(diào)控網(wǎng)絡(luò)的結(jié)構(gòu)學(xué)習(xí)過(guò)程中，為了定量評(píng)價(jià)網(wǎng)絡(luò)構(gòu)建的準(zhǔn)確性，使用靈敏度(sensitivity)和特異度(specificity)作為評(píng)價(jià)指標(biāo)．為此，需要首先計(jì)算以下4個(gè)值．

(1)真陽(yáng)性TP(true positive)：構(gòu)建的網(wǎng)絡(luò)中正確的邊數(shù)．

(2)真陰性 TN(true negative)：構(gòu)建的網(wǎng)絡(luò)中正確的不存在的邊數(shù)．

(3)假陽(yáng)性 FP(false positive)：構(gòu)建的網(wǎng)絡(luò)中錯(cuò)誤的邊數(shù)．

(4)假陰性 FN(false negative)：構(gòu)建的網(wǎng)絡(luò)中錯(cuò)誤的不存在的邊數(shù)．

由此，可以定義靈敏度和特異度的計(jì)算式為

顯然，靈敏度和特異度的值越高，網(wǎng)絡(luò)的構(gòu)建結(jié)果就越好．

圖 3為在含有 14個(gè)基因、24個(gè)時(shí)間點(diǎn)的微陣列數(shù)據(jù)集上用 TSML-GRN算法構(gòu)建的基因調(diào)控網(wǎng)絡(luò)．與KEGG數(shù)據(jù)庫(kù)中的網(wǎng)絡(luò)相比，本文構(gòu)造的網(wǎng)絡(luò)中有11條正確的邊，4條添加的邊．

圖3 TSMI-GRN算法構(gòu)建的14個(gè)基因調(diào)控網(wǎng)絡(luò)Fig.3 Regulation network containing fourteen genes constructed with TSMI-GRN algorithm

本文還在上述數(shù)據(jù)集上測(cè)試了 MH算法[2]和BDE打分爬山搜索算法[2]，表1比較了3種算法構(gòu)建的基因調(diào)控網(wǎng)絡(luò)的靈敏度與特異度．

表1 靈敏度與特異度Tab.1 Sensitivity and specificity %

根據(jù)表 1，可以看出基于時(shí)序數(shù)據(jù)的互信息算法的靈敏度高于其他兩種算法，說(shuō)明 TSMI-GRN 算法能夠更準(zhǔn)確地構(gòu)建基因調(diào)控網(wǎng)絡(luò)．

圖4為在酵母菌周期細(xì)胞的實(shí)驗(yàn)數(shù)據(jù)上，測(cè)試TSMI-GRN算法構(gòu)建的基因調(diào)控網(wǎng)絡(luò)．選取的實(shí)驗(yàn)數(shù)據(jù)集含有 25個(gè)基因、24個(gè)時(shí)間點(diǎn)．運(yùn)用 TSMI-GRN算法在該數(shù)據(jù)集上構(gòu)建的基因調(diào)控網(wǎng)絡(luò)，與數(shù)據(jù)庫(kù)KEGG中的網(wǎng)絡(luò)相比，其中有14條正確的邊，6條方向相反的邊，8條添加的邊．在8條添加的邊中，發(fā)現(xiàn)Cdc28與 Cdc20之間的調(diào)控關(guān)系在文獻(xiàn)[9]實(shí)驗(yàn)中被發(fā)現(xiàn)，Chk1與 Rad9之間的調(diào)控關(guān)系在文獻(xiàn)[10]實(shí)驗(yàn)中被發(fā)現(xiàn)，對(duì)于學(xué)習(xí)到的其他反向邊與添加的邊，對(duì)生物學(xué)實(shí)驗(yàn)具有一定的指導(dǎo)作用．以上的實(shí)驗(yàn)結(jié)果表明，TSMI-GRN算法是有效的．

圖4 TSMI-GRN算法構(gòu)建的包含25個(gè)基因的調(diào)控網(wǎng)絡(luò)Fig.4 Regulation network containing twenty five genes constructed with TSMI-GRN algorithm

4 結(jié) 語(yǔ)

綜合前面的實(shí)驗(yàn)結(jié)果與分析可見(jiàn)，TSMI-GRN算法可以從時(shí)間序列微陣列數(shù)據(jù)中更準(zhǔn)確地構(gòu)建基因調(diào)控網(wǎng)絡(luò)．該算法的優(yōu)點(diǎn)為：①不需要將基因表達(dá)數(shù)據(jù)離散化，與基因表達(dá)數(shù)據(jù)的連續(xù)性相符合，因而可以更準(zhǔn)確地構(gòu)建基因調(diào)控網(wǎng)絡(luò)；②在計(jì)算兩個(gè)基因的互信息時(shí)，使用時(shí)序數(shù)據(jù)的互信息，與傳統(tǒng)的互信息比較，增加了時(shí)間的特性，與時(shí)間序列的數(shù)據(jù)更吻合；③TSMI-GRN算法利用協(xié)方差矩陣計(jì)算兩個(gè)基因間時(shí)序互信息，并且考慮了其他所有基因?qū)@兩個(gè)基因間互信息的影響；④實(shí)驗(yàn)結(jié)果證明，時(shí)序互信息算法的靈敏度高于MH算法以及BDE打分爬山搜索算法的靈敏度；⑤TSMI-GRN算法發(fā)現(xiàn)了 Cdc28與Cdc20、Chk1與Rad9之間的調(diào)控關(guān)系，這在一些生物學(xué)實(shí)驗(yàn)中得到了驗(yàn)證．因而 TSMI-GRN算法是一種有效的算法，在構(gòu)建基因調(diào)控網(wǎng)絡(luò)中發(fā)揮積極作用．

經(jīng)過(guò)實(shí)驗(yàn)驗(yàn)證，TSMI-GRN算法是一種有效的算法，但 TSMI-GRN算法有一定的局限性，它假定微陣列數(shù)據(jù)集服從高斯分布，計(jì)算兩個(gè)基因間的互信息時(shí)用線性熵來(lái)近似代替信息熵．

［1］ 劉萬(wàn)霖，李 棟，朱云平，等．基于微陣列數(shù)據(jù)構(gòu)建基因調(diào)控網(wǎng)絡(luò)[J].遺傳，2007，29(12)：1434-1442．

Liu Wanlin，Li Dong，Zhu Yunping，et al. Constructing gene regulation network from microarray data[J]. Hereditas，2007，29(12)：1434-1442(in Chinese).

［2］ Wu Huihai，Liu Xiaohui. Dynamic Bayesian networks modeling for inferring genetic regulatory networks by search strategy：Comparison between greedy hill climbing and MCMC methods[J]. International Journal of Computational Intellgence，2009，5(2):189-199.

［3］ Zhu Dongxiao，Li Hua. Improved Bayesian network inference using relaxed gene ordering[J]. International Journal of Data Mining and Bioinformatics，2010，4(1)：44-59.

［4］ 王開(kāi)軍，張軍英，趙 峰，等．幾何模式動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)推理基因調(diào)控網(wǎng)絡(luò)[J].西安電子科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2007，34(6)：922-925．

Wang Kaijun，Zhang Junying，Zhao Feng，et al. Geometric-pattern dynamic Bayesian networks reasoning gene regulatory networks[J]. Journal of Xidian University：Natural Science，2007，34(6)：922-925(in Chinese).

［5］ Liu Tiefei，Sung Wingkin，Mittal Ankush. Model gene network by semi-fixed Bayesian network [J]. Expert Systems with Applications，2006，30(1)：42-49.

［6］ Dojer N，Gambin A. Mizera A，et al. Applying dynamic Bayesian networks to perturbed gene expression data[J].BMC Bioinformatics，2006，7：249-259.

［7］ Chan Z S H， Collins L， Kasaboy N. Bayesian learning of sparse gene regulatory networks[J]. Biosystems，2007，87(2/3)：299-306.

［8］ Liang S，F(xiàn)uhrman S，Somogyi R. Reveal，a general reverse engineering algorithm for inference of genetic network architectures[J]. Pacific Symposium on Biocomputing，1998，3：18-29.

［9］ Rudner A D，Murray A W. Phosphorylation by Cdc28 activates the Cdc20-dependent activity of the anaphasepromoting complex[J]. The Journal of Cell Biology，2000，149(7)：1377-1390.

［10］ Loegering D，Arlander S J. Rad9 protects cells from topoisomerase poison-induced cell death[J]. The Journal of Biological Chemistry，2004，279(18)：18641-18647.