基于基因擾動及變分逼近技術的基因調控網絡推斷

董自健 宋鐵成 袁 創

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

基因調控網絡(GRN)結構的推導對理解基因功能和復雜疾病的致病原因至關重要．基因調控網絡的推斷需要識別出基因之間的相互作用,而僅使用實驗方法確定基因之間的相互作用需要進行大量的實驗．由于計算技術和計算能力的快速提高,現在通常不是單純使用實驗手段,而是根據有限的實驗數據,使用統計方法推斷基因調控網絡．

由于引入了擾動信息,基因擾動技術可有效提高GRN推斷的精度．基于基因表達數據,可以把基因調控網絡建模為圖結構[1], 其中基因用節點表示,節點之間的邊表示基因之間的調控關系．另外高斯圖模型[2]、線性回歸模型等[3]都可以用于基因調控網絡的推斷．但這些方法的準確度都不高．結構方程(SEM)模型也可用于推斷基因調控網絡[4-5],如Logsdon等[5]使用SEM對基因網絡進行建模,并使用Adaptive Lasso (AL-based)算法推導網絡參數．通過仿真,作者認為AL-based算法優于已有的其他各種算法．

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

1 基因網絡結構的CFA模型

假設有Ng個基因,Nq個標記,N個測量樣本．基因網絡結構服從如下結構方程模型:

Y=AX+BY+ε

(1)

式中,Y為Ng×N矩陣,其中ykj為基因k在第j次測量中所獲得的表達數據;B為Ng×Ng矩陣,包含所有的網絡參數;X為Nq×NeQTLs數據矩陣,Nq≥Ng；A為Ng×Nq矩陣,表示每個eQTLs對每個基因的影響;ε為Ng×N矩陣,其中εkj表示第j次測量中基因k的表達誤差．

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

式(1)中所有未知參數都在B和A中,因此模型可變為

Y=ΛΩ+ε

(2)

式中,Λ=[BA],Ω=[YX]T．現在模型已變為一類特殊形式的SEM模型,即CFA模型．將式(2)作如下分解:

(3)

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

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

(4)

2 CFA模型的變量逼近推斷

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

將式(4)轉化為

(5)

式中,Ωj為Ω的第j行參數,D=2Ng．現在問題轉化為一個標準的線性回歸問題．如前所述,在貝葉斯推斷中,如果采用MCMC方法,在變量較多的情況下,計算量將非常大．因此本文使用變分逼近替代MCMC方法用于網絡調控參數的推斷．變分逼近技術已得到很多應用,如Logsdon等[6]和Carbonetto等[7]把該方法用于多位點的基因組關聯分析．

2.1 調控參數后驗概率計算

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

(6)

其中,每個分項可表示為如下形式[9]:

(7)

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

(8)

(9)

(10)

2.2 基于重要性抽樣技術的包含概率計算

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

(11)

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

(12)

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

(13)

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

(14)

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

(15)

算法具體過程如下:

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

1 根據式(2)計算Ω;

3 fork=1∶Ng(對調控參數矩陣逐行求解)

form=1∶M

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

while flag

forj=1∶D

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

endj;

end while

endm

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

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

根據式(15)計算調控網絡參數Λk;

endk

3 仿真結果

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

本文主要與AL-based算法在檢測率和誤檢率上進行性能比較,圖1給出了仿真結果．可以看出,在檢測率和誤檢率2個方面,VAM算法均優于AL-based算法．

圖1 VAM算法在30個基因調控網絡中的性能

4 酵母基因調控網絡的推斷

本節使用VAM算法對一酵母基因網絡進行推斷．數據來自對酵母菌株BY4716和RM11-1A進行交叉實驗所獲得的112組基因表達數據,以及基因標記數據[11]．數據中共包含5727個基因．由于測量數據樣本有限,且存在一定的誤差,不可能用這些數據較為準確地推斷出全部5727個基因之間的調控關系,因此需要使用預處理技術,選擇具有較強表達數據差異的基因進行分析[5]．本文采用Logsdon等[5]提供的經過預處理的基因表達數據以及相應的eQTLs數據進行網絡參數推導．這些數據共包含35個基因的相關測量數據,具體酵母基因名稱見圖2．

圖2 部分酵母基因調控網絡結構圖

所推導的調控網絡結構如圖2所示．其中8個基因與其他基因之間沒有相互關系,其他27個基因之間有43條有向邊．實線表示基因之間是正調控,共18條,其中ORC5與NUP60之間有相互作用的正調控;虛線表示負調控,共25條邊．

為了驗證算法推導的結果,本文使用YEASTRACT網站的Generate Regulation Matrix工具[12],得到了包含上述35個基因的調控網絡圖,其中包含3個基因調控對,分別是HAL9調控HIM1[13], PHD1調控YEL073C[13],PHD1調控FYV6[14],這些基因調控關系都經過已有實驗驗證．VAM算法推導出的43條邊中,包含了HAL9調控HIM1和PHD1調控YEL073C這2組調控關系對．

5 結語

本文把基因調控網絡建模為沒有潛變量的驗證性因子分析模型,并把模型分解為線性回歸模型,逐行求解．在處理線性回歸模型時,首先采用變分逼近技術推導參數的后驗概率分布,在此基礎上,使用重要性抽樣技術計算參數包含在線性模型中的后驗概率,并得出網絡參數的加權平均．將VAM算法與AL-based算法進行比較,驗證了VAM算法的有效性和可靠性．最后使用實驗數據,針對一種酵母基因調控網絡進行了推導,得出了一個具有43條邊的參考基因調控網絡結構．

)

[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.