基于基因擾動(dòng)及變分逼近技術(shù)的基因調(diào)控網(wǎng)絡(luò)推斷

董自健 宋鐵成 袁 創(chuàng)

(1東南大學(xué)信息科學(xué)與工程學(xué)院, 南京 210096)(2淮海工學(xué)院電子工程學(xué)院, 連云港 222005)(3香港理工大學(xué)醫(yī)療科技及信息學(xué)系, 香港)

基因調(diào)控網(wǎng)絡(luò)(GRN)結(jié)構(gòu)的推導(dǎo)對(duì)理解基因功能和復(fù)雜疾病的致病原因至關(guān)重要．基因調(diào)控網(wǎng)絡(luò)的推斷需要識(shí)別出基因之間的相互作用,而僅使用實(shí)驗(yàn)方法確定基因之間的相互作用需要進(jìn)行大量的實(shí)驗(yàn)．由于計(jì)算技術(shù)和計(jì)算能力的快速提高,現(xiàn)在通常不是單純使用實(shí)驗(yàn)手段,而是根據(jù)有限的實(shí)驗(yàn)數(shù)據(jù),使用統(tǒng)計(jì)方法推斷基因調(diào)控網(wǎng)絡(luò)．

由于引入了擾動(dòng)信息,基因擾動(dòng)技術(shù)可有效提高GRN推斷的精度．基于基因表達(dá)數(shù)據(jù),可以把基因調(diào)控網(wǎng)絡(luò)建模為圖結(jié)構(gòu)[1], 其中基因用節(jié)點(diǎn)表示,節(jié)點(diǎn)之間的邊表示基因之間的調(diào)控關(guān)系．另外高斯圖模型[2]、線性回歸模型等[3]都可以用于基因調(diào)控網(wǎng)絡(luò)的推斷．但這些方法的準(zhǔn)確度都不高．結(jié)構(gòu)方程(SEM)模型也可用于推斷基因調(diào)控網(wǎng)絡(luò)[4-5],如Logsdon等[5]使用SEM對(duì)基因網(wǎng)絡(luò)進(jìn)行建模,并使用Adaptive Lasso (AL-based)算法推導(dǎo)網(wǎng)絡(luò)參數(shù)．通過仿真,作者認(rèn)為AL-based算法優(yōu)于已有的其他各種算法．

基于基因表達(dá)和eQTLs數(shù)據(jù),本文把基因調(diào)控網(wǎng)絡(luò)建模為SEM結(jié)構(gòu)．由于這2類數(shù)據(jù)均已知,故不存在潛變量,因此本文把SEM模型轉(zhuǎn)化為驗(yàn)證性因子分析(CFA)結(jié)構(gòu),并進(jìn)而簡(jiǎn)化為線性回歸模型．為降低算法對(duì)數(shù)據(jù)的敏感性,本文使用貝葉斯推斷方法代替常用的最大似然以及各種Lasso方法求解線性回歸模型;為了降低計(jì)算量,不使用通常的馬爾科夫-蒙特卡洛(MCMC)方法,而是采用變分逼近方法(VAM)進(jìn)行貝葉斯推斷．仿真結(jié)果表明CFA模型和VAM算法都是有效和可靠的．VAM算法的性能較好,在誤檢率和檢測(cè)率上都遠(yuǎn)高于AL-based算法;與最大似然等方法相比,變分逼近方法無需準(zhǔn)確估計(jì)初始參數(shù)，對(duì)數(shù)據(jù)適應(yīng)性較好．

1 基因網(wǎng)絡(luò)結(jié)構(gòu)的CFA模型

假設(shè)有Ng個(gè)基因,Nq個(gè)標(biāo)記,N個(gè)測(cè)量樣本．基因網(wǎng)絡(luò)結(jié)構(gòu)服從如下結(jié)構(gòu)方程模型:

Y=AX+BY+ε

(1)

式中,Y為Ng×N矩陣,其中ykj為基因k在第j次測(cè)量中所獲得的表達(dá)數(shù)據(jù);B為Ng×Ng矩陣,包含所有的網(wǎng)絡(luò)參數(shù);X為Nq×NeQTLs數(shù)據(jù)矩陣,Nq≥Ng；A為Ng×Nq矩陣,表示每個(gè)eQTLs對(duì)每個(gè)基因的影響;ε為Ng×N矩陣,其中εkj表示第j次測(cè)量中基因k的表達(dá)誤差．

假設(shè)基因自身不存在相互作用,因此矩陣B的對(duì)角元素值為0．另外,假設(shè)所有的eQTLs位點(diǎn)已根據(jù)已有的方法確定出來,但是每個(gè)eQTLs的影響未知,因此A中有Nq個(gè)位置已知的未知參數(shù),剩下的NgNq-Nq個(gè)元素都是0．不失一般性,假設(shè)Ng=Nq．

式(1)中所有未知參數(shù)都在B和A中,因此模型可變?yōu)?/p>

Y=ΛΩ+ε

(2)

式中,Λ=[BA],Ω=[YX]T．現(xiàn)在模型已變?yōu)橐活愄厥庑问降腟EM模型,即CFA模型．將式(2)作如下分解:

(3)

由于已知數(shù)據(jù)均在Ω和Y中,所以可逐行求解Λ中參數(shù),因此問題分解為

Yk=ΛkΩ+Υkk=1,2,…,Ng

(4)

2 CFA模型的變量逼近推斷

假設(shè)Λ矩陣稀疏,即Λ中大部分?jǐn)?shù)據(jù)為0,因此可假設(shè)Λ中所有參數(shù)服從均值為0的高斯分布,且Υk服從N(0,σ2I)．

將式(4)轉(zhuǎn)化為

(5)

式中,Ωj為Ω的第j行參數(shù),D=2Ng．現(xiàn)在問題轉(zhuǎn)化為一個(gè)標(biāo)準(zhǔn)的線性回歸問題．如前所述,在貝葉斯推斷中,如果采用MCMC方法,在變量較多的情況下,計(jì)算量將非常大．因此本文使用變分逼近替代MCMC方法用于網(wǎng)絡(luò)調(diào)控參數(shù)的推斷．變分逼近技術(shù)已得到很多應(yīng)用,如Logsdon等[6]和Carbonetto等[7]把該方法用于多位點(diǎn)的基因組關(guān)聯(lián)分析．

2.1 調(diào)控參數(shù)后驗(yàn)概率計(jì)算

根據(jù)變分理論[8],限制p(θ|Ω,Yk,υ(m))≈q(Λk,Zk),且q(Λk,Zk)可分解為

(6)

其中,每個(gè)分項(xiàng)可表示為如下形式[9]:

(7)

根據(jù)變分逼近理論可得出迭代參數(shù)為[7,10]

(8)

(9)

(10)

2.2 基于重要性抽樣技術(shù)的包含概率計(jì)算

為考察基因k和基因j之間的關(guān)系,還需判斷每個(gè)參數(shù)(Λkj)是否包含在模型中,即計(jì)算

(11)

該積分計(jì)算非常復(fù)雜,需要近似計(jì)算．已知p(zkj=1|Ω,Yk,θ)=αkj,為了計(jì)算Λkj是否包含在模型中,需要根據(jù)式(11)計(jì)算其后驗(yàn)概率．這里使用重要性抽樣,即

(12)

式中,w(υ(m))為重要性抽樣系數(shù),

(13)

式(13)中的p(Yk|Ω,υ)很難計(jì)算,Carbonetto等[7]使用Jessen不等式獲得了其下界:

(14)

式中,αk.×μk表示αk和μk兩個(gè)行向量的點(diǎn)乘．在計(jì)算時(shí)可使用這個(gè)下界代替p(Yk|Ω,υ)．

(15)

算法具體過程如下:

輸入:X,Y,抽樣樣本序列{sπ,sσ,sΛ},以及α(0),μ(0),ξ．

1 根據(jù)式(2)計(jì)算Ω;

3 fork=1∶Ng(對(duì)調(diào)控參數(shù)矩陣逐行求解)

form=1∶M

置υ=υ(m);令flag=1;根據(jù)式(8)計(jì)算s2;

while flag

forj=1∶D

根據(jù)式(9)、(10)分別更新μkj和αkj

endj;

end while

endm

根據(jù)式(13)、(14)獲得重要性系數(shù)w(υ(m));

根據(jù)式(12)獲得后驗(yàn)概率p(zkj=1|Ω,Yk);

根據(jù)式(15)計(jì)算調(diào)控網(wǎng)絡(luò)參數(shù)Λk;

endk

3 仿真結(jié)果

本文分別仿真了含有30個(gè)基因的有向無環(huán)和有向有環(huán)網(wǎng)絡(luò)．每個(gè)基因節(jié)點(diǎn)平均有3條邊,即NE=3．如果從節(jié)點(diǎn)j到節(jié)點(diǎn)k有邊,則從(-1,-0.5)∪(0.5,1)中均勻抽樣;否則Bkj置為0．eQTL中每個(gè)xij以概率(0.25,0.5,0.25)從集合 {1,2,3}中抽取．測(cè)量噪聲服從N(0,0.1)分布．不失一般性,設(shè)A為單位陣．Y通過式 (1)計(jì)算獲得．在仿真中,如果算法返回的p(zkj=1|Ω,Yk)>0.95,且Λkj>0.10,則認(rèn)為從節(jié)點(diǎn)j到節(jié)點(diǎn)k之間有一條邊,其他情況均判為2個(gè)節(jié)點(diǎn)之間無調(diào)控關(guān)系．

本文主要與AL-based算法在檢測(cè)率和誤檢率上進(jìn)行性能比較,圖1給出了仿真結(jié)果．可以看出,在檢測(cè)率和誤檢率2個(gè)方面,VAM算法均優(yōu)于AL-based算法．

圖1 VAM算法在30個(gè)基因調(diào)控網(wǎng)絡(luò)中的性能

4 酵母基因調(diào)控網(wǎng)絡(luò)的推斷

本節(jié)使用VAM算法對(duì)一酵母基因網(wǎng)絡(luò)進(jìn)行推斷．?dāng)?shù)據(jù)來自對(duì)酵母菌株BY4716和RM11-1A進(jìn)行交叉實(shí)驗(yàn)所獲得的112組基因表達(dá)數(shù)據(jù),以及基因標(biāo)記數(shù)據(jù)[11]．?dāng)?shù)據(jù)中共包含5727個(gè)基因．由于測(cè)量數(shù)據(jù)樣本有限,且存在一定的誤差,不可能用這些數(shù)據(jù)較為準(zhǔn)確地推斷出全部5727個(gè)基因之間的調(diào)控關(guān)系,因此需要使用預(yù)處理技術(shù),選擇具有較強(qiáng)表達(dá)數(shù)據(jù)差異的基因進(jìn)行分析[5]．本文采用Logsdon等[5]提供的經(jīng)過預(yù)處理的基因表達(dá)數(shù)據(jù)以及相應(yīng)的eQTLs數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)參數(shù)推導(dǎo)．這些數(shù)據(jù)共包含35個(gè)基因的相關(guān)測(cè)量數(shù)據(jù),具體酵母基因名稱見圖2．

圖2 部分酵母基因調(diào)控網(wǎng)絡(luò)結(jié)構(gòu)圖

所推導(dǎo)的調(diào)控網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示．其中8個(gè)基因與其他基因之間沒有相互關(guān)系,其他27個(gè)基因之間有43條有向邊．實(shí)線表示基因之間是正調(diào)控,共18條,其中ORC5與NUP60之間有相互作用的正調(diào)控;虛線表示負(fù)調(diào)控,共25條邊．

為了驗(yàn)證算法推導(dǎo)的結(jié)果,本文使用YEASTRACT網(wǎng)站的Generate Regulation Matrix工具[12],得到了包含上述35個(gè)基因的調(diào)控網(wǎng)絡(luò)圖,其中包含3個(gè)基因調(diào)控對(duì),分別是HAL9調(diào)控HIM1[13], PHD1調(diào)控YEL073C[13],PHD1調(diào)控FYV6[14],這些基因調(diào)控關(guān)系都經(jīng)過已有實(shí)驗(yàn)驗(yàn)證．VAM算法推導(dǎo)出的43條邊中,包含了HAL9調(diào)控HIM1和PHD1調(diào)控YEL073C這2組調(diào)控關(guān)系對(duì)．

5 結(jié)語(yǔ)

本文把基因調(diào)控網(wǎng)絡(luò)建模為沒有潛變量的驗(yàn)證性因子分析模型,并把模型分解為線性回歸模型,逐行求解．在處理線性回歸模型時(shí),首先采用變分逼近技術(shù)推導(dǎo)參數(shù)的后驗(yàn)概率分布,在此基礎(chǔ)上,使用重要性抽樣技術(shù)計(jì)算參數(shù)包含在線性模型中的后驗(yàn)概率,并得出網(wǎng)絡(luò)參數(shù)的加權(quán)平均．將VAM算法與AL-based算法進(jìn)行比較,驗(yàn)證了VAM算法的有效性和可靠性．最后使用實(shí)驗(yàn)數(shù)據(jù),針對(duì)一種酵母基因調(diào)控網(wǎng)絡(luò)進(jìn)行了推導(dǎo),得出了一個(gè)具有43條邊的參考基因調(diào)控網(wǎng)絡(luò)結(jié)構(gòu)．

)

[1] Friedman N. Inferring cellular networks using probabilistic graphical models [J].Science, 2004,303(5659):799-805.

[2] Schafer J, Strimmer K. An empirical Bayes approach to inferring large-scale gene association networks [J].Bioinformatics, 2005,21(6):754-764.

[3] Schafer J, Strimmer K. A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics [J].StatisticalApplicationsinGeneticsandMolecularBiology, 2005,4(1):1175-1189.

[4] Xiong H, Choe Y. Structural systems identification of genetic regulatory networks [J].Bioinformatics, 2008,24(4):553-560.

[5] Logsdon B A, Mezey J. Gene expression network reconstruction by convex feature selection when incorporating genetic perturbations [J/OL].PLoSComputationalBiology,2010,6(12). http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1001014.

[6] Logsdon B A, Hoffman G E, Mezey J G. A variational Bayes algorithm for fast and accurate multiple locus genome-wide association analysis [J/OL].BMCBioinformatics, 2010,11. http://www.biomedcentral.com/1471-2105/11/58/.

[7] Carbonetto P, Stephens M. Scalable variational inference for Bayesian variable selection in regression, and its accuracy in genetic association studies [J].BayesianAnalysis, 2012,7(1):73-107.

[8] Jordan M I, Ghahramani Z, Jaakkola T S, et al. An introduction to variational methods for graphical models [J].MachineLearning, 1999,37(2):183-233.

[9] Attias H. Independent factor analysis [J].NeuralComputation, 1999,11(4):803-851.

[10] Beal M J. Variational algorithms for approximate Bayesian inference [D]. London: University of London, 2003.

[11] Stranger B E, Forrest M S, Dunning M, et al. Relative impact of nucleotide and copy number variation on gene expression phenotypes [J].Science, 2007,315(5813):848-853.

[12] Knowledge Discovery & Bioinformatics Research Group.YEASTRACT [EB/OL]. [2013-03-10]. http://www.yeastract.com/index.php.

[13] Lee T I, Rinaldi N J, Robert F, et al. Transcriptional regulatory networks in Saccharomyces cerevisiae [J].Science, 2002,298(5594):799-804.

[14] Borneman A R, Leigh-Bell J A, Yu H, et al. Target hub proteins serve as master regulators of development in yeast [J].GenesandDevelopment, 2006,20(4):435-448.