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

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

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

方 法

1.基于靜態(tài)數(shù)據(jù)的隨機(jī)森林算法

經(jīng)典隨機(jī)森林主要是從靜態(tài)基因表達(dá)數(shù)據(jù)中推斷調(diào)控網(wǎng)絡(luò)。我們令X={X1,X2,…,XN}代表穩(wěn)定狀態(tài)下的基因表達(dá)值，其中Xk∈P,k=1,2,…,N代表p個(gè)基因在第k個(gè)狀態(tài)下的表達(dá)值：假設(shè)每個(gè)基因xkj的表達(dá)受其他所有基因的共同調(diào)控，因此第j個(gè)基因的表達(dá)值可由如下函數(shù)表示：

(1)

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

(3)將所有的樹進(jìn)行匯總，根據(jù)樹分裂時(shí)節(jié)點(diǎn)方差平均減少值(IncNodePurity)的大小來判斷特征對(duì)基因j的影響程度：該值越大，說明對(duì)應(yīng)的特征對(duì)基因j存在調(diào)控的可能性越大。

(4)按步驟(2)和(3)分別計(jì)算p個(gè)基因與待選特征之間的節(jié)點(diǎn)方差平均減少值(IncNodePurity)，將p(p-1)個(gè)節(jié)點(diǎn)方差平均減少值從大到小排序：排名越靠前，則越有把握認(rèn)為真實(shí)網(wǎng)絡(luò)中存在對(duì)應(yīng)的調(diào)控關(guān)系。最后通過與金標(biāo)準(zhǔn)的比較來評(píng)價(jià)算法的性能。

2.基于動(dòng)態(tài)數(shù)據(jù)的隨機(jī)森林算法

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

(2)

動(dòng)力學(xué)參數(shù)αj代表基因j的衰減率(decay rate)。當(dāng)時(shí)間序列數(shù)據(jù)中各個(gè)相鄰時(shí)間點(diǎn)的間隔一致時(shí):

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

(3)

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

(4)

(5)

3.結(jié)合靜態(tài)數(shù)據(jù)和動(dòng)態(tài)數(shù)據(jù)的隨機(jī)森林算法

靜態(tài)數(shù)據(jù)和動(dòng)態(tài)數(shù)據(jù)是從不同角度體現(xiàn)網(wǎng)絡(luò)的內(nèi)在結(jié)構(gòu)，同時(shí)利用這兩種數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)推斷能夠更加全面地反映出真實(shí)網(wǎng)絡(luò)狀態(tài)。

(6)

此時(shí)訓(xùn)練樣本可表示為：

(7)

將該訓(xùn)練樣本與算法2中的訓(xùn)練樣本合并：

(8)

基于LSj進(jìn)行網(wǎng)絡(luò)推斷即可同時(shí)利用兩種數(shù)據(jù)中蘊(yùn)含的基因關(guān)聯(lián)信息。

(9)

(10)

4.算法的評(píng)價(jià)標(biāo)準(zhǔn)

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

(11)

(12)

(13)

(14)

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

以上所有算法的網(wǎng)絡(luò)構(gòu)建及性能評(píng)價(jià)均由R-3.5.1軟件實(shí)現(xiàn)：經(jīng)典隨機(jī)森林使用‘randomForest’包構(gòu)建；dynGENIE3[15]和iRafNet[16]算法分別使用對(duì)應(yīng)文獻(xiàn)提供的R包構(gòu)建；AUC和AUPR通過‘ROCR’包計(jì)算所得。

數(shù)據(jù)來源

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

結(jié)果與分析

從DREAM4參賽隊(duì)伍的結(jié)果中發(fā)現(xiàn)，使用基因敲除數(shù)據(jù)推斷網(wǎng)絡(luò)的結(jié)果普遍優(yōu)于從其他靜態(tài)數(shù)據(jù)中推斷的結(jié)果，這可能是因?yàn)榛蚯贸龜?shù)據(jù)中蘊(yùn)含的信息更加豐富[17]。因此，對(duì)于靜態(tài)數(shù)據(jù)我們主要從基因敲除數(shù)據(jù)來推斷網(wǎng)絡(luò)。

表2 比較各種方法在10基因規(guī)模網(wǎng)絡(luò)中的AUC和AUPR值

TS：時(shí)間序列數(shù)據(jù)；KO：基因敲除數(shù)據(jù)

表3 比較各種方法在100基因規(guī)模網(wǎng)絡(luò)中的AUC和AUPR值

TS：時(shí)間序列數(shù)據(jù)；KO：基因敲除數(shù)據(jù)

由表2可知：當(dāng)網(wǎng)絡(luò)結(jié)構(gòu)較為簡(jiǎn)單時(shí)，dynGENIE3和iRafNet算法在5個(gè)網(wǎng)絡(luò)中的ROC和PR曲線線下面積相比于隨機(jī)森林均得到了明顯的提升。同時(shí)發(fā)現(xiàn)dynGENIE3:TS和dynGENIE3:TS+KO曲線下面積基本一致，說明在該情況下，基因敲除數(shù)據(jù)相較于時(shí)間序列數(shù)據(jù)包含較少的額外信息。

由表3可知：當(dāng)網(wǎng)絡(luò)結(jié)構(gòu)變復(fù)雜時(shí)，dynGENIE3和iRafNet算法在各個(gè)網(wǎng)絡(luò)中都能提高PR曲線下面積。比較dynGENIE3:TS和dynGENIE3:TS+KO可發(fā)現(xiàn)當(dāng)結(jié)合基因敲除數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)推斷時(shí)ROC曲線下面積明顯增加：5個(gè)網(wǎng)絡(luò)分別增加了19.40%、29.59%、27.13%、26.97%和15.66%(以ROC曲線下面積等于0.5為基線)，且后者的PR曲線下面積也分別增加了0.0579、0.0388、0.0696、0.0523和0.0286。進(jìn)一步比較dynGENIE3:TS+KO和iRafNet：對(duì)于ROC曲線下面積除了網(wǎng)絡(luò)5，iRafNet均優(yōu)于dynGENIE3:TS+KO；對(duì)于PR曲線下面積iRafNet均明顯大于dynGENIE3:TS+KO。

為了進(jìn)一步驗(yàn)證以上結(jié)果，我們分別用這三種方法對(duì)大腸桿菌的基因表達(dá)數(shù)據(jù)構(gòu)建網(wǎng)絡(luò)并計(jì)算AUC和AUPR值。結(jié)果如下：RF對(duì)應(yīng)的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算法均優(yōu)于dynGENIE3和經(jīng)典隨機(jī)森林算法；并且結(jié)合基因敲除數(shù)據(jù)的dynGENIE3算法能夠有效提高網(wǎng)絡(luò)預(yù)測(cè)的準(zhǔn)確性。與模擬數(shù)據(jù)的結(jié)果基本一致。

討 論

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

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