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

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

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

原理與方法

1.基本思想

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

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

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

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

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

DiffRank-RF計算局部結構改變測量指標ΔC和全局結構改變測量指標ΔBC的公式分別為

圖1 網絡示例圖

(1)

(2)

(3)

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

模擬研究

1.模擬實驗目的：通過模擬實驗評價DiffRank-RF算法在不同樣本量情況下篩選差異位點的準確性和穩(wěn)定性，同時與SAM、PLS方法進行比較，探討DiffRank-RF算法最優(yōu)的適用范圍和λ參數(shù)設置。

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

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

圖2 有向網絡模擬實驗設置條件

4.模擬實驗結果：表1模擬實驗結果顯示，在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值都越高，當樣本量大于200時效果趨于平緩，PRE指標在DiffRank-RF方法λ=0.75時要優(yōu)于其他情況(圖3B)。

表1 DiffRank-RF差異網絡分析與SAM、PLS比較結果

圖3 DiffRank-RF差異網絡分析與SAM、PLS的準確性比較

實例分析

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

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

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

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

圖4 DiffRank-RF分析結果在GeneMINIA中的關系示意圖

討 論

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

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

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

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