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

基于bootstrap方法的貝葉斯網(wǎng)絡結構學習算法在構建基因調控網(wǎng)絡中的應用*

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

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

貝葉斯網(wǎng)絡 結構學習 bootstrap

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

原理與方法

1.貝葉斯網(wǎng)絡結構學習

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

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

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

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

2.bootstrap方法置信度估計

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

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

使用bootstrap估計置信度的方法是通過對數(shù)據(jù)集有放回地重抽樣，然后通過對多個bootstrap數(shù)據(jù)集進行學習得出多個網(wǎng)絡，在這多個網(wǎng)絡中任意兩節(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指導的學習算法得出得分最高的網(wǎng)絡結構Gi

end

模擬實驗

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

1.模擬實驗1

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

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

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

圖1 模擬實驗1的真實網(wǎng)絡關系圖

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

2.模擬實驗2

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

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

實例分析

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

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

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

圖3 利用bootstrap方法構建的卵巢癌基因調控網(wǎng)絡

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

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

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

討 論

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

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

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

3.虞慧婷，吳騁，柳偉偉，等.基于貝葉斯網(wǎng)絡的原發(fā)性肝癌預后影響因素相互關系研究.中國衛(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.

（責任編輯：鄧 妍）

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

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

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