蔡雨晴 李軼群 徐 歡 宋 微 楊 凱 王文杰 李 康△

1.哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計學(xué)教研室(150081)2.哈爾濱工業(yè)大學(xué)生命科學(xué)與技術(shù)學(xué)院

差異表達(dá)分析常被用于各種疾病標(biāo)志物的篩選研究中，如傳統(tǒng)的t檢驗、顯著性分析(significance of microarrays,SAM)檢驗、偏最小二乘(pa least square,PLS)等方法。然而，這些方法主要是通過比較不同分組之間基因表達(dá)均值的差異篩選標(biāo)記物，忽視了物質(zhì)之間的相互調(diào)控關(guān)系，致使研究結(jié)果不夠穩(wěn)定或檢驗效率低。在組學(xué)研究中，由于基因調(diào)控和蛋白質(zhì)的互相作用，很有可能在表達(dá)量上還沒有呈現(xiàn)出明顯差別時，在調(diào)控關(guān)系上已經(jīng)發(fā)生了一定的改變。差異網(wǎng)絡(luò)分析方法更加注重不同分組情況下調(diào)控關(guān)系和網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的差別，并由此篩選出具有潛在生物學(xué)意義的標(biāo)記物。本文提出DiffRank-RF差異網(wǎng)絡(luò)分析方法，通過模擬實(shí)驗評價該方法的準(zhǔn)確性和適用條件，并與傳統(tǒng)的變量篩選方法進(jìn)行比較，最后應(yīng)用于乳腺癌實(shí)際數(shù)據(jù)，得到相應(yīng)的分析結(jié)果。

原理與方法

1.基本思想

隨機(jī)森林方法提高了預(yù)測精度，對多重共線性不敏感。利用隨機(jī)森林(random forest,RF)回歸模型，可以建立任一變量Xk對其它變量的回歸模型：

Xk=RF(X1,X2,…,Xk-1,Xk+1,…,Xm)+ε

網(wǎng)絡(luò)共有m個變量，其中ε為模型的殘差。根據(jù)衡量變量重要的VIM值作為有向連接兩節(jié)點(diǎn)的權(quán)重，可以建立RF網(wǎng)絡(luò)[1]。利用R包randomForest即可實(shí)現(xiàn)通過隨機(jī)森林回歸構(gòu)建網(wǎng)絡(luò)。

差異網(wǎng)絡(luò)分析使用DiffRank[2]算法。首先根據(jù)隨機(jī)森林(RF)構(gòu)建網(wǎng)絡(luò)，再結(jié)合網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的局部指標(biāo)連接權(quán)重(connectivity)、度(degree)以及全局指標(biāo)最短路徑(shortest path)等統(tǒng)計量發(fā)現(xiàn)導(dǎo)致網(wǎng)絡(luò)差異的重要變量。連接權(quán)重即變量之間的關(guān)聯(lián)強(qiáng)弱，可用RF建網(wǎng)得到的VIM值度量，并用連接邊線的粗細(xì)表示權(quán)重大小(見圖1)。度是在網(wǎng)絡(luò)中某一變量的連接邊數(shù)量，圖1中可見變量G1的度為5。DiffRank-RF算法將被分析節(jié)點(diǎn)的所有直接連接點(diǎn)的權(quán)重進(jìn)行相加得到網(wǎng)絡(luò)局部測量指標(biāo)ΔC。最短路徑是指變量間權(quán)重之和最小的一條連接路徑，DiffRank-RF計算經(jīng)過節(jié)點(diǎn)的最短路徑數(shù)量占所有最短路徑數(shù)量的比值來表示節(jié)點(diǎn)的中介中心性(between centrality，BC)，可以分析網(wǎng)絡(luò)中所有節(jié)點(diǎn)(包括直接連接點(diǎn)和間接連接點(diǎn))對被分析節(jié)點(diǎn)的影響。當(dāng)節(jié)點(diǎn)的度或連接權(quán)重較小，卻經(jīng)過網(wǎng)絡(luò)的多數(shù)最短路徑時，仍可認(rèn)為該節(jié)點(diǎn)是網(wǎng)絡(luò)中的重要節(jié)點(diǎn)，ΔBC值能夠反映這一現(xiàn)象。

2.統(tǒng)計量計算

DiffRank-RF計算局部結(jié)構(gòu)改變測量指標(biāo)ΔC和全局結(jié)構(gòu)改變測量指標(biāo)ΔBC的公式分別為

圖1 網(wǎng)絡(luò)示例圖

(1)

(2)

(3)

其中，A和B分別代表兩個不同分組情況下隨機(jī)森林回歸所構(gòu)建的網(wǎng)絡(luò)，分別包含N個變量。VIM是隨機(jī)森林得到的變量重要性評分，表示變量v與其它相連變量的連接權(quán)重。πvi為變量v在網(wǎng)絡(luò)中第i次迭代的差異評分，用參數(shù)λ結(jié)合兩部分指標(biāo)，λ取值范圍為[0，1]，可根據(jù)模擬試驗選取不同情況下合適的λ值。任一變量的π初始值可設(shè)為1/N，結(jié)果收斂時循環(huán)停止。SPv(s,t)可表示為通過變量v的一個N×N矩陣，在網(wǎng)絡(luò)中任意兩變量s、t的最短路徑若通過變量v，則在矩陣中用1表示，否則用0表示。ΔBC(v)計算通過變量v的最短路徑數(shù)量來反映變量v在網(wǎng)絡(luò)中的中介中心性?；诿恳蛔兞康牟町愒u分π給所有變量排序，π越大表示在差異網(wǎng)絡(luò)中貢獻(xiàn)最大，即所篩選的差異位點(diǎn)。

模擬研究

1.模擬實(shí)驗?zāi)康模和ㄟ^模擬實(shí)驗評價DiffRank-RF算法在不同樣本量情況下篩選差異位點(diǎn)的準(zhǔn)確性和穩(wěn)定性，同時與SAM、PLS方法進(jìn)行比較，探討DiffRank-RF算法最優(yōu)的適用范圍和λ參數(shù)設(shè)置。

2.模擬實(shí)驗設(shè)置：有向模擬網(wǎng)絡(luò)設(shè)置20個變量和25條有向邊(見圖2)，包括變量間的線性調(diào)控和非線性調(diào)控關(guān)系和交互作用，其中線性關(guān)系由線性方程產(chǎn)生，相關(guān)系數(shù)為隨機(jī)產(chǎn)生的固定值，誤差從正態(tài)分布中隨機(jī)抽樣，非線性關(guān)系在線性基礎(chǔ)上指數(shù)形式產(chǎn)生。實(shí)驗設(shè)置樣本量分別為50，100，200，500和1000。對樣本數(shù)據(jù)應(yīng)用隨機(jī)森林回歸方法構(gòu)建兩個網(wǎng)絡(luò)，通過DiffRank-RF進(jìn)行差異網(wǎng)絡(luò)分析，分別使用AUC值及預(yù)測準(zhǔn)確率(PRE)指標(biāo)與SAM和PLS方法進(jìn)行比較。以上過程隨機(jī)重復(fù)100次。

3.閾值選擇：隨機(jī)森林構(gòu)建網(wǎng)絡(luò)時，VIM值通過置換檢驗可以得到其均值的隨機(jī)分布，選取95%分位數(shù)為閾值以判斷節(jié)點(diǎn)之間是否存在真實(shí)邊。在進(jìn)行預(yù)測準(zhǔn)確率比較時，選取PLS結(jié)果中VIP、SAM得分、DiffRank-RF結(jié)果秩次排在前5位的變量為預(yù)測差異變量。

圖2 有向網(wǎng)絡(luò)模擬實(shí)驗設(shè)置條件

4.模擬實(shí)驗結(jié)果：表1模擬實(shí)驗結(jié)果顯示，在AUC評價中，DiffRank-RF方法在λ=0.5時隨樣本量增加AUC值增加最明顯，但穩(wěn)定性較差(見圖3A)，λ=1時穩(wěn)定性最優(yōu)，綜合看來λ=0.75效果最好，且DiffRank-RF不管λ取何值時，效果都優(yōu)于SAM和PLS方法。隨樣本量逐漸增加，DiffRank-RF、SAM和PLS方法AUC值都越高，當(dāng)樣本量大于200時效果趨于平緩，PRE指標(biāo)在DiffRank-RF方法λ=0.75時要優(yōu)于其他情況(圖3B)。

表1 DiffRank-RF差異網(wǎng)絡(luò)分析與SAM、PLS比較結(jié)果

圖3 DiffRank-RF差異網(wǎng)絡(luò)分析與SAM、PLS的準(zhǔn)確性比較

實(shí)例分析

數(shù)據(jù)來源：TCGA數(shù)據(jù)庫中531例乳腺癌患者及63例對照的mRNA基因表達(dá)數(shù)據(jù)，選取p53信號通路進(jìn)行分析。分別選取λ=0、0.75和1，對這條通路內(nèi)所有基因進(jìn)行DiffRank-RF差異網(wǎng)絡(luò)分析，分析結(jié)果見表2。

結(jié)果顯示，DiffRank-RF差異網(wǎng)絡(luò)分析方法λ取0和0.75時篩選的變量有較大重疊，而與λ=1時相比差別較大；同時可以看到DiffRank-RF方法篩選的變量與傳統(tǒng)的SAM和PLS相比差別較大，幾乎無重疊。SAM和PLS兩種方法之間篩選出的結(jié)果則十分相近。

表2 乳腺癌與對照數(shù)據(jù)使用三種方法篩選變量的結(jié)果(排序前10)

通過文獻(xiàn)查閱，CDK4是細(xì)胞周期中G1-S期調(diào)控的中心基因，已發(fā)現(xiàn)CDK4的高表達(dá)廣泛存在于人類的多種腫瘤中，CDK4的異常表達(dá)與腫瘤的發(fā)生密切相關(guān)。CDK4、CDKN2A(p16)和CDK2同屬于CDK家族與細(xì)胞周期調(diào)控有關(guān)的基因，其中CDKN2A是CDK4的抑制因子，阻止細(xì)胞進(jìn)入S期，同時對CDK2也有抑制作用[3]，有研究表明CDKN2A改變會影響乳腺癌患者的生存和預(yù)后[4]。PTEN是繼p53后另一個較為廣泛地與腫瘤發(fā)生關(guān)系密切的基因，對細(xì)胞周期進(jìn)展和細(xì)胞凋亡有重要作用，同時，PTEN與CDK2抑制劑(CDKN1A)對卵巢癌細(xì)胞生長抑制具有協(xié)同作用[5]。在細(xì)胞凋亡的調(diào)控過程中，CASP3和CASP8發(fā)揮了關(guān)鍵作用，其中CASP3的高表達(dá)與乳腺癌生存時間有顯著性關(guān)系[6]。使用GeneMINIA[7]基因/蛋白互作網(wǎng)絡(luò)數(shù)據(jù)庫可以將篩選出的基因畫出網(wǎng)絡(luò)圖，圖4給出了DiffRank-RF方法在λ=0.75時的網(wǎng)絡(luò)示意圖。

圖4 DiffRank-RF分析結(jié)果在GeneMINIA中的關(guān)系示意圖

討 論

傳統(tǒng)的差異基因篩選方法主要是根據(jù)基因表達(dá)量在不同分組中的差異進(jìn)行篩選。實(shí)際的基因網(wǎng)絡(luò)有可能其表達(dá)量改變不大，但其調(diào)控關(guān)系發(fā)生變化，此時傳統(tǒng)方法有較低的檢驗效率，本文給出的DiffRank-RF方法則能夠充分反映不同組間調(diào)控網(wǎng)絡(luò)的差異，篩選出重要的基因。

已有的多種網(wǎng)絡(luò)構(gòu)建方法中，隨機(jī)森林方法能夠識別變量之間的非線性關(guān)系和交互作用，且隨機(jī)森林可以構(gòu)建有向網(wǎng)絡(luò)。由于基因之間的調(diào)控通常為有向的，因此DiffRank-RF方法具有明顯的優(yōu)勢。

DiffRank-RF算法根據(jù)λ不同取值能夠發(fā)現(xiàn)網(wǎng)絡(luò)中不同功能的基因，當(dāng)λ=1時，基因排序靠前，表明該基因與直接關(guān)聯(lián)基因的調(diào)控關(guān)系較強(qiáng)或直接關(guān)聯(lián)基因數(shù)量較多，即在網(wǎng)絡(luò)局部作用較大；當(dāng)λ=0時，基因排序靠前，表明其在網(wǎng)絡(luò)中的中介中心性較高，可被視為網(wǎng)絡(luò)的中心基因，參與網(wǎng)絡(luò)的全局調(diào)控。需要注意：當(dāng)變量數(shù)目較少時，網(wǎng)絡(luò)中的最短路徑數(shù)量也會相對減少，此時全局指標(biāo)(最短路徑算法)應(yīng)用有限，應(yīng)更多的利用連接權(quán)重進(jìn)行差異網(wǎng)絡(luò)分析，λ可適當(dāng)取較大的值；而當(dāng)變量數(shù)目較多時，結(jié)合全局指標(biāo)能夠納入更多生物學(xué)信息，此時建議λ取值0.75。

本文在篩選變量時，主要根據(jù)評價統(tǒng)計量值的大小排序選擇最前面的基因。為了能夠?qū)ζ溥M(jìn)行檢驗，可以使用置換檢驗的方法，根據(jù)檢驗的P值進(jìn)行篩選。