哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室（150086） 侯 艷 楊 凱 李 康

基于隨機(jī)森林回歸的網(wǎng)絡(luò)構(gòu)建方法及應(yīng)用*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室（150086） 侯 艷 楊 凱 李 康△

目的探討基于隨機(jī)森林（RF）回歸估計(jì)因果關(guān)系網(wǎng)絡(luò)的效果。方法通過(guò)模擬實(shí)驗(yàn)設(shè)定因果關(guān)系網(wǎng)絡(luò)，對(duì)數(shù)據(jù)標(biāo)準(zhǔn)化后，利用全條件RF回歸對(duì)其進(jìn)行估計(jì)并評(píng)價(jià)其準(zhǔn)確性。另外將該方法用于卵巢癌基因表達(dá)譜數(shù)據(jù)，并對(duì)分析結(jié)果進(jìn)行驗(yàn)證。結(jié)果模擬實(shí)驗(yàn)結(jié)果表明RF回歸對(duì)于預(yù)先設(shè)定網(wǎng)絡(luò)關(guān)系的識(shí)別能力明顯優(yōu)于貝葉斯網(wǎng)絡(luò)方法。當(dāng)選擇合適的閾值時(shí)，隨著樣本含量的增加基于隨機(jī)森林回歸方法構(gòu)建的網(wǎng)絡(luò)準(zhǔn)確性不斷提高，但傳統(tǒng)經(jīng)典的貝葉斯方法效果基本保持不變；實(shí)例分析結(jié)果驗(yàn)證，基于RF回歸方法能夠得到與現(xiàn)有數(shù)據(jù)庫(kù)的網(wǎng)絡(luò)結(jié)構(gòu)。結(jié)論應(yīng)用基于RF回歸方法估計(jì)的網(wǎng)絡(luò)，能夠在樣本量較少的情況下得出準(zhǔn)確度較高的網(wǎng)絡(luò)。

調(diào)控網(wǎng)絡(luò) 隨機(jī)森林回歸 貝葉斯網(wǎng)絡(luò)

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

原理與方法

1．基本思想

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

2．網(wǎng)絡(luò)變量關(guān)聯(lián)計(jì)算

網(wǎng)絡(luò)變量關(guān)聯(lián)可以用通過(guò)RF回歸得到的變量重要性測(cè)量值（VIM）來(lái)衡量。基于樹(shù)方法的變量重要性測(cè)量方法有多種，本文使用方差改變量法，即對(duì)回歸樹(shù)的每一個(gè)節(jié)點(diǎn)t，計(jì)算由變量分裂導(dǎo)致的輸出變量方差的總減少量，定義為

n（t）表示到達(dá)節(jié)點(diǎn)t的樣本量，Var（t）為應(yīng)變量的樣本方差；n（tL）和n（tR）分別表示由節(jié)點(diǎn)t分裂的兩新節(jié)點(diǎn)的樣本量，即有n（t）＝n（tL）＋n（tR）；Var（tL）和Var（tR）為兩新節(jié)點(diǎn)的樣本方差。對(duì)于一棵樹(shù)，一個(gè)變量的重要性可以通過(guò)用這個(gè)變量分裂的所有節(jié)點(diǎn)的I（·）值相加獲得，RF回歸變量Xk的重要性可以用b棵樹(shù)I（·）的平均值進(jìn)行衡量，用公式表示為

3．網(wǎng)絡(luò)連接權(quán)重的計(jì)算

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

規(guī)定上式j(luò)＝k時(shí)VIMkj＝0，即對(duì)角線元素為零。

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

其中S是用來(lái)估計(jì)VIM值的樣本數(shù)據(jù)。因此在數(shù)據(jù)標(biāo)準(zhǔn)化的情況下，不同RF回歸模型的VIM值具有一定的可比性，從而可以對(duì)矩陣VIMkj所有的值排序，再通過(guò)設(shè)定的閾值c（VIM）（或通過(guò)確定連接的邊數(shù)），取排序中｛VIMkj≥c（VIM）｝估計(jì)出網(wǎng)絡(luò)連接的邊和方向。如何給出合適的閾值是一個(gè)需要研究的問(wèn)題。

上述算法實(shí)現(xiàn)可以使用random Forest R軟件包。

模擬研究

1．無(wú)噪聲情況的網(wǎng)絡(luò)估計(jì)

（1）模擬實(shí)驗(yàn)?zāi)康暮蜅l件設(shè)置：研究網(wǎng)絡(luò)變量之間存在線性關(guān)系及非線性相關(guān)關(guān)系時(shí)，基于RF方法構(gòu)建相關(guān)網(wǎng)絡(luò)的準(zhǔn)確性，同時(shí)與經(jīng)典的貝葉斯網(wǎng)絡(luò)構(gòu)建方法進(jìn)行比較。實(shí)驗(yàn)設(shè)置15個(gè)變量和16條有向邊（圖1），選取的網(wǎng)絡(luò)中包括：一個(gè)變量調(diào)控多個(gè)變量，如G12→｛G14，G17｝；多個(gè)變量調(diào)控一個(gè)變量｛G4，G9，G13｝→G8，線性級(jí)聯(lián)G11→G9→G10，前饋環(huán)路｛G2，G3，G4｝。

（2）網(wǎng)絡(luò)參數(shù)及數(shù)據(jù)：變量之間的調(diào)控分別具有線性關(guān)系和交互關(guān)系，即｛G13，G9｝→G8和｛G2，G4｝→G3為交互作用，其他為線性調(diào)控關(guān)系。誤差服從正態(tài)分布。樣本量分別設(shè)為50、100、200、500和1000。

圖1 模擬實(shí)驗(yàn)設(shè)置的網(wǎng)絡(luò)結(jié)構(gòu)

（3）網(wǎng)絡(luò)數(shù)據(jù)模型：基于RF回歸構(gòu)建網(wǎng)絡(luò)，分別設(shè)定VIM值為前15、20和25條邊作為陽(yáng)性邊，剩余的邊作為陰性邊；對(duì)于貝葉斯網(wǎng)絡(luò)，取結(jié)果為1的邊作為陽(yáng)性邊，結(jié)果為0的邊作為陰性邊，分別使用AUC值及預(yù)測(cè)準(zhǔn)確率（PRE）進(jìn)行比較，對(duì)兩種網(wǎng)絡(luò)構(gòu)建方法構(gòu)建的網(wǎng)絡(luò)進(jìn)行評(píng)價(jià)。以上過(guò)程隨機(jī)重復(fù)100次。

（4）模擬試驗(yàn)結(jié)果

表1模擬結(jié)果顯示，當(dāng)選擇隨機(jī)森林的邊數(shù)（閾值）為15時(shí)（接近真實(shí)邊的數(shù)目），在準(zhǔn)確性評(píng)價(jià)中，AUC和PRE兩項(xiàng)評(píng)價(jià)指標(biāo)明顯優(yōu)于其他兩種情況，同時(shí)明顯優(yōu)于貝葉斯方法；但是，當(dāng)選定的邊數(shù)（閾值）明顯大于真實(shí)邊數(shù)時(shí)（25條邊），其準(zhǔn)確性明顯下降。當(dāng)樣本量小于200時(shí)，RF和貝葉斯建模的整體準(zhǔn)確性不夠穩(wěn)定，在大于200時(shí)隨樣本量的增加趨于平緩（圖2）。

2．具有噪聲情況下的網(wǎng)絡(luò)估計(jì)

模擬實(shí)驗(yàn)?zāi)康暮蜅l件設(shè)置：研究基于隨機(jī)森林回歸方法在網(wǎng)絡(luò)變量間存在噪聲變量時(shí)的識(shí)別能力，同時(shí)與不存在噪聲變量的模擬數(shù)據(jù)分析效果進(jìn)行比較。真實(shí)網(wǎng)絡(luò)變量的條件設(shè)置同前，另加入200個(gè)服從正態(tài)分布的噪聲變量。模擬隨機(jī)重復(fù)100次。

模擬結(jié)果見(jiàn)圖3。結(jié)果表明：利用隨機(jī)森林回歸估計(jì)網(wǎng)絡(luò)的方法與不加入噪聲變量的結(jié)果相近，即加入噪聲變量后仍然能夠很好地識(shí)別出變量之間的調(diào)控關(guān)系，并且保持較好的準(zhǔn)確性。

3．閾值選擇

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

表1 基于隨機(jī)森林回歸網(wǎng)絡(luò)構(gòu)建方法與貝葉斯網(wǎng)絡(luò)對(duì)網(wǎng)絡(luò)的識(shí)別效果（±S）

表1 基于隨機(jī)森林回歸網(wǎng)絡(luò)構(gòu)建方法與貝葉斯網(wǎng)絡(luò)對(duì)網(wǎng)絡(luò)的識(shí)別效果（±S）

評(píng)價(jià)指標(biāo)樣本量隨機(jī)森林網(wǎng)絡(luò)構(gòu)建方法15條邊20條邊25條邊貝葉斯網(wǎng)絡(luò)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 基于隨機(jī)森林回歸網(wǎng)絡(luò)構(gòu)建方法與貝葉斯網(wǎng)絡(luò)估計(jì)網(wǎng)絡(luò)的準(zhǔn)確性

圖3 基于隨機(jī)森林回歸的網(wǎng)絡(luò)構(gòu)建方法抗噪能力情況

實(shí)例分析

數(shù)據(jù)來(lái)源：通過(guò)TCGA數(shù)據(jù)庫(kù)下載570例卵巢癌患者以及8例健康對(duì)照數(shù)據(jù)的基因表達(dá)譜數(shù)據(jù)［7］，其中包含12042個(gè)基因的表達(dá)值。

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

進(jìn)而，通過(guò)查詢已有的基因／蛋白互作網(wǎng)絡(luò)數(shù)據(jù)庫(kù)GeneMANIA［8］，發(fā)現(xiàn)這12個(gè)基因中，有9條邊的因果關(guān)系出現(xiàn)在GeneMANIA的網(wǎng)絡(luò)中（圖5）。在這9條邊中，基因CCNB1、CCNB2、CHEK1和CHEK2參與了比較多的調(diào)控，這個(gè)結(jié)果與隨機(jī)森林回歸所構(gòu)建網(wǎng)絡(luò)的結(jié)構(gòu)相一致，說(shuō)明這幾個(gè)基因在p53信號(hào)通路中起到了很重要的作用。在剩余的7條邊中，有4條邊是與基因GTSE1（即PRKAR1A）有關(guān)，而在GeneMANIA中并沒(méi)有相應(yīng)的調(diào)控關(guān)系，提示這4條邊需要進(jìn)一步研究。

圖4 卵巢癌患者p53信號(hào)通路中基因網(wǎng)絡(luò)結(jié)構(gòu)

圖5 12個(gè)基因在GeneMANIA中的網(wǎng)絡(luò)關(guān)系

討 論

1．基于RF回歸估計(jì)網(wǎng)絡(luò)的方法，兩變量連接的權(quán)重wkj主要是利用全自變量條件下標(biāo)準(zhǔn)化的VIMkj得到。實(shí)際中，基于樹(shù)方法的變量重要性測(cè)量方法有多種，本文使用方差改變量法。另外，VIM值的計(jì)算還有隨機(jī)置換法，這種方法與方差改變量法相比在理論上有一定的優(yōu)勢(shì)，但在實(shí)際應(yīng)用中得到的結(jié)果可能相近，而且需要更多的計(jì)算。

2．理論上，RF回歸估計(jì)網(wǎng)絡(luò)對(duì)數(shù)據(jù)的分布和變量之間的關(guān)系不需要做出任何假定，能夠挖掘出具有交互作用的變量和各種復(fù)雜的非線性關(guān)系；而貝葉斯網(wǎng)絡(luò)模型更適合單調(diào)或線性的調(diào)控關(guān)系。貝葉斯模型的優(yōu)點(diǎn)之一是能夠給出調(diào)控關(guān)系的方向，但模擬實(shí)驗(yàn)表明RF回歸估計(jì)網(wǎng)絡(luò)方法同樣能夠給出調(diào)控的方向，其原理是對(duì)任一變量，在全自變量條件下，能夠?qū)ζ浞讲钭龀鲚^多解釋的變量最有可能是調(diào)控其變化的變量。相對(duì)而言，對(duì)于一個(gè)變量調(diào)控多個(gè)變量更容易識(shí)別，如圖1中的G12→｛G17，G14｝；而對(duì)于多個(gè)變量調(diào)控一個(gè)變量則相對(duì)容易判錯(cuò)方向，如｛G4，G9，G13｝→G8。

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

4．關(guān)于RF回歸參數(shù)的設(shè)置，在隨機(jī)森林中，生成每一棵分類樹(shù)時(shí)，所應(yīng)用的自助樣本集從原始的訓(xùn)練樣本集中隨機(jī)選取，每一棵樹(shù)所應(yīng)用的變量也是從所有變量mall中隨機(jī)選取，并用袋外數(shù)據(jù)來(lái)衡量回歸的效果。隨機(jī)森林中最重要的參數(shù)有兩個(gè)，即建立回歸樹(shù)的個(gè)數(shù)ntree和建立每個(gè)節(jié)點(diǎn)所取的變量?jī)蓴?shù)目ntry。理論上講，ntree越大得到的回歸森林越穩(wěn)定，本文取ntree＝1000；ntry的選取主要有和ntry＝mall－1兩種方法，因本文變量數(shù)目不是很大，取ntry＝mall－1。

5．實(shí)例分析結(jié)果表明，網(wǎng)絡(luò)中多數(shù)的連接邊能夠被識(shí)別出來(lái)，并且能夠給出調(diào)控的方向。然而，需要注意的是，實(shí)際各種變量之間的關(guān)系可能有多種，如基因共表達(dá)、基因調(diào)控、代謝通路等，其中共表達(dá)為無(wú)方向性，而RF回歸方法同樣能夠?qū)Υ诉M(jìn)行分析，如VIMij≈VIMji時(shí)其方向不明確。需要注意的是，對(duì)于斷面數(shù)據(jù)得到的因果網(wǎng)絡(luò)關(guān)系最終需要實(shí)驗(yàn)的方法加以證實(shí)。

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

（責(zé)任編輯：郭海強(qiáng)）

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

國(guó)家自然科學(xué)基金（81473072）；中國(guó)博士后面上項(xiàng)目（2015M571445）

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