基于隨機森林回歸的網絡構建方法及應用*

哈爾濱醫科大學衛生統計學教研室（150086） 侯 艷 楊 凱 李 康

基于隨機森林回歸的網絡構建方法及應用*

哈爾濱醫科大學衛生統計學教研室（150086） 侯 艷 楊 凱 李 康△

目的探討基于隨機森林（RF）回歸估計因果關系網絡的效果。方法通過模擬實驗設定因果關系網絡，對數據標準化后，利用全條件RF回歸對其進行估計并評價其準確性。另外將該方法用于卵巢癌基因表達譜數據，并對分析結果進行驗證。結果模擬實驗結果表明RF回歸對于預先設定網絡關系的識別能力明顯優于貝葉斯網絡方法。當選擇合適的閾值時，隨著樣本含量的增加基于隨機森林回歸方法構建的網絡準確性不斷提高，但傳統經典的貝葉斯方法效果基本保持不變；實例分析結果驗證，基于RF回歸方法能夠得到與現有數據庫的網絡結構。結論應用基于RF回歸方法估計的網絡，能夠在樣本量較少的情況下得出準確度較高的網絡。

調控網絡 隨機森林回歸 貝葉斯網絡

目前有多種方法對基因調控、蛋白互作及代謝通路等網絡進行估計，例如通過計算變量間偏相關系數推斷變量間的條件相關關系、基于概率的方法確定信息網絡邊［1］、基于圖形理論和信息傳遞算法獲得網絡中的直接信息流［2－3］、利用因果推斷模型來獲取變量的相關關系，如貝葉斯網絡模型［4］等。當前的一些研究顯示，貝葉斯方法可能是揭示復雜的細胞調控網絡的最有效方法［5］。然而，貝葉斯方法使用的變量不能過多，且局限于低階交互作用，因此在實際應用過程中受到一定的限制。本文在簡要介紹隨機森林（RF）回歸構建網絡方法的基礎上，通過模擬實驗研究這種方法的適用性和特點，與普通的貝葉斯方法進行比較，最后給出應用實例。

原理與方法

1．基本思想

RF回歸是一種推斷網絡結構的新算法［6］，其基本思想是對網絡隨機向量Xnet中m個變量中每個變量做應變量，其余m－1個變量作為自變量做RF回歸分析，在全條件下得到任意兩變量之間連接的權重，排序后按照一定的閾值取連接的邊和方向，推斷出需要估計的網絡。由于RF回歸不會出現過擬合的問題，因此可以在高維情況下進行分析；同時RF回歸對數據的分布和變量之間的關系不需要做出任何假定，能夠挖掘出具有交互作用的變量和各種復雜的非線性關系［9］。

2．網絡變量關聯計算

網絡變量關聯可以用通過RF回歸得到的變量重要性測量值（VIM）來衡量?；跇浞椒ǖ淖兞恐匾詼y量方法有多種，本文使用方差改變量法，即對回歸樹的每一個節點t，計算由變量分裂導致的輸出變量方差的總減少量，定義為

n（t）表示到達節點t的樣本量，Var（t）為應變量的樣本方差；n（tL）和n（tR）分別表示由節點t分裂的兩新節點的樣本量，即有n（t）＝n（tL）＋n（tR）；Var（tL）和Var（tR）為兩新節點的樣本方差。對于一棵樹，一個變量的重要性可以通過用這個變量分裂的所有節點的I（·）值相加獲得，RF回歸變量Xk的重要性可以用b棵樹I（·）的平均值進行衡量，用公式表示為

3．網絡連接權重的計算

VIMkj表示在全自變量X－k條件下應變量Xk（k＝1，2，…，m）與自變量Xj（j≠k）的關聯度，因此兩變量連接的權重wkj可以通過關聯度VIMkj得到，即wkj＝VIMkj。為了使m個不同RF回歸模型中的VIM值具有可比性，需要預先對原始數據進行標準化，使所有變量有相等的方差，從而使從不同模型得到的權重具有可比性。對數據作m次RF回歸，得到如下矩陣：

規定上式j＝k時VIMkj＝0，即對角線元素為零。

可以證明，在RF回歸中所有變量重要性的總和等于因變量中被其解釋的總方差，在未修剪樹的情況下則通常非常接近因變量的總方差，可用公式表示為

其中S是用來估計VIM值的樣本數據。因此在數據標準化的情況下，不同RF回歸模型的VIM值具有一定的可比性，從而可以對矩陣VIMkj所有的值排序，再通過設定的閾值c（VIM）（或通過確定連接的邊數），取排序中｛VIMkj≥c（VIM）｝估計出網絡連接的邊和方向。如何給出合適的閾值是一個需要研究的問題。

上述算法實現可以使用random Forest R軟件包。

模擬研究

1．無噪聲情況的網絡估計

（1）模擬實驗目的和條件設置：研究網絡變量之間存在線性關系及非線性相關關系時，基于RF方法構建相關網絡的準確性，同時與經典的貝葉斯網絡構建方法進行比較。實驗設置15個變量和16條有向邊（圖1），選取的網絡中包括：一個變量調控多個變量，如G12→｛G14，G17｝；多個變量調控一個變量｛G4，G9，G13｝→G8，線性級聯G11→G9→G10，前饋環路｛G2，G3，G4｝。

（2）網絡參數及數據：變量之間的調控分別具有線性關系和交互關系，即｛G13，G9｝→G8和｛G2，G4｝→G3為交互作用，其他為線性調控關系。誤差服從正態分布。樣本量分別設為50、100、200、500和1000。

圖1 模擬實驗設置的網絡結構

（3）網絡數據模型：基于RF回歸構建網絡，分別設定VIM值為前15、20和25條邊作為陽性邊，剩余的邊作為陰性邊；對于貝葉斯網絡，取結果為1的邊作為陽性邊，結果為0的邊作為陰性邊，分別使用AUC值及預測準確率（PRE）進行比較，對兩種網絡構建方法構建的網絡進行評價。以上過程隨機重復100次。

（4）模擬試驗結果

表1模擬結果顯示，當選擇隨機森林的邊數（閾值）為15時（接近真實邊的數目），在準確性評價中，AUC和PRE兩項評價指標明顯優于其他兩種情況，同時明顯優于貝葉斯方法；但是，當選定的邊數（閾值）明顯大于真實邊數時（25條邊），其準確性明顯下降。當樣本量小于200時，RF和貝葉斯建模的整體準確性不夠穩定，在大于200時隨樣本量的增加趨于平緩（圖2）。

2．具有噪聲情況下的網絡估計

模擬實驗目的和條件設置：研究基于隨機森林回歸方法在網絡變量間存在噪聲變量時的識別能力，同時與不存在噪聲變量的模擬數據分析效果進行比較。真實網絡變量的條件設置同前，另加入200個服從正態分布的噪聲變量。模擬隨機重復100次。

模擬結果見圖3。結果表明：利用隨機森林回歸估計網絡的方法與不加入噪聲變量的結果相近，即加入噪聲變量后仍然能夠很好地識別出變量之間的調控關系，并且保持較好的準確性。

3．閾值選擇

上述兩個模擬實驗網絡邊數的選擇，使用了選擇15、20和25條邊三種情況的VIM閾值。為了能夠自動選擇最合適的閾值，現使用隨機置換的方法，在樣本量為500例的情況下，打亂上述實驗設置的15個變量的模擬數值（無噪聲情況），使用隨機森林回歸計算相應的VIM值，進行100次置換。置換后所得到的VIM值可以得到隨機情況下VIM值的分布，其99%分位數為0．1034，即選擇閾值c（VIM）＝0．1034。100次模擬數據中，變量調控關系大于閾值0．1034的平均有15條，這一結果與實驗設置16條邊的真實情況極為相近。由此可以得出，在實際應用中，在不知道真實調控關系的數量時，可以通過隨機置換的方法估計真實調控關系的數量。需要注意的是，置換檢驗閾值的選擇，可根據具體情況決定，如果需要控制“邊”的數量，可選擇比較嚴格的標準（如P99）；若想適當放寬選入邊的數量，則可以選擇P95作為閾值。

表1 基于隨機森林回歸網絡構建方法與貝葉斯網絡對網絡的識別效果（±S）

表1 基于隨機森林回歸網絡構建方法與貝葉斯網絡對網絡的識別效果（±S）

評價指標樣本量隨機森林網絡構建方法15條邊20條邊25條邊貝葉斯網絡AUC 50 0．760±0．058 0．727±0．048 0．705±0．038 0．745±0．045 100 0．838±0．048 0．801±0．036 0．757±0．025 0．800±0．033 200 0．896±0．034 0．839±0．027 0．782±0．017 0．829±0．034 500 0．945±0．022 0．865±0．016 0．792±0．012 0．842±0．027 1000 0．952±0．017 0．868±0．013 0．795±0．009 0．852±0．023 PRE 50 0．583±0．104 0．508±0．082 0．454±0．062 0．520±0．087 100 0．715±0．086 0．624±0．063 0．529±0．045 0．626±0．065 200 0．816±0．062 0．687±0．048 0．570±0．032 0．681±0．067 500 0．902±0．043 0．731±0．030 0．594±0．027 0．705±0．053 1000 0．916±0．034 0．744±0．029 0．612±0．025 0．724±0．045

圖2 基于隨機森林回歸網絡構建方法與貝葉斯網絡估計網絡的準確性

圖3 基于隨機森林回歸的網絡構建方法抗噪能力情況

實例分析

數據來源：通過TCGA數據庫下載570例卵巢癌患者以及8例健康對照數據的基因表達譜數據［7］，其中包含12042個基因的表達值。

為了簡化分析過程，首先篩選出與卵巢癌相關的基因，再對這部分基因構建網絡。本研究使用基于W ilcoxon秩和檢驗的置換檢驗，進行1000次置換，篩選出P＜0．05（校正后）的基因一共744個，通過對這部分基因進行KEGG通路富集分析，結果有12個基因顯著富集在p53信號通路中。通過置換檢驗方法，100次隨機置換后的VIM值的99%分位數為0．1157。確定邊數選擇的閾值c（VIM）＝0．1157，獲得16條可能具有調控關系的邊，網絡如圖4所示。

進而，通過查詢已有的基因／蛋白互作網絡數據庫GeneMANIA［8］，發現這12個基因中，有9條邊的因果關系出現在GeneMANIA的網絡中（圖5）。在這9條邊中，基因CCNB1、CCNB2、CHEK1和CHEK2參與了比較多的調控，這個結果與隨機森林回歸所構建網絡的結構相一致，說明這幾個基因在p53信號通路中起到了很重要的作用。在剩余的7條邊中，有4條邊是與基因GTSE1（即PRKAR1A）有關，而在GeneMANIA中并沒有相應的調控關系，提示這4條邊需要進一步研究。

圖4 卵巢癌患者p53信號通路中基因網絡結構

圖5 12個基因在GeneMANIA中的網絡關系

討 論

1．基于RF回歸估計網絡的方法，兩變量連接的權重wkj主要是利用全自變量條件下標準化的VIMkj得到。實際中，基于樹方法的變量重要性測量方法有多種，本文使用方差改變量法。另外，VIM值的計算還有隨機置換法，這種方法與方差改變量法相比在理論上有一定的優勢，但在實際應用中得到的結果可能相近，而且需要更多的計算。

2．理論上，RF回歸估計網絡對數據的分布和變量之間的關系不需要做出任何假定，能夠挖掘出具有交互作用的變量和各種復雜的非線性關系；而貝葉斯網絡模型更適合單調或線性的調控關系。貝葉斯模型的優點之一是能夠給出調控關系的方向，但模擬實驗表明RF回歸估計網絡方法同樣能夠給出調控的方向，其原理是對任一變量，在全自變量條件下，能夠對其方差做出較多解釋的變量最有可能是調控其變化的變量。相對而言，對于一個變量調控多個變量更容易識別，如圖1中的G12→｛G17，G14｝；而對于多個變量調控一個變量則相對容易判錯方向，如｛G4，G9，G13｝→G8。

3．RF可以在高維（如m＞500）情況下構建網絡，因此可以用于組學（如基因組／蛋白組）研究。同時RF回歸方法可以分析非線性復雜關系的網絡，附加的模擬實驗結果表明，在變量之間呈指數、多項式等關系時，使用這種方法仍能夠獲得理想的結果。另外，這種方法也可以很容易推廣到時間序列數據的分析。

4．關于RF回歸參數的設置，在隨機森林中，生成每一棵分類樹時，所應用的自助樣本集從原始的訓練樣本集中隨機選取，每一棵樹所應用的變量也是從所有變量mall中隨機選取，并用袋外數據來衡量回歸的效果。隨機森林中最重要的參數有兩個，即建立回歸樹的個數ntree和建立每個節點所取的變量兩數目ntry。理論上講，ntree越大得到的回歸森林越穩定，本文取ntree＝1000；ntry的選取主要有和ntry＝mall－1兩種方法，因本文變量數目不是很大，取ntry＝mall－1。

5．實例分析結果表明，網絡中多數的連接邊能夠被識別出來，并且能夠給出調控的方向。然而，需要注意的是，實際各種變量之間的關系可能有多種，如基因共表達、基因調控、代謝通路等，其中共表達為無方向性，而RF回歸方法同樣能夠對此進行分析，如VIMij≈VIMji時其方向不明確。需要注意的是，對于斷面數據得到的因果網絡關系最終需要實驗的方法加以證實。

［1］Butte AJ，Kohane IS．Mutual information relevance networks：functional genomic clustering using pairwise entropy measurements．Pac Symp Biocomput，2000，418-429．

［2］Pesch R，Lysenko A，Hindle M，et al．Graph-based sequence annotation using a data integration approach．Journal of integrative bioinformatics，2008，5：2

［3］Harley E，Bonner A，Goodman N．Uniform integration of genome mapping data using intersection graphs．Bioinformatics，2001，17（6）：487-494．

［4］Mani S，Cooper G F．A Bayesian local causal discovery algorithm，2004：731-735．

［5］Margolin AA．ARACNE：an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context．BMC bioinformatics，2006，7（suppl．1）：S7．

［6］Breiman．L．Random Forests．Machine learning，2001，45：5-32．

［7］Integratedgenom ic analyses of ovarian carcinoma．Nature，2011，474（7353）：609-615．

［8］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．

（責任編輯：郭海強）

Network Reconstruction with Random Forest Regression and its Application

Hou Yan，Yang Kai，Li Kang（DepartmentofHealthStatistics，SchoolofPublicHealth，HarbinMedicalUniversity（150086），Harbin）

ObjectiveTo investigate the performance of network reconstruction based on random forest regression．MethodsSimulation studies were performed to evaluate the accuracy for network reconstruction with standardized data and conditional random forest regression．ResultsSimulation studies demonstrated that the network reconstruction performance with random forest regression is better than thatwith Bayesian network．In particular，when the thresholds are selected appropriately，the performance for network reconstruction based on random forest regression could improve with the increase of sample size while the traditional Bayesian network w ill remain stable．Besides，we applied this approach to the realexample and achieved satisfactory performance．ConclusionThe proposed method in this paper could achieve satisfactory performance for network reconstruction in small sample size．

Regulatory network；Random forest regression；Bayesian network

國家自然科學基金（81473072）；中國博士后面上項目（2015M571445）

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