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

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

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

原理與方法

1.高斯圖模型

(1)

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

在處理高維組學(xué)數(shù)據(jù)時(shí)，往往會(huì)遇到“m?n”的問(wèn)題，即變量數(shù)m遠(yuǎn)大于樣本數(shù)n，此時(shí)協(xié)方差矩陣不是唯一的，因此傳統(tǒng)的統(tǒng)計(jì)模型無(wú)法求出偏相關(guān)系數(shù)。針對(duì)此種情況，有學(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個(gè)變量的數(shù)據(jù)按K個(gè)分類(lèi)標(biāo)簽分成不同數(shù)據(jù)集。對(duì)于公式(2)，記Y=Xi。根據(jù)極大似然和平均場(chǎng)近似估計(jì)理論，可以直接得出K個(gè)獨(dú)立網(wǎng)絡(luò)結(jié)構(gòu)滿(mǎn)足：

(6)

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

(7)

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

(8)

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

模擬實(shí)驗(yàn)

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

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

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

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

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

圖2 模擬實(shí)驗(yàn)各評(píng)價(jià)指標(biāo)隨樣本量不同的變化情況

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

實(shí)例分析

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

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

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

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

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

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

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

討 論

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

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

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

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