基于bootstrap方法的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法在構(gòu)建基因調(diào)控網(wǎng)絡(luò)中的應(yīng)用*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計教研室（150081） 李海龍 侯 艷 柯朝甫 李 康

基于bootstrap方法的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法在構(gòu)建基因調(diào)控網(wǎng)絡(luò)中的應(yīng)用*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計教研室（150081） 李海龍 侯 艷 柯朝甫 李 康△

目的探討基于bootstrap重抽樣方法的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法構(gòu)建網(wǎng)絡(luò)的性能，并將其應(yīng)用于卵巢癌基因表達譜數(shù)據(jù)分析。方法通過模擬實驗和實例驗證本文給出的算法構(gòu)建網(wǎng)絡(luò)的有效性，同時將這種算法應(yīng)用于構(gòu)建基因調(diào)控網(wǎng)絡(luò)。結(jié)果模擬實驗顯示，在樣本量較小的情況下，基于bootstrap算法構(gòu)建的貝葉斯網(wǎng)絡(luò)明顯優(yōu)于普通貝葉斯方法構(gòu)建的網(wǎng)絡(luò)；實例分析結(jié)果也表明，應(yīng)用本文的方法能夠得到有價值的網(wǎng)絡(luò)結(jié)構(gòu)。結(jié)論應(yīng)用本文給出的算法能夠在樣本量較少的情況下得出準(zhǔn)確度較高的網(wǎng)絡(luò)，同時能夠給出網(wǎng)絡(luò)結(jié)構(gòu)中各條邊置信度的估計值。

貝葉斯網(wǎng)絡(luò) 結(jié)構(gòu)學(xué)習(xí) bootstrap

貝葉斯網(wǎng)絡(luò)是一種概率圖形模型，它能夠發(fā)現(xiàn)變量之間潛在的依賴關(guān)系。其模型構(gòu)建可分為三個步驟：①網(wǎng)絡(luò)變量的確定；②網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)；③參數(shù)估計［1－2］。貝葉斯網(wǎng)絡(luò)的結(jié)構(gòu)學(xué)習(xí)是根據(jù)原始數(shù)據(jù)，通過一定的搜索策略找到得分最高的網(wǎng)絡(luò)結(jié)構(gòu)，得分高說明網(wǎng)絡(luò)結(jié)構(gòu)能夠很好地代表數(shù)據(jù)中變量間的調(diào)控關(guān)系［3］。然而，實際中由于樣本量不足，常常出現(xiàn)一些結(jié)構(gòu)不同而得分相近的網(wǎng)絡(luò)，難以從得分相近的網(wǎng)絡(luò)中分辨出哪一種結(jié)構(gòu)更接近真實網(wǎng)絡(luò)［4－5］。此外，一般的結(jié)構(gòu)學(xué)習(xí)方法難以根據(jù)評價指標(biāo)來評價網(wǎng)絡(luò)結(jié)構(gòu)的可靠程度。本文將貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)方法與bootstrap重抽樣方法相結(jié)合，通過設(shè)定閾值得到包含高置信度邊的網(wǎng)絡(luò)，并與一般的結(jié)構(gòu)學(xué)習(xí)方法相比較，考察其有效性。最后運用本文給出的方法對卵巢癌基因表達譜數(shù)據(jù)進行分析，做出生物學(xué)解釋。

原理與方法

1.貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)

貝葉斯網(wǎng)絡(luò)是一個有向無環(huán)圖，可以表示成一組隨機變量的聯(lián)合概率分布。形式上一組隨機變量X＝｛X1，…，Xn｝的貝葉斯網(wǎng)絡(luò)可以用B＝（G，θ）表示，其中第一個成分θ表示一個有向無環(huán)圖，圖中節(jié)點代表隨機變量，節(jié)點之間的邊代表變量之間的直接依賴關(guān)系。第二個成分θ，代表一組量化網(wǎng)絡(luò)的參數(shù)θ＝（θ1，θ1，…，θm′），m′＞m以條件概率分布的形式表示，即θi＝PB（Xi｜pa（Xi）），其中pa（Xi）表示變量Xi在圖G中的父節(jié)點集。貝葉斯網(wǎng)絡(luò)B給一組變量X定義的聯(lián)合概率分布：

貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)可以歸結(jié)為：對于給定的數(shù)據(jù)訓(xùn)練集D，尋找一個網(wǎng)絡(luò)B使之能與數(shù)據(jù)集D最匹配。解決這個問題最常用的方法就是引入一個得分函數(shù)來評價對應(yīng)訓(xùn)練集所得網(wǎng)絡(luò)的擬合程度，然后根據(jù)得分搜索到最優(yōu)網(wǎng)絡(luò)。本文采用BIC得分函數(shù)，并運用貪婪爬山法（greed hill-climbing）結(jié)合隨機重搜索得分最高的網(wǎng)絡(luò)結(jié)構(gòu)，這種方法能夠避免陷入局部最優(yōu)。

網(wǎng)絡(luò)的得分使用BIC準(zhǔn)則確定，BIC得分越大，構(gòu)建的網(wǎng)絡(luò)越好，其計算公式為

其中N為數(shù)據(jù)的總例數(shù)，d為網(wǎng)絡(luò)的參數(shù)個數(shù)。

2.bootstrap方法置信度估計

對于網(wǎng)絡(luò)G的結(jié)構(gòu)，感興趣的特征可以是某條有向邊X→Y，也可以是無向邊X－Y。總之，可以將這些邊用字母fij來表示，并通過網(wǎng)絡(luò)結(jié)構(gòu)的函數(shù)轉(zhuǎn)換成集合｛0，1｝表示，fij＝0表示節(jié)點Xi和節(jié)點Xj不連接，fij＝1表示兩節(jié)點連接，簡記為f。

PN（f）表示貝葉斯網(wǎng)絡(luò)B中抽到一個任意兩節(jié)點是否相連網(wǎng)絡(luò)的概率。如果結(jié)構(gòu)學(xué)習(xí)過程一致，則希望當(dāng)樣本量N足夠大時，pN（f）會收斂于f（G）。也就是說，如果真實網(wǎng)絡(luò)結(jié)構(gòu)G中確實存在節(jié)點相連特征f，則它的置信度應(yīng)該接近1，相反如果不存在則應(yīng)該接近于0。

使用bootstrap估計置信度的方法是通過對數(shù)據(jù)集有放回地重抽樣，然后通過對多個bootstrap數(shù)據(jù)集進行學(xué)習(xí)得出多個網(wǎng)絡(luò)，在這多個網(wǎng)絡(luò)中任意兩節(jié)點相連接（包括方向）的頻率就是其置信度估計。算法的過程如下：

M（bootstrap重抽樣次數(shù)），F(xiàn)s（得分函數(shù)），t（閾值）

Output：G，包含概率大于閾值有向邊Xi→Xj的圖形Fori＝1 toMdo

有放回地從數(shù)據(jù)D中抽取N個觀測得到數(shù)據(jù)集Di

根據(jù)Di通過得分函數(shù)Fs指導(dǎo)的學(xué)習(xí)算法得出得分最高的網(wǎng)絡(luò)結(jié)構(gòu)Gi

end

模擬實驗

模擬數(shù)據(jù)來源于已知的真實網(wǎng)絡(luò)結(jié)構(gòu)，目的是檢驗bootstrap平均模型的有效性，即將bootstrap方法中高置信度特征與真實網(wǎng)絡(luò)中的特征進行比較，若bootstrap平均模型包含了大部分原真實網(wǎng)絡(luò)中存在的邊，并具有較高置信度，則能說明該方法的有效性。

1.模擬實驗1

根據(jù)已給定的網(wǎng)絡(luò)結(jié)構(gòu)產(chǎn)生相應(yīng)的模擬數(shù)據(jù)（參見圖1）。網(wǎng)絡(luò)包含7個隨機變量（節(jié)點）和7條邊，根據(jù)此網(wǎng)絡(luò)結(jié)構(gòu)產(chǎn)生10000個觀測，其中變量均服從正態(tài)分布。每次從數(shù)據(jù)集中隨機抽出100例樣本作為結(jié)構(gòu)學(xué)習(xí)的數(shù)據(jù)集，重復(fù)抽樣得到100個數(shù)據(jù)集，分別用典型的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)方法以及貝葉斯網(wǎng)絡(luò)的bootstrap方法分別對一個數(shù)據(jù)集進行結(jié)構(gòu)學(xué)習(xí)得出網(wǎng)絡(luò)，重復(fù)實驗100次。

采用基于信息準(zhǔn)則的BIC得分函數(shù)確定最優(yōu)網(wǎng)絡(luò)［4］，搜索過程采用貪婪爬山法，bootstrap重抽樣次數(shù)為300次。為了避免陷入局部最優(yōu)，在搜索過程中結(jié)合隨機重啟搜索。通過這個過程嘗試尋找能使得分提高最多的網(wǎng)絡(luò)結(jié)構(gòu)，直到結(jié)構(gòu)的改變無法繼續(xù)提高得分為止。一旦爬山法陷入局部最優(yōu)，算法將隨機擾動網(wǎng)絡(luò)結(jié)構(gòu)中的邊（添加、刪除和反向）并重新開始搜索。在重啟一定次數(shù)后終止搜索，選出得分最高的網(wǎng)絡(luò)作為結(jié)果。最后，根據(jù)設(shè)定的三個不同的閾值t＝0.5，0.7，0.9，將pN（f）≥t的所有連接邊輸出得到最終結(jié)果網(wǎng)絡(luò)。模擬使用R軟件包bnlearn［6］和編程實現(xiàn)。

評價構(gòu)建網(wǎng)絡(luò)的指標(biāo)分別使用真陽性數(shù)目、假陽性數(shù)目、假陰性數(shù)目、真陰性數(shù)目、靈敏度、特異度和準(zhǔn)確度，其中準(zhǔn)確度為真陽性邊占陽性邊的比例，相當(dāng)于診斷試驗中的陽性預(yù)測值。計算這些指標(biāo)時需要對100次實驗的結(jié)果取平均值，表1給出了使用不同方法和取不同閾值的網(wǎng)絡(luò)評價結(jié)果，即分別使用普通的貝葉斯方法（origin）和取不同閾值t的基于bootstrap的貝葉斯方法。

圖1 模擬實驗1的真實網(wǎng)絡(luò)關(guān)系圖

表1 使用不同方法和取不同閾值的網(wǎng)絡(luò)評價結(jié)果

2.模擬實驗2

使用ICU-Alarm網(wǎng)絡(luò)模擬數(shù)據(jù)。該數(shù)據(jù)產(chǎn)生于ICU-Alarm網(wǎng)絡(luò)模型，此模型是機器學(xué)習(xí)中網(wǎng)絡(luò)學(xué)習(xí)問題的經(jīng)典模型，廣泛應(yīng)用于評價網(wǎng)絡(luò)學(xué)習(xí)方法。ICU-Alarm網(wǎng)絡(luò)模型包含37個隨機變量，46條邊。在樣本量N＝100，300，600，1000下比較網(wǎng)絡(luò)學(xué)習(xí)的結(jié)果，每個樣本量下抽取100個數(shù)據(jù)集，重復(fù)實驗100次。然后，分別使用普通的貝葉斯網(wǎng)絡(luò)模型和基于bootstrap的貝葉斯方法構(gòu)建網(wǎng)絡(luò)，并對其進行評價。

模擬實驗評價結(jié)果見圖2。結(jié)果顯示：使用基于bootstrap的貝葉斯方法得到網(wǎng)絡(luò)模型明顯優(yōu)于使用普通貝葉斯網(wǎng)絡(luò)模型。同時可以看到，當(dāng)樣本量增加時，構(gòu)建的網(wǎng)絡(luò)的結(jié)構(gòu)學(xué)習(xí)越來越準(zhǔn)確，即真陽性邊增加，假陽性和假陰性的邊減少；另外，提高閾值，真陽性和假陽性邊減少，但容易漏掉真實邊，說明合理設(shè)定閾值的必要性。由于真實邊（46條）相對于網(wǎng)絡(luò)所有可能邊（1332條）要少很多，因此不同方法的特異度均接近于1，假陽性率均低于3%。

實例分析

為了研究卵巢癌的分子生物學(xué)機制，本研究通過對卵巢癌患者基因表達數(shù)據(jù)進行分析并構(gòu)建貝葉斯網(wǎng)絡(luò)，從網(wǎng)絡(luò)中得出基因之間的調(diào)控關(guān)系，并結(jié)合生物功能和通路數(shù)據(jù)庫查詢以及查閱文獻，對網(wǎng)絡(luò)進行生物學(xué)解釋，從基因組學(xué)的角度為卵巢癌的發(fā)病機制提供線索［7］。

本研究從TCGA數(shù)據(jù)庫下載570例卵巢癌患者基因表達譜數(shù)據(jù)，以及8例健康對照數(shù)據(jù)［8］。全基因組表達譜數(shù)據(jù)一共測得12042個基因的表達值，由于基因的數(shù)目過多，需要先篩選出與卵巢癌相關(guān)的基因，再對這部分基因構(gòu)建貝葉斯網(wǎng)絡(luò)。對分析變量的篩選不僅能提高建模的效率，也使構(gòu)建的網(wǎng)絡(luò)更加合理，有助于對其進行解釋。本研究使用基于Wilcoxon秩和檢驗的置換檢驗［8］，進行1000次置換，篩選出P＜0.05（校正后）的基因一共744個。繼而，對這部分基因進行KEGG通路富集分析，結(jié)果有12個基因顯著富集在p53信號通路中。

將映射上這個通路的12個基因的表達數(shù)據(jù)整理出來，并對數(shù)據(jù)構(gòu)建貝葉斯網(wǎng)絡(luò)。貝葉斯網(wǎng)絡(luò)的搜索過程采用貪婪爬山搜索法，再結(jié)合bootstrap重抽樣方法對網(wǎng)絡(luò)特征進行置信度估計，重抽樣次數(shù)設(shè)為1000次以保證結(jié)果的穩(wěn)定性。將閾值設(shè)定為0.8，結(jié)果如圖3所示，其中節(jié)點代表富集于p53信號通路中的基因，灰色的節(jié)點代表樞紐基因（Hub Gene），信度大于0.8小于1的邊用虛線表示，信度等于1的邊用實線表示。

圖3 利用bootstrap方法構(gòu)建的卵巢癌基因調(diào)控網(wǎng)絡(luò)

為了驗證bootstrap方法置信度評價的可信程度，本研究通過隨機重排列每個基因的測量值產(chǎn)生一個新的數(shù)據(jù)集。在這樣一個數(shù)據(jù)集中，基因彼此之間是獨立的，所以我們并不期望能從中找出真實的邊。具體做法如下，對每個基因下的所有測量值隨機打亂順序產(chǎn)生新的數(shù)據(jù)集，對此數(shù)據(jù)進行學(xué)習(xí)構(gòu)建貝葉斯網(wǎng)絡(luò)，結(jié)合bootstrap方法得出每條邊的置信度，重復(fù)100次這樣的實驗。結(jié)果見圖4，圖中實線為真實數(shù)據(jù)下不同置信度水平下有向邊的數(shù)目，虛線表示隨機重排列數(shù)據(jù)下有向邊的數(shù)目，橫軸表示置信度閾值，縱軸表示大于等于對應(yīng)置信度閾值的有向邊數(shù)目。正如預(yù)期，對隨機數(shù)據(jù)集構(gòu)建的網(wǎng)絡(luò)中邊的置信度普遍比較低。如圖4所示，比較原始數(shù)據(jù)集和隨機數(shù)據(jù)集在不同置信度下的邊數(shù)，可以看出原始數(shù)據(jù)的邊數(shù)分布在高置信度區(qū)域有更長更重的尾部。當(dāng)置信度大于0.2時，兩條分布曲線出現(xiàn)間隔，隨著置信度的增大間隔也越來越大，即在原始數(shù)據(jù)上得到的網(wǎng)絡(luò)關(guān)系具有一定的可信度，說明貝葉斯網(wǎng)絡(luò)的bootstrap估計方法確實能夠發(fā)現(xiàn)大量的網(wǎng)絡(luò)關(guān)系。

圖3中構(gòu)建的貝葉斯網(wǎng)絡(luò)反映了基因之間的調(diào)控關(guān)系，通過查詢已有的基因／蛋白互作網(wǎng)絡(luò)數(shù)據(jù)庫［9－10］（如STRING，GENEMANIA等），貝葉斯網(wǎng)絡(luò)中基因調(diào)控關(guān)系80%以上得到支持。圖3中RRM 2基因受多個基因調(diào)控，可以將它定義為樞紐基因。已有大量文獻報道該基因與卵巢癌的診斷，預(yù)后以及化療有關(guān)［11－12］，它所編碼的蛋白構(gòu)成氧化還原酶，能催化核苷酸還原成脫氧核苷酸的反應(yīng)，為DNA合成提供前體準(zhǔn)備。圖3中另外一個樞紐基因CHEK1調(diào)控多個基因，該基因與卵巢癌有密切聯(lián)系［13－14］，其編碼的蛋白屬于絲氨酸／蘇氨酸蛋白激酶家族，在DNA損傷反應(yīng)中起著重要作用。

圖4 100次重復(fù)實驗不同閾值下網(wǎng)絡(luò)中邊的數(shù)目比較

討 論

本研究的目的是驗證貝葉斯網(wǎng)絡(luò)的bootstrap估計方法的性能，并將其應(yīng)用于癌癥基因組數(shù)據(jù)的分析，揭示基因之間的調(diào)控關(guān)系。通過模擬實驗將貝葉斯網(wǎng)絡(luò)bootstrap方法與一般結(jié)構(gòu)學(xué)習(xí)方法的結(jié)果進行比較，驗證了改進后方法的性能。此外，貝葉斯網(wǎng)絡(luò)bootstrap方法能給出網(wǎng)絡(luò)中邊的置信度，這可以為研究者提供更多的信息。通過模擬實驗我們檢驗了置信度可以作為評價特征真實性的度量。本研究得出以下幾點結(jié)論：①bootstrap估計是謹(jǐn)慎可靠的，在高置信度的情況下網(wǎng)絡(luò)幾乎不包含假陽性。②當(dāng)數(shù)據(jù)集樣本量較少（相對其所要推斷的模型復(fù)雜度）時，本文給出的bootstrap方法相比原始方法新方法能夠得出更準(zhǔn)確的結(jié)果。③閾值的設(shè)定十分重要，它直接影響最終結(jié)果，要根據(jù)實際情況設(shè)置，如果實際中想發(fā)現(xiàn)更多的網(wǎng)絡(luò)關(guān)系，可選較小的閾值（如t＝0.3），如果想得到更可信的網(wǎng)絡(luò)關(guān)系，則應(yīng)選取大的值（如t＝0.7）。總之，建立生物學(xué)網(wǎng)絡(luò)可以更好地驗證差異變量，揭示變量之間的因果關(guān)系，本文將bootstrap方法應(yīng)用于貝葉斯網(wǎng)絡(luò)估計，獲得了較為理想的結(jié)果，更深入的問題有待進一步研究。

1.游項云，李康.貝葉斯網(wǎng)絡(luò)方法在基因調(diào)控研究中的應(yīng)用.中國衛(wèi)生統(tǒng)計，2009，26（1）：83-86.

2.范麗珺，游頂云，張旺，等.貝葉斯因果關(guān)系網(wǎng)絡(luò)模型在斷面調(diào)查數(shù)據(jù)中的應(yīng)用.中國醫(yī)院統(tǒng)計，2010，17（2）：97-100.

3.虞慧婷，吳騁，柳偉偉，等.基于貝葉斯網(wǎng)絡(luò)的原發(fā)性肝癌預(yù)后影響因素相互關(guān)系研究.中國衛(wèi)生統(tǒng)計，2008，25（1）：10-14.

4.Friedman N，Goldszmidt M，Wyner A.Data analysis with Bayesian networks：A bootstrap approach.Proceedings of the Fifteenth conference on Uncertainty in artificial intelligence.Morgan Kaufmann Publishers Inc.，1999：196-205.

5.Broom BM，Do KA，Subramanian D.Model averaging strategies for structure learning in Bayesian networks with limited data.BMC Bioinformatics，2012，13（Suppl 13）：S10.

6.Scutari M.Learning Bayesian networks with the bnlearn R package. arXiv preprint arXiv：0908.3817，2009.

7.Friedman N，Linial M，Nachman I，Pe’er D.Using Bayesian networks to analyze expression data.J Comput Biol，2000，7（3-4）：601-620

8.Bell D，Berchuck A，Birrer M，et al.Integrated genomic analyses of ovarian carcinoma.Nature，2011，474（7353）：609-615.

9.Warde-Farley D，Donaldson SL，Comes O，et al.The GeneMANIA prediction server：biological network integration for gene prioritization and predicting gene function.Nucleic Acids Res，2010，38（Web Server issue）：W 214-220.

10.Franceschini A，Szklarczyk D，F(xiàn)rankild S，et al.STRING v9.1：proteinprotein interaction networks，with increased coverage and integration. Nucleic Acids Res，2013，41（Database issue）：D808-815.

11.Ferrandina G，Mey V，Nannizzi S，et al.Expression of nucleoside transporters，deoxycitidine kinase，ribonucleotide reductase regulatory subunits，and gemcitabine catabolic enzymes in primary ovarian cancer. Cancer Chemother Pharmacol，2010，65（4）：679-686.

12.Zhang M，Wang J，Yao R，et al.Small interfering RNA（siRNA）-mediated silencing of the M2 subunit of ribonucleotide reductase：a novel therapeutic strategy in ovarian cancer.International Journal of Gynecological Cancer，2013，23（4）：659-666.

13.Connell CM，Shibata A，Tookman LA，et al.Genomic DNA damage and ATR-Chk1 signaling determine oncolytic adenoviral efficacy in human ovarian cancer cells.J Clin Invest，2011，121（4）：1283-1297.

14.Kumar G，Breen EJ，Ranganathan S.Identification of ovarian cancer associated genes using an integrated approach in a Boolean framework. BMC Syst Biol，2013，7：12.

（責(zé)任編輯：鄧 妍）

The Application of Bayes Network Structure Learning Algorithm Based on Bootstrap Method to the Construction of Gene Regulatory Networks

Li Hailong，Hou Yan，Ke Chaofu，et al.（Department of Medical Statistics，Harbin Medical University（150081），Harbin）

ObjectiveTo explore the performance of Bayes network structure learning algorithm based on bootstrap method in network construction，and to apply it to the analysis of ovarian cancer gene expression data.MethodsThe efficiency of the algorithm given in this article was testified with simulation data and gene expression data，and meanwhile this algorithm was used to construct gene regulatory networks.ResultsBayes network structure learning based on bootstrap method performed better than the general Bayes network in the case of small sample sizes，as shown in simulation tests；the results of gene expression data analysis also indicated that this algorithm could provide valuable network structures.ConclusionBayes network structure learning algorithm based on bootstrap method can establish highly precise network models even with small sample sizes，and meanwhile provide the confidence estimates of each edge in the network.

Bayes network；Structure learning；Bootstrap

*高等學(xué)校博士學(xué)科專項基金（2012230711004）；國家自然科學(xué)基金（81172767）

△通信作者：李康，E-mail：likang＠ems.hrbmu.edu.cn