隨機森林在基因調控網絡推斷上的比較與應用*

張銘智 尤東方 何文靜 張汝陽,2 陳 峰,3,4 胡志斌 趙 楊,6,7△

基因調控網絡(gene regulatory network,GRN)是當前功能基因組學所研究的重要內容之一，作為一種描述基因間相互作用關系的方式，在推斷復雜疾病的致病原理過程中發揮著重要的作用。通過對基因調控網絡的分析，我們能夠更加深入地了解各基因的生物學功能、理解基因間的調控機理并推斷出未知基因的功能，這對疾病診斷、臨床實踐、藥物研發等方面有指導性的意義[1-2]。近年來，隨著高通量生物實驗技術的快速發展和計算機技術的進步，使得從高維生物組學數據中推斷復雜的基因調控網絡成為可能。大量推斷方法應運而生：基于信息論方法[3-5]、基于模型方法[6-8]、基于監督或非監督學習方法[9-10]等等。這些方法主要是從基因表達數據中獲取功能關系信息，用以模擬真實的網絡結構，從局部甚至是整個基因組的意義上揭示基因間可能存在的相互關系[11]。

隨機森林最早是由Leo Breiman和Adele Cutler共同提出的一種集成算法[12]。該算法不但能夠很好地處理組學數據中存在的數據不完整、受噪聲影響、高維低樣本和基因間的非線性調控等問題，而且能夠大幅度縮短運算時間。因此，越來越多的學者傾向于使用隨機森林算法進行網絡構建。如Huynh-Thu等[13]、Marit Ackermann等[14]均證明了隨機森林的有效性和易擴展性。dynGENIE3(dynamical GENIE3)[15]和iRafNet(integrative random forest for gene regulatory network inference)[16]是兩種具有代表性的擴展方法，在已有文獻中未見兩種方法和經典隨機森林三者之間的比較，故本文旨在比較上述三種算法的網絡重建能力，從而為提高推斷性能提供參考性建議。

方 法

1.基于靜態數據的隨機森林算法

經典隨機森林主要是從靜態基因表達數據中推斷調控網絡。我們令X={X1,X2,…,XN}代表穩定狀態下的基因表達值，其中Xk∈P,k=1,2,…,N代表p個基因在第k個狀態下的表達值：假設每個基因xkj的表達受其他所有基因的共同調控，因此第j個基因的表達值可由如下函數表示：

(1)

(1)采用Bootstrap方法從訓練樣本Xj={(xkj,Xk-j),k=1,2,…,N}中有放回地隨機抽取N條觀測得到一個子樣本集合；重復上述抽樣方法B次，得到B個子樣本用于構建B棵獨立的回歸決策樹。

(3)將所有的樹進行匯總，根據樹分裂時節點方差平均減少值(IncNodePurity)的大小來判斷特征對基因j的影響程度：該值越大，說明對應的特征對基因j存在調控的可能性越大。

(4)按步驟(2)和(3)分別計算p個基因與待選特征之間的節點方差平均減少值(IncNodePurity)，將p(p-1)個節點方差平均減少值從大到小排序：排名越靠前，則越有把握認為真實網絡中存在對應的調控關系。最后通過與金標準的比較來評價算法的性能。

2.基于動態數據的隨機森林算法

動態數據包含了時間歷程，顯示了網絡在受到外界干擾時做出的反應以及干擾移除后回歸穩態水平的動態變化過程，理論上更能說明基因調控的方向性和因果關系。因此有學者提出了一種利用時間序列數據進行網絡推斷的方法dynGENIE3。該算法與經典隨機森林算法的區別在于訓練樣本的構造不同：令DTS={X(t1),X(t2),…,X(tN)}代表時間序列數據中各時間點的表達值(TS:time-series)，其中X(tk)∈p,k=1,2,…,N代表p個基因在第k個時間節點的表達值：X(tk)=(X(tk)1,X(tk)2,…,X(tk)p)T。假設基因j在tk時刻的變化率受當前時刻所有基因的影響，即滿足常微分方程：

(2)

動力學參數αj代表基因j的衰減率(decay rate)。當時間序列數據中各個相鄰時間點的間隔一致時:

αj=1/(t(k+1)-tk)

(3)

通過(2)和(3)可得到基因j在t(k+1)時刻表達值的估計:

(4)

(5)

3.結合靜態數據和動態數據的隨機森林算法

靜態數據和動態數據是從不同角度體現網絡的內在結構，同時利用這兩種數據進行網絡推斷能夠更加全面地反映出真實網絡狀態。

(6)

此時訓練樣本可表示為：

(7)

將該訓練樣本與算法2中的訓練樣本合并：

(8)

基于LSj進行網絡推斷即可同時利用兩種數據中蘊含的基因關聯信息。

(9)

(10)

4.算法的評價標準

將得到的結果與金標準進行比較，把推斷結果分為四類：真陽性(true positive,TP)、假陽性(false positive,FP)、真陰性(true negative,TN)以及假陰性(false negative,FN)，如表1所示。則特異度(specificity,Sp)、靈敏度(sensitivity,Se)、查全率(recall,R)和查準率(precision,P)的計算公式如下：

(11)

(12)

(13)

(14)

本文主要通過受試者工作特征曲線線下面積(Area under ROC,AUC)以及精度召回率曲線線下面積(Area under PR,AUPR)的大小來衡量不同算法的優劣性。

以上所有算法的網絡構建及性能評價均由R-3.5.1軟件實現：經典隨機森林使用‘randomForest’包構建；dynGENIE3[15]和iRafNet[16]算法分別使用對應文獻提供的R包構建；AUC和AUPR通過‘ROCR’包計算所得。

數據來源

本文模擬數據來自Dialogue for Reverse Engineering Assessments and Methods第四次競賽(DREAM4)，包含5個10基因規模和5個100基因規模的網絡調控數據，每個網絡均含有靜態數據和動態數據，且所有數據都經過歸一化處理。其中靜態數據包括野生型數據(wild type)、基因敲除數據(knock out)和基因敲低數據(knock down)；野生型數據為穩定狀態下各基因的表達量、敲除數據是抑制某基因時其余基因的表達情況、敲低數據則是在某個基因表達強度減半時其余基因的表達情況。動態數據為時間序列數據(time series)，反映了穩態網絡在受到外界干擾及干擾去除后所有基因表達量的變化情況，測量的時間間隔為50秒，共測量21次。同時還提供了這10個網絡中真實存在的調控關系(金標準)，用于評價算法的優劣性。所有的模擬數據和金標準可在DREAM4官網上下載(http://dreamchallenges.org/project-list/dream4-2009/)。驗證數據為DREAM5提供的大腸桿菌基因表達微陣列數據，包括805個微陣列芯片，4511個基因的表達值。基因表達數據及每個微陣列芯片的詳細描述均可在DREAM5官網中獲取(https://www.synapse.org/#!Synapse:syn2787209/wiki/70349)。

結果與分析

從DREAM4參賽隊伍的結果中發現，使用基因敲除數據推斷網絡的結果普遍優于從其他靜態數據中推斷的結果，這可能是因為基因敲除數據中蘊含的信息更加豐富[17]。因此，對于靜態數據我們主要從基因敲除數據來推斷網絡。

表2 比較各種方法在10基因規模網絡中的AUC和AUPR值

TS：時間序列數據；KO：基因敲除數據

表3 比較各種方法在100基因規模網絡中的AUC和AUPR值

TS：時間序列數據；KO：基因敲除數據

由表2可知：當網絡結構較為簡單時，dynGENIE3和iRafNet算法在5個網絡中的ROC和PR曲線線下面積相比于隨機森林均得到了明顯的提升。同時發現dynGENIE3:TS和dynGENIE3:TS+KO曲線下面積基本一致，說明在該情況下，基因敲除數據相較于時間序列數據包含較少的額外信息。

由表3可知：當網絡結構變復雜時，dynGENIE3和iRafNet算法在各個網絡中都能提高PR曲線下面積。比較dynGENIE3:TS和dynGENIE3:TS+KO可發現當結合基因敲除數據進行網絡推斷時ROC曲線下面積明顯增加：5個網絡分別增加了19.40%、29.59%、27.13%、26.97%和15.66%(以ROC曲線下面積等于0.5為基線)，且后者的PR曲線下面積也分別增加了0.0579、0.0388、0.0696、0.0523和0.0286。進一步比較dynGENIE3:TS+KO和iRafNet：對于ROC曲線下面積除了網絡5，iRafNet均優于dynGENIE3:TS+KO；對于PR曲線下面積iRafNet均明顯大于dynGENIE3:TS+KO。

為了進一步驗證以上結果，我們分別用這三種方法對大腸桿菌的基因表達數據構建網絡并計算AUC和AUPR值。結果如下：RF對應的AUC和AUPR值分別為0.5633和0.030；dynGENIE3:TS為0.5579和0.019；dynGENIE3:TS+KO為0.5837和0.031；iRafNet為0.6482和0.102。無論是AUC還是AUPR值，iRafNet算法均優于dynGENIE3和經典隨機森林算法；并且結合基因敲除數據的dynGENIE3算法能夠有效提高網絡預測的準確性。與模擬數據的結果基本一致。

討 論

基因調控網絡的構建已經成為當前研究的熱點領域之一，是對不同組學數據進行處理和挖掘的過程，將表達數據轉換成由若干節點和邊組成的圖形，為進一步深入了解節點間復雜的調控關系和作用機理提供了參考模型。本文基于模擬數據和大腸桿菌基因表達數據對經典隨機森林(RF)、dynGENIE3和iRafNet三種方法的網絡推斷能力進行了比較，這三種方法均將含有p個基因的網絡推斷問題轉化為p個不同的特征選擇問題。從100基因大小網絡的推斷中我們發現：RF算法對應的ROC曲線下面積大于dynGENIE3:TS，與前文提到的利用動態數據推斷網絡更加有效這一說法相矛盾。主要原因是時間序列數據中只包含21個時間點的表達值，不能夠充分地反映各基因間存在的調控關系。而dynGENIE3:TS對應的PR曲線下面積更大，說明其推斷結果的覆蓋率和準確率比經典隨機森林更高。故當時間序列數據提供更多時間點的基因表達測量值時，理論上其推斷性能將優于靜態數據。對于dynGENIE3算法，在結合靜態數據后ROC和PR曲線下面積均得到了提升，這表明同時結合多種數據進行網絡推斷確實能夠提高算法的性能。雖然iRafNet和dynGENIE3:TS+KO都是結合動靜態數據網絡推斷算法，但兩者的性能卻存在著差異：dynGENIE3直接將基因敲除數據作為訓練樣本的一部分，每棵樹都利用部分信息；而RafNet則是將基因敲除數據轉換成抽樣權重，使得每次分裂時都能充分利用數據中蘊含的信息。因此，如何充分提取數據中的信息將是提高算法性能的另一重要途徑。

綜上所述，隨機森林是一種高效便捷的機器學習方法，不但能識別出基因間存在的非線性關系，也能同時整合多種不同類型的數據，在未來的研究中還可結合單核苷酸多態性(single nucleotide polymorphisms，SNP)、甲基化表達等數據進一步完善基因調控網絡。隨機森林本身也存在著不足，該算法無法很好地分離出網絡中的直接調控效應和間接調控效應：當真實網絡中存在i→k→j通路時，很可能會錯誤地認為基因j是受基因i直接調控的：i→j，使推斷結果中假陽性比率上升，整體精確率下降；而在i→k→j這一通路的基礎上如果同時存在i→j，又會降低i→j、k→j的推斷效能，使查全率下降[18]，共同導致PR曲線下面積的減少。另一方面，如果基因i同時對基因j和k有著調控關系，即j←i→k，將會導致基因j和k之間產生虛假關聯：j←k或j→k。因此，如何將直接效應和間接效應分離，如何對算法產生的結果進行假設檢驗并剔除虛假關聯和混雜將是提高算法性能的重要研究方向。