哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計學(xué)教研室(150081) 宋 微 王文杰 徐 歡 蔡雨晴 李 康

【提 要】 目的 探討基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法的性能，并將其應(yīng)用于卵巢癌基因表達譜數(shù)據(jù)分析。方法 通過模擬實驗評價其識別差異邊及差異節(jié)點的準(zhǔn)確性，并與傳統(tǒng)方法做對比。同時應(yīng)用上皮性卵巢癌基因組學(xué)數(shù)據(jù)，構(gòu)建差異網(wǎng)絡(luò)模型。結(jié)果 模擬試驗結(jié)果表明，基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法識別差異邊能力明顯優(yōu)于高斯圖模型方法；實例分析結(jié)果表明，本文方法構(gòu)建的差異網(wǎng)絡(luò)模型具有實際意義。結(jié)論 應(yīng)用基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法能得出準(zhǔn)確度較高的差異網(wǎng)絡(luò)，效果優(yōu)于傳統(tǒng)方法。

基因的作用通常不是獨立的，而是基因之間相互影響、相互制約，形成復(fù)雜的基因調(diào)控網(wǎng)絡(luò)。某些疾病發(fā)生時，基因表達量變化不明顯，而基因網(wǎng)絡(luò)中調(diào)控關(guān)系卻發(fā)生明顯改變，這種情況傳統(tǒng)的差異分析方法通常無法識別。而在實驗條件下，由于獨立誤差、人為因素、試驗條件限制等原因，有些差異調(diào)控關(guān)系往往不能被研究者發(fā)現(xiàn)。然而，在大數(shù)據(jù)網(wǎng)絡(luò)結(jié)構(gòu)分析中便可被發(fā)現(xiàn)具有統(tǒng)計學(xué)顯著性。因此，在疾病研究領(lǐng)域，差異網(wǎng)絡(luò)的研究已越來越受到研究者的重視。目前最常用的基因調(diào)控網(wǎng)絡(luò)模型有布爾網(wǎng)絡(luò)模型、線性組合模型、加權(quán)矩陣模型、高斯圖模型和貝葉斯網(wǎng)絡(luò)模型等[1-2]。其中圖模型是在變量條件獨立的情況下構(gòu)建基因網(wǎng)絡(luò)，排除了其他變量的影響，因此能在復(fù)雜的相互關(guān)系中準(zhǔn)確提取出兩變量間的真實調(diào)控關(guān)系。貝葉斯因子(Bayes factor,BF)是基于貝葉斯準(zhǔn)則研究得出的方法，廣泛用于模型比較，通過計算BF可得出差異模型的直接證據(jù)。本文將高斯圖模型與貝葉斯準(zhǔn)則相結(jié)合，建立差異網(wǎng)絡(luò)模型，推斷出網(wǎng)絡(luò)中差異的調(diào)控關(guān)系，并與一般的高斯圖模型方法相比較，考察其有效性。最后運用本文給出的方法對卵巢癌基因表達譜數(shù)據(jù)進行分析，做出生物學(xué)解釋。

原理與方法

1.高斯圖模型

(1)

基于以上高斯理論，當(dāng)ρXiXj|X-i,-j≠0時，變量i與變量j之間的邊存在于高斯圖模型中。

在處理高維組學(xué)數(shù)據(jù)時，往往會遇到“m?n”的問題，即變量數(shù)m遠大于樣本數(shù)n，此時協(xié)方差矩陣不是唯一的，因此傳統(tǒng)的統(tǒng)計模型無法求出偏相關(guān)系數(shù)。針對此種情況，有學(xué)者提出應(yīng)用調(diào)整回歸模型推導(dǎo)兩變量的偏相關(guān)系數(shù)，將變量Xi作為因變量，剩余變量Xj(j≠i)作為自變量，構(gòu)建回歸方程：

(2)

(3)

(4)

ρXiXj|X-i,-j=

(5)

2.差異網(wǎng)絡(luò)分析

為推斷不同條件下網(wǎng)絡(luò)中的差異結(jié)構(gòu)，將含有m個變量的數(shù)據(jù)按K個分類標(biāo)簽分成不同數(shù)據(jù)集。對于公式(2)，記Y=Xi。根據(jù)極大似然和平均場近似估計理論，可以直接得出K個獨立網(wǎng)絡(luò)結(jié)構(gòu)滿足：

(6)

對于給定的模型M(分類)，基于模型參數(shù)的變分推斷方法，結(jié)合變分下限和可逆跳躍馬爾科夫鏈蒙特卡洛方法(RJMCMC)[4]，可以求得模型證據(jù)p(X|M)的估計值：

(7)

公式中X表示模型中所有變量的觀察值，θ代表所有模型參數(shù)的集合。對于K=1和K=2兩個模型，通過計算貝葉斯因子

(8)

判斷支持哪個模型成立的證據(jù)更為充分，即各變量(如基因)對其他變量調(diào)控關(guān)系在兩類中是否有顯著的差異。對每個變量X1,……，Xm計算貝葉斯因子，如果大于給定閾值，則進一步根據(jù)偏相關(guān)系數(shù)判斷差異邊，最后根據(jù)給定的假發(fā)現(xiàn)率FDR值判定網(wǎng)絡(luò)中的調(diào)控邊，得出差異網(wǎng)絡(luò)。

模擬實驗

研究通過對設(shè)置的網(wǎng)絡(luò)結(jié)構(gòu)偏相關(guān)系數(shù)矩陣的運算產(chǎn)生不同樣本量的模擬數(shù)據(jù)[2]。檢驗基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法(記為Bayes)，將貝葉斯因子選擇的差異網(wǎng)絡(luò)模型與設(shè)置的真實差異網(wǎng)絡(luò)結(jié)構(gòu)進行比較，若分析結(jié)果能夠獲得更多的已知差異網(wǎng)絡(luò)中存在的邊，則說明該方法有更好的分析效果。根據(jù)前期研究結(jié)果和大量文獻經(jīng)驗結(jié)果，本模擬取BF(i)>2和FDR<0.025。

模擬產(chǎn)生兩組數(shù)據(jù)，網(wǎng)絡(luò)中共包含50個節(jié)點，其中有6個為差異節(jié)點(G30、G38、G1、G37、G4、G42)，共有49條共同邊和8條差異邊(黑色)，49條共同邊偏相關(guān)系數(shù)在(-1,1)隨機抽取，8條差異邊的偏相關(guān)系數(shù)在(0.2,0.5)隨機抽取。黑色邊上的數(shù)字“1”和“2”分別表示兩種不同情況下網(wǎng)絡(luò)中的差異邊。數(shù)字“1”為僅在第一組數(shù)據(jù)中有調(diào)控關(guān)系，數(shù)字“2”為僅在第二組數(shù)據(jù)中有調(diào)控關(guān)系(參見圖1和表1)。模擬選取樣本量n1=n2={30,60,90,120,150}，模擬實驗重復(fù)100次。

圖1 模擬實驗1設(shè)置的真實差異網(wǎng)絡(luò)關(guān)系圖

差異邊節(jié)點1節(jié)點2第一組第二組1G30G1102G30G4103G30G37104G30G42105G38G1016G38G4017G38G37018G38G4201

對數(shù)據(jù)分別使用Bayes、GeneNet、FastGGM三種方法進行分析，其中后兩種方法是兩種常用的高斯圖模型方法[2,5]。模型的評價指標(biāo)選用準(zhǔn)確度、假發(fā)現(xiàn)率、靈敏度和特異度，其中準(zhǔn)確度為真陽性邊占檢查出的陽性邊的比例，相當(dāng)于診斷試驗中的陽性預(yù)測值，假發(fā)現(xiàn)率為假陽性邊占檢查出的陽性邊的比例。Bayes方法模擬實驗結(jié)果顯示在圖2中，可以看到：每組樣本量為60時，準(zhǔn)確度維持在0.6左右；在每組樣本量為90時，差異節(jié)點的靈敏度達到1，差異邊的靈敏度在0.75附近。圖3顯示了三種不同方法在不同樣本量下的準(zhǔn)確度和假發(fā)現(xiàn)率，可以看到，Bayes方法明顯優(yōu)于另外兩種常用的方法。

圖2 模擬實驗各評價指標(biāo)隨樣本量不同的變化情況

圖3 三種方法在不同樣本量下準(zhǔn)確度和假發(fā)現(xiàn)率箱式圖

實例分析

本研究通過對卵巢癌患者復(fù)發(fā)情況及基因表達數(shù)據(jù)進行分析，應(yīng)用基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法構(gòu)建差異網(wǎng)絡(luò)模型，得出基因間差異調(diào)控關(guān)系，及存在差異調(diào)控關(guān)系的基因。通過結(jié)合生物學(xué)知識、通路數(shù)據(jù)庫、文獻查詢，對差異網(wǎng)絡(luò)結(jié)果進行生物學(xué)解釋，從基因組學(xué)角度為卵巢癌復(fù)發(fā)機制提供線索。

本研究從R包curatedOvarianData (version:1.8.0)中下載Ⅱ-Ⅳ期上皮性卵巢癌患者基因表達譜數(shù)據(jù)(GSE49997)，排除了臨床信息中缺失復(fù)發(fā)時間以及復(fù)發(fā)時間<90天的患者，基因表達數(shù)據(jù)應(yīng)用Z-score法進行標(biāo)準(zhǔn)化。將患者根據(jù)生存情況分為復(fù)發(fā)組124例和非復(fù)發(fā)組70例，全基因組表達譜數(shù)據(jù)一共測得16024個基因的表達值。由于基因數(shù)目過多，需要先篩選出與上皮性卵巢癌復(fù)發(fā)相關(guān)的基因通路，再對通路中的全部基因構(gòu)建差異網(wǎng)絡(luò)模型。變量篩選有助于提高建模效率，使得差異網(wǎng)絡(luò)模型更加合理。本研究使用基于LASSO的Cox比例風(fēng)險回歸模型，對通路富集后的基因進行得分計算，并通過得分矩陣對通路進行篩選[6]，最終給出與卵巢癌復(fù)發(fā)有相關(guān)的12條通路。其中，F(xiàn)oxO信號通路經(jīng)動物實驗證明有抑癌作用[7]，同時，研究表明FoxO蛋白表達量與卵巢癌生存期相關(guān)[8]?；谝陨辖Y(jié)果，本研究選取FoxO信號通路中全部基因，在復(fù)發(fā)組和非復(fù)發(fā)組中進行差異網(wǎng)絡(luò)模型的構(gòu)建與分析。

將映射在FoxO信號通路中的119個基因的表達數(shù)據(jù)整理出來，標(biāo)準(zhǔn)化處理后結(jié)合復(fù)發(fā)狀態(tài)數(shù)據(jù)進行差異網(wǎng)絡(luò)模型的構(gòu)建。應(yīng)用基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法構(gòu)建差異網(wǎng)絡(luò)模型，結(jié)果如圖4所示，其中節(jié)點代表FoxO信號通路中的基因，黑色虛線代表復(fù)發(fā)與非復(fù)發(fā)狀態(tài)下的共同邊，13條黑色虛線代表只存在于復(fù)發(fā)患者中的差異邊，9條黑色實線則代表只存在于非復(fù)發(fā)患者中的差異邊(表2)。存在差異調(diào)控關(guān)系的基因如表3所示，其對應(yīng)的差異調(diào)控關(guān)系頻數(shù)越大，越可能成為差異網(wǎng)絡(luò)中的樞紐基因(Hub Gene)。

圖4 卵巢癌患者復(fù)發(fā)與非復(fù)發(fā)狀態(tài)下的差異基因調(diào)控網(wǎng)絡(luò)

通過查詢GeneMANIA及KEGG基因數(shù)據(jù)庫，發(fā)現(xiàn)有8條差異邊出現(xiàn)在數(shù)據(jù)庫中，例如NLK與PRKAA1、FOXO4和SGK2三個基因間的調(diào)控關(guān)系。表3中NLK、STK4、HOMER1、PRKAA1、EGFR五個基因與多個基因間存在差異調(diào)控關(guān)系，可以將他們視為差異網(wǎng)絡(luò)中的樞紐基因。有文獻報道NLK與卵巢癌分期、分級、化療和預(yù)后有關(guān)，NLK基因的高表達加速了順鉑治療中的細胞凋亡過程，從而延長病人生存期[9]。另一樞紐基因EGFR與腫瘤細胞的增殖、血管生成、腫瘤侵襲、轉(zhuǎn)移及細胞凋亡的抑制有關(guān)，EGFR的表達水平可作為上皮性卵巢癌患者生存期的預(yù)后因子，其高表達與生存期的降低相關(guān)[10]。

表2 卵巢癌患者復(fù)發(fā)狀態(tài)差異網(wǎng)絡(luò)中的差異調(diào)控關(guān)系

表3 卵巢癌患者復(fù)發(fā)狀態(tài)差異網(wǎng)絡(luò)中的差異基因及其對應(yīng)的差異調(diào)控關(guān)系頻數(shù)

討 論

傳統(tǒng)變量篩選方法主要針對不同條件下基因表達量的差異，無法準(zhǔn)確識別表達量未發(fā)生明顯變化但調(diào)控關(guān)系改變的基因。基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法是在高斯圖模型基礎(chǔ)上，給出的一種能夠識別不同條件下調(diào)控關(guān)系改變的變量和調(diào)控邊的一種方法。該方法利用貝葉斯原理，在假定的先驗分布上對模型的參數(shù)進行估計，獲得偏相關(guān)系數(shù)，在此基礎(chǔ)上從整體網(wǎng)絡(luò)的角度進行差異推斷。

模擬實驗結(jié)果表明，差異貝葉斯方法具有較好的識別差異邊的能力，本研究得到以下幾點結(jié)論：①基于貝葉斯準(zhǔn)則的差異網(wǎng)絡(luò)分析方法結(jié)果可靠，即在保證高準(zhǔn)確度的同時，還具有較低假發(fā)現(xiàn)率的特點；②該方法受樣本量的影響較小，在樣本量小于變量數(shù)的情況下也能識別差異網(wǎng)絡(luò)中的差異調(diào)控關(guān)系；③傳統(tǒng)的高斯圖模型推斷差異網(wǎng)絡(luò)時，需要在不同情況下單獨建網(wǎng)后比較差異邊，相比之下，加入貝葉斯因子進行模型選擇，提高了高斯圖模型識別差異邊的能力。

通過實例分析，對上皮性卵巢癌復(fù)發(fā)和非復(fù)發(fā)的調(diào)控網(wǎng)絡(luò)進行了對比，得出兩組在FOXO信號通路中差異調(diào)控關(guān)系及差異網(wǎng)絡(luò)的樞紐基因，為卵巢癌復(fù)發(fā)的基因調(diào)整網(wǎng)絡(luò)機制的研究提供了線索。

應(yīng)用該方法理論上對變量的數(shù)目沒有限制，因此適用任何差異基因調(diào)控網(wǎng)絡(luò)的比較，并能準(zhǔn)確地進行差異調(diào)控邊的定位。這種方法可以通過貝葉斯因子得出差異模型證據(jù)，進而從數(shù)值上檢驗差異存在的可能性大小。然而，生物網(wǎng)絡(luò)通常是非常復(fù)雜的，要確證不同條件下的真實差異調(diào)控關(guān)系，還必須通過生物學(xué)實驗進行驗證。