Aderhold, A.

Dondelinger, F.

Grzegorczyk, M.

Husmeier, D.

Lebre, S.

English

Enlighten: Publications

A Bayesian regression and multiple changepoint model for systems biology

Enlighten

         Dondelinger, F., Aderhold, A., Lebre, S., Grzegorczyk, M., and Husmeier, D. (2011) A Bayesian regression and multiple changepoint model for systems biology. In: Conesa, D., Forte, A., Lopez-Quilez, A. and Munoz, F. (eds.) International Workshop on Statistical Modelling. Copiformes S.L., Valencia, Spain, pp. 189-194. ISBN 9788469451298   Copyright © 2011 The Authors   http://eprints.gla.ac.uk/69384        Deposited on: 15 February 2013                    Enlighten – Research publications by members of the University of Glasgow http://eprints.gla.ac.uk A Bayesian regression and multiplechangepoint model for systems biologyFrank Dondelinger1,2,6, Andrej Aderhold1, Sophie Le`bre3,Marco Grzegorczyk4,5, and Dirk Husmeier11 Biomathematics and Statistics Scotland, JCMB, Edinburgh, EH9 3JZ, UK.2 Institute for Adaptive and Neural Computation, School of Informatics, Univer-sity of Edinburgh, Edinburgh, EH8 9AB, UK3 LSIIT - UMR 7005, Universite´ de Strasbourg, 67412 Illkirch, France.4 Department of Statistics, TU Dortmund, 44221 Dortmund, Germany5 Department of Mathematics, Carl von Ossietzky University Oldenburg, Ger-many6 Communicating Author. Email: frankd@bioss.ac.ukAbstract: We propose a Bayesian regression and multiple changepoint modelfor reverse engineering gene regulatory networks from high-throughput gene ex-pression proﬁles. We report results from a recently held international gene net-work reconstruction competition, in which our method was objectively assessedin a blind study. While we did not win the competition, the scores indicate thatthe proposed method favourably compares with the majority of competing ap-proaches and clearly belongs to the group of highest-ranked performers.Keywords: Systems biology; gene regulatory network inference; Bayesian mul-tiple changepoint model; RJMCMC; DREAM1 IntroductionThe objective of the highly topical ﬁeld of systems biology is the reverseengineering of molecular regulatory networks and signalling pathways fromhigh-throughput post-genomic data, and a ﬂurry of activities in the statis-tics and machine learning communities are currently aimed at solving thisproblem. A variety of methods from statistics and machine learning havebeen applied to this end. See e.g. Grzegorczyk et al. (2008) and Cantoneet al. (2009) for brief reviews. In the present paper, we propose a Bayesianregression and multiple changepoint model, with Bayesian inference basedon reversible jump Markov chain Monte Carlo (RJMCMC) (Green, 1995).We participated in a recently held gene regulatory network prediction com-petition (DREAM 5), which assures that the comparative evaluation withother methods was done objectively.2 Bayesian model for gene regulatory networks2 ModelMultiple changepoints: Let 푝 be the number of target genes, whose ex-pression values 푦 = {푦푖(푡)}1≤푖≤푝,1≤푡≤푁 are measured on 푁 separate chips.ℳ푖 is the set of parents (regulators) associated with target gene 푖 in thegene regulatory network. We model the diﬀerences in the regulatory re-lationships measured by diﬀerent chips (assumed to be in some naturalorder, e.g. a time series) with a multiple changepoint process. For eachtarget gene 푖, an unknown number 푘푖 of changepoints deﬁne 푘푖 + 1 non-overlapping segments. Segment ℎ ∈ {1, .., 푘푖+1} starts at changepoint 휉ℎ−1푖and stops before 휉ℎ푖 , so that 휉푖 = (휉0푖 , ..., 휉ℎ−1푖 , 휉ℎ푖 , ..., 휉푘푖+1푖 ) with 휉ℎ−1푖 < 휉ℎ푖 .This changepoint process induces a partition of the chip ordering, 푦ℎ푖 =(푦푖(푡))휉ℎ−1푖 ≤푡<휉ℎ푖. The network structure ℳ푖 remains the same for each seg-ment ℎ, but the other parameters of the model can vary.Regression model: For all genes 푖, the random variable 푌푖(푡) refers to theexpression of gene 푖 on chip 푡. Within any segment ℎ, the expression of gene푖 at chip 푡 depends on the gene expression values on chip 푡 of a set 푅푖 of푚 potential regulator genes (parents), with 푖 /∈ 푅푖. We deﬁne a regressionmodel by (a) the set of 푠푖 parents denoted by ℳ푖 = {푗1, ..., 푗푠푖} ⊆ 푅푖, and(b) a set of parameters ((푎ℎ푖푗)푗∈푅푖 , 휎ℎ푖 ); 푎ℎ푖푗 ∈ ℝ, 휎ℎ푖 > 0. For all 푗 ∕= 0, 푎ℎ푖푗 = 0if 푗 /∈ ℳ푖. For all genes 푖, for all chips 푡 in segment ℎ (휉ℎ−1푖 ≤ 푡 < 휉ℎ푖 ), therandom variable 푌푖(푡) depends on the 푚 variables {푌푗(푡)}푗∈푅푖 according to푌푖(푡) = 푎ℎ푖0 +∑푗∈ℳ푖푎ℎ푖푗 푌푗(푡) + 휀푖(푡) (1)where the noise 휀푖(푡) is assumed to be Gaussian with mean 0 and variance(휎ℎ푖 )2, 휀푖(푡) ∼ 푁(0, (휎ℎ푖 )2). We deﬁne 푎ℎ푖 = (푎ℎ푖푗)푗∈푅푖 .Prior: The 푘푖+1 segments are delimited by 푘푖 changepoints, where 푘푖 is dis-tributed a priori as a truncated Poisson random variable with mean 휆 andmaximum 푘 = 푁−2: 푃 (푘푖∣휆) ∝ 휆푘푖푘푖! 1l{푘푖≤푘} . Conditional on 푘푖 changepoints,the changepoint positions vector 휉푖 = (휉0푖 , 휉1푖 , ..., 휉푘푖+1푖 ) takes non-overlappinginteger values, which we take to be uniformly distributed a priori. For allgenes 푖, the number 푠푖 of parents for node 푖 follows a truncated Poissondistribution with mean Λ and maximum 푠 = 5: 푃 (푠푖∣Λ) ∝ Λ푠푖푠푖! 1l{푠푖≤푠}. Con-ditional on 푠푖, the prior for the parent setℳ푖 is a uniform distribution overall parent sets with cardinality 푠푖: 푃 (ℳ푖 ∣∣ℳ푖∣ = 푠푖) = 1/( 푝푠푖). The overallprior on the network structures is given by marginalization:푃 (ℳ푖∣Λ) =∑푠푠푖=1푃 (ℳ푖∣푠푖)푃 (푠푖∣Λ) (2)Conditional on the parent setℳ푖 of size 푠푖, we assume for the prior distri-bution 푃 (푎ℎ푖 ∣ℳ푖, 휎ℎ푖 ) of the 푠푖+1 regression coeﬃcients for each segment ℎa zero-mean multivariate Gaussian with covariance matrix (휎ℎ푖 )2Σ푎ℎ푖, wherefollowing Andrieu and Doucet (1999) we set Σ푎ℎ푖 = 훿−2퐷†푎ℎ푖(푦)퐷푎ℎ푖(푦), and퐷푎ℎ푖(푦) is the (휉ℎ푖 − 휉ℎ−1푖 ) × (푠푖 + 1) matrix whose ﬁrst column is a vectorof 1 (for the constant in model (1)) and each (푗 + 1)푡ℎ column containsthe observed values (푦푗(푡))휉ℎ−1푖 −1≤푡<휉ℎ푖 −1for all regulatory genes 푗 in ℳ푖.Dondelinger et al. 3Finally, the conjugate prior for the variance (휎ℎ푖 )2 is the inverse gamma dis-tribution, 푃 ((휎ℎ푖 )2) = ℐ풢(휐0, 훾0). Following Le`bre et al. (2010), we set thehyperparameters for shape, 휐0 = 0.5, and scale, 훾0 = 0.05, to ﬁxed valuesthat give a vague distribution. The terms 휆 and Λ can be interpreted asthe expected number of changepoints and parents, respectively, and 훿2 isthe expected signal-to-noise ratio. These hyperparameters are drawn fromvague conjugate hyperpriors, which are in the (inverse) gamma distributionfamily: 푃 (Λ) = 푃 (휆) = 풢푎(0.5, 1) and 푃 (훿2) = ℐ풢(2, 0.2).Posterior: Equation (1) implies that푃 (푦ℎ푖 ∣휉ℎ−1푖 , 휉ℎ푖 ,ℳ푖, 푎ℎ푖 , 휎ℎ푖 ) ∝ exp⎛⎝− (푦ℎ푖 −퐷푎ℎ푖 (푦)푎ℎ푖 )† (푦ℎ푖 −퐷푎ℎ푖 (푦)푎ℎ푖 )2(휎ℎ푖 )2⎞⎠ (3)From Bayes theorem, the posterior is given by the following equation:푃 (푘, 휉,ℳ, 푎, 휎, 휆,Λ, 훿2∣푦) ∝ 푃 (훿2)푃 (휆)푃 (Λ)푝∏푖=1푃 (푘푖∣휆)푃 (휉푖∣푘푖)푃 (ℳ푖∣Λ) (4)푘푖∏ℎ=1푃 ([휎ℎ푖 ]2)푃 (푎ℎ푖 ∣ℳ푖, [휎ℎ푖 ]2, 훿2)푃 (푦ℎ푖 ∣휉ℎ−1푖 , 휉ℎ푖 ,ℳ푖, 푎ℎ푖 , [휎ℎ푖 ]2)Inference: An attractive feature of the chosen model is that the marginal-ization over the parameters 푎 and 휎 in the posterior distribution of (4)is analytically tractable: 푃 (푘,휉,ℳ,휆,Λ,훿2∣푦) = ∫ 푃 (푘,휉,ℳ,푎,휎,휆,Λ,훿2∣푦)푑푎푑휎See Andrieu and Doucet (1999), Le`bre et al. (2010) for details and an ex-plicit expression. The number of changepoints and their location, 푘, 휉, thenetwork structureℳ and the hyperparameters 휆, Λ and 훿2 can be sampledfrom the posterior 푃 (푘, 휉,ℳ, 휆,Λ, 훿2∣푦) with RJMCMC. A detailed descrip-tion can be found in Le`bre et al. (2010). The posterior probabilities of thegene interactions submitted to DREAM are obtained from the posteriorsample of network structures ℳ by marginalization.0 5 10 15 20 25 3000.050.10.150.20.250.30.35Team Prediction Performance (AUPRC)Team NumberAUPRC Score0 5 10 15 20 25 300.50.550.60.650.70.750.80.850.9Team Prediction Performance (AUROC)Team NumberAUROC ScoreFIGURE 1. Areas under the precision recall (left) and ROC (right) curves ob-tained on an in silico data set by all teams participating in the DREAM 5 com-petition. The circles indicate the performance of our proposed method.3 Simulations and ResultsTo assess the performance of the proposed method we participated in acompetition organised by the DREAM (Dialogue for Reverse Engineering4 Bayesian model for gene regulatory networksTABLE 1. This table summarises the information about the DREAM 5 NetworkInference Challenge data sets. For each data set, we show which organism it camefrom, how many genes were measured, how many of those genes were identiﬁed astranscription factors (possibly regulatory genes) and how many chips (datapoints)were included.Data Set Organism Genes Transcription Factors Chips1 Synthetic 1643 195 8062 S. Aureus 2810 99 1603 E. Coli 4511 334 8054 S. Cerevisiae 5950 333 536Assessments and Methods) consortium in autumn of 2010. The goal was toreverse engineer gene regulatory networks from gene expression data sets.Participants were given four microarray compendia and were challenged toinfer the structure of the underlying transcriptional regulatory networks.The ﬁrst compendium was based on an in-silico (i.e. simulated) network,the other three compendia were obtained from microorganisms. Each com-pendium consisted of hundreds of microarray experiments, which includeda wide range of genetic, drug, and environmental perturbations. More infor-mation is available in Table 1 and at http://wiki.c2b2.columbia.edu/dream/index.php/The_DREAM_Project. Network predictions were evalu-ated by the organisers on a subset of known interactions for each organism,or on the known network for the in-silico case (which is more objective).Our method assumes an ordering of the microarray chips. While this con-dition is naturally met for time course experiments, it does not hold forthe varying experimental conditions of the DREAM data. We therefore re-sorted to the heuristic pre-processing step of mapping the high-dimensionalgene expression proﬁles onto a one-dimensional self-organising map (SOM)initialized by the ﬁrst principal component. We applied the software pack-age som in R with default parameter settings. To reduce the computa-tional complexity of the RJMCMC simulations we applied a pre-ﬁlteringstep based on TESLA (Ahmed and Xing, 2009), a time-varying networkinference method based on L1-regularised linear regression. For each genewe identiﬁed a set of 20 potential candidate regulators, based on the 20regression coeﬃcients with the largest modulus.We assessed the convergence of our simulations with standard diagnosticsbased on Gelman-Rubin potential scale reduction factors (PSRF). Owingto unexpected downtime of the computer cluster we were using, only thesimulations on the ﬁrst two data sets showed a suﬃcient degree of conver-gence (PSRF≤ 1.2); for the latter data sets we submitted the results fromTESLA. The second data set was later removed from the evaluation by theorganisers. Figure 1 shows the results for the in silico data set obtained fromthe rankings of interactions submitted by all participating teams, using twoDondelinger et al. 5criteria: the area under the precision-recall curve (AUPRC), and the areaunder the receiver-operator characteristic (AUROC) curve. As discussedin Davis and Goadrich (2006), AUPRC gives a more faithful indicationof the network reconstruction accuracy than AUROC, and it is thus seenthat our method clearly lies in the group of the 5 top-ranked models. Thissuggests that it compares favourably with the majority of existing schemesand provides a useful tool for contemporary research in systems biology.Acknowledgments: Marco Grzegorczyk is supported by the GraduateSchool Statistische Modellbildung of the Department of Statistics, TUDortmund University. Dirk Husmeier is supported by the Scottish Govern-ment Rural and Environment Research and Analysis Directorate (RERAD)and under the EU FP7 project TiMet’. Andrej Aderhold’s involvement inthis project was funded by RERAD. Frank Dondelinger’s research is fundedby RERAD and the UK Engineering and Physical Sciences Research Coun-cil (EPSRC). We are grateful to Tony Travis for granting us user time onthe Beowulf cluster at the Rowett Institute in Aberdeen and for providingsupport with respect to a parallelisation of the processes.ReferencesAhmed, A., and Xing, E.P. (2009). Recovering time-varying networks ofdependencies in social and biological studies. Proceedings of the Na-tional Academy of Sciences, 106:29, 11878-11883.Andrieu, C., and Doucet, A. (1999). Joint Bayesian model selection andestimation of noisy sinusoids via reversible jumpMCMC. IEEE Trans-actions on Signal Processing, 47, 2667-2676.Cantone, I., Marucci, L., Iorio, F., Ricci, M. A., Belcastro, V., Bansal, M.,Santini, S., di Bernardo, M., di Bernardo, D., and Cosma, M. P.(2009). A yeast synthetic network for in vivo assessment of reverse-engineering and modeling approaches. Cell, 4:130.Davis, J. and Goadrich, M. (2006). The relationship between precision-recalland ROC curves. Proc. of the 23rd Int. Conf. on Machine LearningGreen, P.J. (1995). Reversible jump Markov chain Monte Carlo computa-tion and Bayesian model determination. Biometrika, 82:4, 711-732.Grzegorczyk, M., Husmeier, D. and Werhli, A.V. (2008). Reverse engineer-ing gene regulatory networks with various machine learning methods.In: Analysis of Microarray Data: A Network-Based Approach, WileyOnline Library.Le`bre, S., Becq, J., Devaux, F., Stumpf, M.P.H. and Lelandais, G. (2010).Statistical inference of the time-varying structure of gene regulationnetworks. BMC Systems Biology, 137, 172-181.

A yeast synthetic network for in vivo assessment of reverseengineering and modeling approaches.

Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC.

Recovering time-varying networks of dependencies in social and biological studies.

Reverse engineering gene regulatory networks with various machine learning methods. In: Analysis of Microarray Data: A Network-Based Approach,

Statistical inference of the time-varying structure of gene regulation networks.

The relationship between precision-recall and ROC curves.

http://eprints.gla.ac.uk/69384/1/69384.pdf

A Bayesian regression and multiple changepoint model for systems biology

Abstract

Similar works

Full text

Available Versions

Enlighten: Publications

Enlighten