Skip to main content

Integrative analysis of time course metabolic data and biomarker discovery



Metabolomics time-course experiments provide the opportunity to understand the changes to an organism by observing the evolution of metabolic profiles in response to internal or external stimuli. Along with other omic longitudinal profiling technologies, these techniques have great potential to uncover complex relations between variations across diverse omic variables and provide unique insights into the underlying biology of the system. However, many statistical methods currently used to analyse short time-series omic data are i) prone to overfitting, ii) do not fully take into account the experimental design or iii) do not make full use of the multivariate information intrinsic to the data or iv) are unable to uncover multiple associations between different omic data. The model we propose is an attempt to i) overcome overfitting by using a weakly informative Bayesian model, ii) capture experimental design conditions through a mixed-effects model, iii) model interdependencies between variables by augmenting the mixed-effects model with a conditional auto-regressive (CAR) component and iv) identify potential associations between heterogeneous omic variables by using a horseshoe prior.


We assess the performance of our model on synthetic and real datasets and show that it can outperform comparable models for metabolomic longitudinal data analysis. In addition, our proposed method provides the analyst with new insights on the data as it is able to identify metabolic biomarkers related to treatment, infer perturbed pathways as a result of treatment and find significant associations with additional omic variables. We also show through simulation that our model is fairly robust against inaccuracies in metabolite assignments. On real data, we demonstrate that the number of profiled metabolites slightly affects the predictive ability of the model.


Our single model approach to longitudinal analysis of metabolomics data provides an approach simultaneously for integrative analysis and biomarker discovery. In addition, it lends better interpretation by allowing analysis at the pathway level. An accompanying R package for the model has been developed using the probabilistic programming language Stan. The package offers user-friendly functions for simulating data, fitting the model, assessing model fit and postprocessing the results. The main aim of the R package is to offer freely accessible resources for integrative longitudinal analysis for metabolomics scientists and various visualization functions easy-to-use for applied researchers to interpret results.


Over the past few years, there has been a significant development in high-throughput omics technologies e.g. metabolomics, transcriptomics, genomics, epigenomics and proteomics along with a growing interest into joint modeling of multi-omic data [1, 2]. In metabolomics, several approaches are used to understand the response of a biological system as a function of an internal or external perturbation by monitoring “the chemical fingerprints that specific cellular processes leave behind" [3]. These chemical fingerprints are most commonly interrogated in terms of metabolite (i.e. low weight molecules) concentration, structure and transformation pathways (i.e set of chemical reactions) in order to identify biomarkers related to the studied process. Biomarker discovery consists of identifying a metabolite that has significant association patterns with a particular phenotype (disease, clinical variables, physical trait, etc) and that can be thus used as an indicator of that specific phenotype. Typical experimental platforms use analytical techniques such as nuclear magnetic resonance spectroscopy (NMR) [4] and mass spectrometry (MS) [5] to generate appropriate spectral metabolomic profiles of the studied biological system.

Metabolomic datasets are characterized by high correlation structures at different levels. For example, chromatographic correlation between two spectral peaks often results from adduct or isotopic effects whereas other correlation structures can result from peaks that represent related molecules operating within networks of chemical reactions in multiple injections. In addition, further correlation structure is present in longitudinal metabolomic studies due to repeated measurements of observations over time. Additional challenges include not only the low number of time points and samples compared to the number of metabolic variables, but also integration of a different omic data with the metabolomic data.

Metabolomic time series are often short due to experimental costs or ethical considerations. Typically, less than 10 time points are available compared to a large number of metabolic variables profiled at each time point e.g. hundreds of metabolic variables for targeted experiments and thousands for untargeted experiments. For this reason, the number of temporal patterns that can arise is limited (due to the limited number of degrees of freedom). Some temporal patterns may be repeated and thus these patterns can be induced by random processes. Models fitted to a small number of data points are prone to overfitting i.e. the model is very sensitive to small fluctuations. This can lead to a poor fitting to new data and thus a high generalisation error. It is also important to consider the number of parameters of the statistical model and make use of the simplest models in order to avoid overparametrisation. Finally, monitoring variables within and between multiple omic types can substantially enhance the understanding of the underlying biological mechanism and provide a systems biology approach as these omic variables represent entities that are often involved in related cellular processes [1, 2]. For all these reasons, metabolomics scientists need robust models which allow cautious interpretation of the data. Models are needed which integrate heterogeneous omic data and take into account both experimental conditions and biological variation.

There is a growing interest in longitudinal experiments for heterogeneous omics data and statistical models to infer biomarkers of a particular treatment or disease over time. Some approaches aim to infer influential or significant metabolites using dynamic metabolomic data under the assumption that metabolites are independent. These models include fitting smooth splines mixed effects models (SME) to time curves [6] and linear mixed effects models augmented with a variable selection [7]. However metabolomic data exhibit rich correlation structure, which is biologically relevant and should be modelled.

Seemingly unrelated regression accounts for metabolite correlation by using correlated regression errors and can be used to identify biologically significant metabolites [8, 9]. In gene expression data analysis, [10] recently proposed to use confirmatory factor analysis to capture gene-pathway relationships and a conditional autoregressive model to capture relationships between a set of pathways where a network has been constructed based on KEGG [11] pathways. The latter accounts for biological variation in the data and aids interpretation.

In the metabolomics literature, traditional frameworks for metabolomic data analysis use dimensionality reduction techniques, namely principal component analysis (PCA), partial least squares (PLS) [12] and PLS derived models (OPLS [13], O2PLS [14], OnPLS [15]) to take into account high correlations between metabolites. Extension of PLS to O2PLS and OnPLS allows for integrative analysis of heterogeneous omic data. One of the interests of PCA (PLS) derived models is to be able to visually assess whether or not there is a time effect in the data and identify metabolites that change over time by looking into time trajectories of each metabolite [16].

Extensions to PCA and PLS for longitudinal analysis include lagged PCA (PLS) and dynamic PCA (PLS) where a backshift matrix is introduced to take into account time dependency [17, 18]. Similarly, [19] used a set of piecewise orthogonal projections latent structures to describe changes between neighbouring time points. PARAFAC [20] is a multi-linear unsupervised decomposition method that can account for the multi-way variation seen in dynamic metabolomics data. Similar models such as ASCA [21] and APCA [22] also seek to account for temporal variation by combining analysis of variance (ANOVA) with PCA. Recently, dynamic probabilistic PCA (DPPCA) was proposed in [23] as a generative probabilistic model of longitudinal metabolomic data where a stochastic volatility model is used for the latent variables. The main inconvenience of these approaches is that further techniques such as variable importance scores have to be separately applied to the data in order to identify biomarkers and they do not take heterogeneous data (i.e. data from different omics techniques) into account.

In this paper, we are interested in dose-response time course experiments where additional omic variables (bacteria, genes, transcripts, etc) are monitored along with metabolites in the context of biomarker discovery. The main contribution of our work is a single probabilistic generative model that i) can overcome overfitting via the use of weakly informative priors ii) makes use of mixed effects models to model the experimental design iii) models metabolite interactions using pathway information through a conditional auto-regressive (CAR) component and iv) uncovers multiple associations between metabolites and other omic variables by using a horseshoe prior. An additional benefit of our approach is that it naturally yields a list of perturbed metabolic pathways.



We denote a metabolomics data set by \(\boldsymbol {X} \in \mathbb {R}^{N \times T \times M}\) where N is the number of individuals, T the number of time points and M the number of metabolite variables (henceforth termed metabolites for simplicity). \(\boldsymbol {Y} \in \mathbb {R}^{N \times T \times K}\) is an additional continuous omic data measured along with X where K is the number of associated omic variables. The set of N individuals consists of a set of cases and controls. Throughout the paper, index i always runs through individuals, index t runs through time points, index m runs through metabolites and index k runs through Y variables. Vector quantities are written in bold. Matrices are written in bold capitals. Our goal is to build a simple model that can identify biomarkers of a specific treatment over time, taking into account the multiple sources of variation in the data.

The model is built on three levels (see “Methods” section for details on full model): First, a CAR component to capture interaction between metabolites. Second, a variable selection model to uncover associations between metabolites and Y data. Third, a mixed effects component to model the experimental design. This yields the following hierarchical model:

$$\begin{array}{@{}rcl@{}} \boldsymbol{x}_{it}^{e} \vert \boldsymbol{\mu}_{it}, \boldsymbol{C}, \sigma \sim N \left(\boldsymbol{\mu}_{it}, \left(\boldsymbol{I}_{M}-\boldsymbol{C \left(\boldsymbol{\phi}^{e}\right)} \right)^{-1} \sigma^{2} \right) \end{array} $$
$$\begin{array}{@{}rcl@{}} \mu_{itm} = \alpha_{m} + \gamma_{im} + \boldsymbol{\beta}_{m} \boldsymbol{y}_{it} + \nu_{itm} \end{array} $$

where the mean metabolite level μit is a function of covariates of sample i at time point t, influence of metabolite j on metabolite m is captured through coefficients cmj elements of the matrix \(\boldsymbol {C} \left (\boldsymbol {\phi } \right) = \sum _{p=1}^{P}\phi _{p} \boldsymbol {G_{p}}\boldsymbol {A_{p}}\) where P is the number of pathways (see “Methods” section for details). αm represents treatment effect for metabolite m, γim represents individual perturbations for metabolite m, νitm represents temporal effects for metabolite m of individual i at time point t and βm quantifies interactions between metabolite m and other omic variables. The indicator e stands for {controls,cases} and ϕ quantifies pathway perturbation. Particularly, by specifying different dependence parameters for metabolite interactions in cases and controls, the model is able to identify perturbed pathways by comparing ϕcases and ϕcontrols.

We refer to our model as “ iCARH ” model for “ integrative CAR Horseshoe ” model. The model we propose is implemented in the iCARH package and freely available from the Comprehensive R Archive Network (CRAN).

In the following sections, we perform experiments on both synthetic and real data to investigate whether our algorithm gives reasonable solutions. We first try our method on a simulated dataset in “Simulation study” section to get an understanding of the performance of our method. In “Case study” section, we test our method on a dataset comprising metabolomic and bacterial composition profiles in a drug treatment experiment. A fully worked reproducible example using the iCARH package on a publicly available dataset is available in Additional File 1. All experiments were run on a computer with an Intel i7 processor running at 2.8 GHz using 16 Gb of RAM.

Simulation study

To get better understanding of our method and test its applicability, we first perform our approach on synthetic datasets. We will mainly focus on the ability of our model to infer pathway perturbation.

In this simulation, we first extract from the KEGG database the proportion of pathways in which a single metabolite is involved. We then use these proportions to randomly generate a binary matrix Z with dimensions M×P indicating random assignments of metabolites to pathways where M is the number of metabolites, P is the number of pathways. Each design matrix Ap is then equal to zpzpT where zp is the pth column of Z. We simulate pathway perturbation according to an indicator variable ω as follows: If ω=1 then pathways are perturbed hence \(\phi _{p}^{\text {controls}}\) and \(\phi _{p}^{\text {cases}}\) are simulated from normal distributions with different means. If ω=0 then \(\phi _{p}^{\text {controls}}=\phi _{p}^{\text {cases}}\) and pathways are not perturbed.

The rest of the parameters is set as follows : number of bacterial variables K=1, number of metabolites M=40, number of time points T=7, number of samples N=22, number of pathways P=11, global parameter τ fixed to 1.2 (to induce a medium degree of sparsity), parameters νitm,γim,μitm and multi-omics time course data \( \boldsymbol {x}_{it}^{e}\) simulated according to Eqs. 1, 2, 11, and 12, respectively (See “Methods” section). Finally, we generated 10 datasets according to the simulations above in order to assess how our model infers perturbed pathways.

We set non-informative uniform priors on \(\alpha _{m}, \sigma _{\gamma _{im}}, \theta _{m}, \sigma _{\mu _{m}}^{2}\). We set an informative prior on \(\sigma ^{2}_{\gamma _{m}} \sim \text {inverse-gamma} \left (1, 0.1 \right)\) as we expect low variability amongst biological samples of the same group.

Assessing uniform and beta-like priors

We first compare inference of the model under a uniform prior for \(\phi ^{e}_{p}\) and the prior in Eq. 5 (See Fig. 1). Inference is done using 2000 iterations of Hamiltonian Monte Carlo sampling and 1000 warm-up iterations.

Fig. 1

Left : Boxplots of Area Under the Curve (AUC) for pathway perturbation inference with uniform prior on ϕ compared to a beta like prior on ϕ for 10 simulated datasets. We infer perturbations based on the posterior distribution of ϕcontrolsϕcases. Right : Distribution of ϕcontrolsϕcases under the uniform prior (white) and the beta like prior (shaded). The true perturbed pathways are printed in bold on the x axis

We here show, through simulation, that an informative beta like prior compares better than a non-informative uniform prior in inferring significant pathways. The left plot in Fig. 1 shows the boxplots of the 95% confidence interval of the Area Under the Curve (AUC) for pathway perturbation inference for 10 simulated datasets with uniform prior on ϕ and a beta like prior on ϕ. We infer perturbations based on the posterior probability that ϕcontrols and ϕcases are different i.e. the 95% credible interval of ϕcontrolsϕcases does not contain zero. The AUC values for the beta like distribution is significantly higher than the AUC values for the uniform distribution. On average pathway peturbation inference under the uniform distribution reduces to a random guess with an average AUC of 0.53. This is likely due to the lack of variance of the uniform distribution. The right plot in Fig. 1 shows the posterior distributions of ϕcontrolsϕcases under the uniform prior (white) and the beta like prior (shaded) for each pathway. The true perturbed pathways are printed in bold on the x axis. In the following section, we assess pathway inference against design inaccuracies.

Assessing pathway inference against design inaccuracies

It is very common in metabolomics data to find metabolites that are correlated but not in the same KEGG pathway. In the following simulation we assess how inaccuracies in the covariance structure between metabolites and the design matrices Ap affect the iCARH model. We used the 10 datasets from the previous simulation and perturbed the design matrices by selecting a random fraction of metabolites in each pathway. We then randomly (falsely) assign these metabolites to no pathway, or to different pathways. We similarly run the model for 2000 iterations of Hamiltonian Monte Carlo sampling and 1000 warm-up iterations for each of the fractions {0,0.18,0.35,0.44,0.5,0.62} of perturbed metabolites. Finally, in the same fashion, we assess perturbations based on the 95% credible interval of ϕcontrolsϕcases. Figure 2 is a series of average Receiver Operating Characteristic (ROC) curves across 10 datasets for each of the fractions {0,0.18,0.35,0.44,0.5,0.62} of perturbed metabolites. On average, the performance of our model reduces to a random guess (AUC of 0.5) if 50% of the metabolites in each pathway is perturbed. The AUC of our model reaches 0.97 if no metabolites are perturbed and is about 0.88 if 18% of the metabolites in each pathway are perturbed.

Fig. 2

Average Receiver Operating Characteristic (ROC) curves for pathway perturbation inference across 10 datasets for different factions of “falsely” assigned metabolites. Boxplots show the 95% confidence intervals of the true positive rate

Case study

In this section, we test our model on an actual metabolomic data and 16S data for bacterial profiles. In this study we are interested in the influence of a diabetes drug metformin on a non-diabetic model. Metformin is the first-line medicine to treat type 2 diabetes. It has also been suggested that metformin has anti-cancer [24], cardiovascular [25] and anti-aging [26] effects. Because of their very large metabolic capacity, the gut bacteria can influence toxicity and metabolism of drugs. Here, we are particularly looking for metabolic biomarkers indicative of microbiota changes as a result of treatment.

The study design is as follows: metabolic profiles of 24 animals are acquired on 9 equally spaced timepoints using different mass spectrometry techniques from plasma samples. Bacterial profiles are acquired using Illumina MiSeq [27]. The study has allowed for two groups of 12 animals where the drug has been administrated to the second group (timepoints 3 to 7) allowing for an acclimatation period (timepoints 1 and 2) and a recovery period (timepoints 8 and 9). As metabolites are mapped to pathways, prior filtering and/or metabolite annotation needs to be performed beforehand. After data processing and metabolite identification, a total of 56 metabolites and 6 bacteria species are further analyzed using our model. Preliminary investigation of the data shows observable associations between correlations and inter-pathway and intra-pathway metabolites in Fig. 3 which motivates fitting the iCARH model the data. Inference is done using 2000 iterations of Hamiltonian Monte Carlo sampling and 1000 warm-up iterations.

Fig. 3

Figure shows boxplots of correlations between metabolites that are part of the same pathway (intra-pathway correlations) compared to correlations between metabolites that belong to different pathways (inter-pathway correlations) in the metformin study (See “Case study” section). Although correlations are not very close to 1. There is an empirically observed correlation signature related to whether metabolites are in the same pathway or not

We assess performance of our model for different values of τ using the Watanabe-Akaike information criterion (WAIC). Tested values of τ comprise 1, 1.2, 5, 10 with corresponding WAIC values of 7317.296, 7322.798, 7317.457, 7316.476 respectively. WAIC values are very similar for different values of τ which suggests to use the most selective model with τ=10 as it is the simplest i.e with the smallest number of selected variables.

Assessing model fit

In order to assess our model fit, we perform posterior predictive checks of our model compared to DPPCA [23]. The DPPCA model is a multivariate model using PCA, where PCA scores are modelled via a stochastic volatility model. In the Bayesian framework, posterior predictive checks consist in comparing data simulated from the posterior predictive distribution with the observed data. The mean absolute deviations (MADs) are computed between the observed covariance matrix and the covariance matrix of simulated data. The experiment was repeated for different numbers of metabolites selected randomly from the whole dataset. As metabolites are also predictors under the CAR model, the performance is expected to improve when the number of metabolites increases. The process was also replicated for inference using the DPPCA model [23]. Figure 4 shows MADs of our model and the DPPCA model. Although MADs for the DPPCA model decrease when the number of metabolites increases, it is still slightly higher than MADs for the iCARH model. Overall, our model clearly outperforms the DPPCA model.

Fig. 4

Posterior predictive checks for mean absolute deviation (MAD) compared to DPPCA for different numbers of metabolites included. Vertical bars show the 95% confidence intervals of MADs.The MADs decrease as the number of metabolites increases. Our model performs clearly better than the DPPCA model

In addition to posterior predictive checks, normality checks are another way to assess if the observed results are not mainly a product of misspecified priors. Specifically, goodness of fit was checked by using \(\Psi ^{-1}_{e} \left (\boldsymbol {x}_{it} - \boldsymbol {\mu }_{it} \right) \sim N\left (0, \boldsymbol {I}_{M}\right)\) where Ψe denotes the Cholesky factor of (IMC(ϕe))−1σ2. Zero-mean and normality were thus checked for \(\Psi ^{-1}_{e} \left (\boldsymbol {x}_{it} - \boldsymbol {\mu }_{it} \right)\) (See Fig. 5).

Fig. 5

Right and left panels show model fit assessment for controls and cases for metformin data. Left : quantile-quantile normal plot of \(\Psi ^{-1}_{\text {cases}} \left (\boldsymbol {x}_{it} - \boldsymbol {\mu }_{it} \right)\). Right : quantile-quantile normal plot of \(\Psi ^{-1}_{\text {controls}} \left (\boldsymbol {x}_{it} - \boldsymbol {\mu }_{it} \right)\)

Data results

In the standard model, αm represents an indicator variable for the treatment effect. The treatment variable can also be continuous as in this data example (drug measurements) and modeled by \(\boldsymbol {\alpha }_{m} = \beta ^{\alpha }_{m} \boldsymbol {y}_{\text {drug}}\). The treatment effect can now be simply summarized by \(\beta ^{\alpha }_{m}\). Figure 6 is a series of boxplots of 95% credible intervals of posterior means of \(\beta ^{\alpha }_{m}\) for metabolites 13 to 31. We are mainly interested in “metabolite 27” as it is associated with bacteria species 2.

Fig. 6

Estimates of effects of treatment on metabolite profiles are captured by \(\beta ^{\alpha }_{m}\) for the metformin data. The figure depicts boxplots of 95% credible intervals of posterior means of \(\beta ^{\alpha }_{m}\). Only part of the data is plotted as we are mainly interested in “metabolite 27”

Figures 7 and 8 show posterior distributions of ϕe for each pathway and estimates of effects of bacteria on metabolites. Results in “Simulation study” section suggest to compare the covariance structure of metabolites in the observed data with the covariance induced by the design matrices in order to have an a priori idea on the robustness of pathway inference (See Fig. 2). For a correlation threshold of 0.3, about 25% of the metabolites are misspecified in the design matrices which corresponds to an AUC around 0.8 according to Fig. 2. If we set a higher correlation threshold, a lower number of metabolites are misspecified. For example, for a correlation threshold of 0.5, only 8% of the metabolites are misspecified. This supports the use of the iCARH model for pathway perturbation inference for this data.

Fig. 7

Posterior distributions of ϕe for each pathway for each treatment group for the metformin dataset. Posterior distributions of ϕe for some pathways are flatter for the controls than the cases which might be indicative mild pathway alterations. These pathways correspond as well to the top 11 pathways inferred by MetaboAnalyst [30]. The p-values are the p-values returned by MetaboAnalyst for each pathway

Fig. 8

Metformin data: Estimates of effects of bacteria on metabolite profiles are captured by βm. The figure depicts boxplots of 95% credible intervals of posterior means of βm. Different metabolites present significant changes along with some bacterial profiles. For example, “metabolite 27” is positively associated with bacteria species 2 but negatively associated with bacteria species 5 and can be considered as indicator of changes in these bacteria species abundance

Estimates of effects of bacteria on metabolite profiles are captured by βm. Some metabolites present significant changes along with the bacterial profiles. For example, “metabolite 27”, a hydroxy fatty acid, is associated with alterations in abundance of 4 bacteria species. Figure 7 shows that, as a result of treatment, KEGG pathways are not significantly altered. However, distributions of ϕcontrols for “fatty acids biosynthesis” and “biosynthesis of unsaturated fatty acids” KEGG pathways are remarkably flatter than the distributions of ϕcases. These pathways involve the previously identified hydroxy fatty acid metabolite. Our analysis confirms previously reported studies that hydroxy fatty acids might be produced by the gut microbiome [28, 29]. On the other hand, results from MetaboAnalyst [30] give p-values between 6.3×10−2 and 1.8×10−10 indicating that all pathways are significant to changes of the treatment except one; glycerophospholipid metabolism. We think that this discrepancy in the results is due to the way iCARH and metaboAnalyst model pathway perturbation i.e. whilst metaboAnalyst considers that mean level changes of one or more metabolite concentrations involved in a pathway indicate perturbation of the latter pathway. iCARH considers that changes in the covariance structure of metabolites in the same pathway are indicative of pathway perturbation.


Identifying biomarkers in time course metabolic data and inferring significant associations with heterogeneous omic variables is extremely challenging due to the several sources of variations of the data. In addition, existing methods developed to analyse such data are very scarce and have the limitations of i) overfitting to the few available data points or ii) confounding the experimental and longitudinal variation or iii) ignoring the metabolite interactions or iv) ignoring effects of other omic variables. In this paper, the model we have developed combines several approaches to take into account the different aspects of the data namely the number of time points, the experimental variation captured by μit, interactions between metabolites captured by ϕ and interactions with additional omic variables captured by βm.

Our results demonstrate that our model successfully addresses the main questions of a metabolomic study. Most importantly, our model is able to identify metabolic biomarkers related to treatment, infer perturbed pathways as a result of treatment and find significant associations with additional omic variables. We have shown that providing an informative prior on metabolic pathways and an informative prior over the parameter ϕ is a significant improvement over the DPPCA model. Particularly, our model is more robust to slight variations usually observed in short time series data thanks to the small number of covariance parameters (in the covariance matrix) needed to estimate compared to DPPCA. We have also shown through simulation that an informative beta like prior compares better than a non-informative uniform prior in inferring significant pathways. On different real data, we have investigated how the number of profiled metabolites can affect the predictive ability of the model and carried out a fully reproducible application of iCARH.

Several potential extensions arise naturally from our model. In terms of the metabolite interactions component, many research questions can arise. Alternative strategies to modeling metabolite interactions can be examined such as modeling the non-zero elements of the adjacency matrix C of each pathway as random variables. This strategy was adopted in the CAR literature by [31, 32] to take into account step changes in spatial variation. Step changes can potentially be useful to model changes in metabolites correlations as a result of treatment. Lee [33] provide an overview of different CAR models used in spatial modeling. The proposed models can be adapted to fit into the metabolomics literature.

Another potential extension concerns the source for pathway annotation and modeling. In this research paper, the KEGG pathways are used but this poses certain limitations regarding the performance of the model, due to some shortcomings in the database (missing compounds, inaccuracies in the database etc). In this sense extending the tool to support software formats and applications that enable assembly of superpathways (SBML [34], OWL [35], KEGGConverter [36] or KENEV [37]) would probably increase its performance in terms of accuracy and computational effeciency, as the adjacency matrix would be known a priori and its sparsity would be stronger.

From a practical point of view, the model has been fitted using HMC sampling but takes a large amount of time (about 1 h) mostly because of the variable selection procedure and metabolites interdependence. This could be addressed by using variational Bayes. In fact, variational Bayes inference procedures offer cost-effective inference by means of principled approximations and appealing computational time for high dimensional data. A variational bayes inference of CAR models was proposed by [38] for high dimensional data, and a variational bayes approach for variable selection was recently proposed by [39].


Metabolomics longitudinal profiling techniques are imperative to understand the effect of a drug or a disease across time and can provide enhanced understanding of the underlying biology of the system. In a data integration framework, we have illustrated the use of the CAR model to incorporate metabolites interactions in the model and the horseshoe prior to identify association with heterogeneous omic variables obtained by other omic techniques. The combination of the CAR and horseshoe levels yields the “integrative CAR Horseshoe” (iCARH) model which we presented in this article. Our model is accompanied by an R package with various visualization functions easy-to-use for applied researchers.

The iCARH model has various appealing features such that it is able to identify metabolic biomarkers related to treatment, infer perturbed pathways as a result of treatment and identify potential associations between heterogeneous omic variables. Clearly, these appealing features open up further research topics.


In this section we describe theoretical details behind the three levels of our iCARH model: Metabolite dependencies, integrative analysis with other omics data and experimental design.

Metabolite dependencies

In any integrative biological model, it is useful to be able to interpret the model at a systems level, e.g. according to functional groups of biological molecules, rather than attempting to interpret results for individual molecules. Metabolic pathways are the most widely used groupings for this type of analysis in metabolomics, and have been widely used to interpret experimental data, usually by performing over-representation or enrichment analyses [4044]. Since pathways are regulated in a coordinated fashion, it is natural to assume that the levels of metabolites which are members of the same pathway may be correlated. This dependence, though weak, is observable in associations between correlations and network distance [4547], and also observable in the real data used in our study (See Fig. 3). We therefore incorporate a pathway-based correlation component into our model via a CAR approach, in a similar fashion to [10]. The extent of pathway specific correlations will vary according to the experimental system and assay, and may in some cases provide little extra information. Nonetheless, including such a pathway based component can greatly increase the interpretability of the resulting model beyond one which does not include such grouping information. In this context, some metabolite peaks in the data need to be identified in order to be mapped to pathways. As it will be later clear in the CAR model we will use, if metabolites are not identified and hence can not be mapped to pathways, no pathway-induced correlation will be assumed.

We assume that the concentration of each metabolite is linearly influenced by concentration levels of metabolites in the same pathway. Linear dependencies have been investigated in genomics in order to uncover functional modules [45]. Modeling linear associations is appealing as it captures the overall trend and also less prone to overfitting small amounts of data. Let \(\boldsymbol {C} \in \mathbb {R}^{M \times M}\) be the design matrix quantifying metabolite interactions such that matrix elements cmm=0,cmj≠0 if metabolites m and j are in the same pathway and 0 otherwise. Thus, metabolite levels can be expressed as:

$$ {}x_{itm} \vert \boldsymbol{x}_{it,-m}, \boldsymbol{\mu}_{it}, \boldsymbol{C}, \sigma \sim N\left(\mu_{itm} + \sum_{\substack{j=1\\ j \neq m}}^{M} c_{mj} (x_{itj}-\mu_{itj}), \sigma^{2}\right) $$

where xit,−m represents measurements of metabolites of sample i at time point t excluding metabolite m, and μitm is a function of covariates of sample i for metabolite m at time point t taking into account additional variation in the data (See “Integrative analysis” and “Experimental design” sections). If we define IM the Mth order identity matrix, the joint distribution of xit can be explicitly written as [48]:

$$\begin{array}{@{}rcl@{}} \boldsymbol{x}_{it} \vert \boldsymbol{\mu}_{it}, \boldsymbol{C}, \sigma \sim N \left(\boldsymbol{\mu}_{it}, \left(\boldsymbol{I}_{M}-\boldsymbol{C} \right)^{-1} \sigma^{2} \right) \end{array} $$

An important output of our modeling procedure is identification of which pathways are “on" or “off" as an effect of treatment. In the CAR literature, the design matrix C can be modeled as a scaled product of a diagonal weight matrix and an adjacency matrix. In order to infer which pathways are perturbed we construct the distance matrix based on the individual contribution of each pathway. To be precise, we define \(\boldsymbol {C} \left (\boldsymbol {\phi } \right) = \sum _{p=1}^{P}\phi _{p} \boldsymbol {G_{p}}\boldsymbol {A_{p}}\) where P is the number of pathways. The distance matrices Ap are a zero-diagonal symmetric adjacency matrices with elements \(a_{mj}^{p}\) equal to the inverse of the length of the shortest path between metabolites m and j if they are in pathway p and 0 otherwise. A path between two metabolites consists of the number of reactions that lead from one metabolite to the other, and the shortest path is the path that contains the smallest number of reactions. The diagonal matrices Gp comprise the reciprocal of the number of neighbors of each metabolite in pathway p i.e \( \left (g_{mm}^{p} \right)^{-1}= \sum _{j=1}^{M} (a_{mj}>0) \) so that the squared partial correlation between two metabolites \(\text {cor} \left (x_{itm}, x_{itj} \vert \boldsymbol {x}_{it,-(m,j)} \right)^{2} \propto \phi _{p}^{2} g_{mm}^{p} g_{jj}^{p} \) is reduced when more metabolites from the same pathway are profiled [48]. The vector of coefficients \(\boldsymbol {\phi } = \lbrace \phi _{p} \rbrace _{p=1}^{P}\) is estimated from the data. It is referred to as spatial-dependence parameter in the CAR literature. In the context of this work, the vector of coefficients \(\boldsymbol {\phi } = \lbrace \phi _{p} \rbrace _{p=1}^{P}\) quantifies pathway contribution, for example ϕ1=0 indicates no contribution.

Under the CAR setting, we turn the reader attention that if pathway information is not available (i.e. all/some metabolites are not identified) then no pathway-induced correlation is assumed in the data and inference will be performed such that the design matrix C in Eq. 4 is a zero matrix. Hence, metabolites are assumed to be independent as the covariance matrix between metabolites is diagonal in this case.

The model needs to comply with the condition that IMC(ϕ) is positive definite. If we assume that pathways are a priori equally perturbed, ϕp must fall in the interval \(\left (\frac {1}{P\xi _{p}^{1}}, \frac {1}{P\xi _{p}^{2}} \right)\) where \(\xi _{p}^{1}\) and \(\xi _{p}^{2}\) are the minimum and maximum eigenvalues of GpAp, respectively. In practice, strong interaction between observed metabolites of pathway p is reproduced in CAR models only when the scaling parameter ϕp is quite close to one of the boundaries \(\frac {1}{P\xi _{p}^{1}}, \frac {1}{P\xi _{p}^{2}}\). Hence, we use a beta-type prior for ϕp that places substantial mass on large values of |ϕp| [49]:

$$\begin{array}{@{}rcl@{}} \mathrm{p} \left(\phi_{p} \right) = \frac{1}{B \left(\frac{1}{2}, \frac{1}{2} \right)} \left(\phi_{p}- \frac{1}{P\xi_{p}^{1}} \right)^{-\frac{1}{2}} \left(\frac{1}{P\xi_{p}^{2}} - \phi_{p} \right)^{-\frac{1}{2}} \end{array} $$

where B is the beta function. The parameter σ2 captures variance heterogeneity in metabolite intensities and is given an inverse gamma prior G(ψ,ψ−1). This prior provides 2ψ pseudo-observations in addition to NT available observations. In order to build a reasonably informative prior we set ψ=N×T/4.

Integrative analysis

In this section, we turn our attention to modeling the association between heterogeneous omic variables such as transcripts and metabolites. Association between omic variables involves complex processes where often only few variables are significant which motivates the use of shrinkage priors for integrative analysis and cross-omics biomarker discovery (See [50] for a review on shrinkage priors). Recently, [51] proposed the “horseshoe” prior as a prior based on a scale mixture of normals where scale parameters are modeled as the product of a global shrinkage (scale) parameter and a local shrinkage (scale) parameter. This definition allows for an additional flexibility where sparsity can be controlled at a global level for each metabolite (i.e. how many non-zero coefficients?) and a local level for each metabolite (i.e. which coefficients are non-zero?). The horseshoe prior has been widely recognized and extended by the statistical community since its introduction by [51] as it benefits from various desirable properties such as simple analytic form, easy computation and preservation of significant coefficients (no over-shrinkage) [52]. In order to model the association between heterogeneous omic variables, we extend the horseshoe prior via the following hierarchical shrinkage model by introducing an additional variable τ to control the overall sparsity level for all metabolites:

$$\begin{array}{@{}rcl@{}} \mu_{itm} = \alpha_{m} + \gamma_{im} + \boldsymbol{\beta}_{m} \boldsymbol{y}_{it} + \nu_{itm} \end{array} $$
$$\begin{array}{@{}rcl@{}} \beta_{mk} \vert \lambda_{mk}, \sigma_{\beta_{m}} \sim N \left(0, \lambda_{mk}^{2} \sigma^{2}_{\beta_{m}}\right) \end{array} $$
$$\begin{array}{@{}rcl@{}} \lambda_{mk} \vert \tau \sim \text{St}^{+} \left(\tau,0,1 \right) \end{array} $$

where αm represents the treatment effect for metabolite m, \( \gamma _{im} \sim N\left (0, \sigma _{\gamma _{m}}^{2} \right)\) represents individual perturbations for metabolite m, \(\nu _{itm} \vert \nu _{i,t-1,m} \sim N\left (\theta _{m} \nu _{i,t-1,m}, \sigma _{\nu _{m}}^{2}\right)\) follows and auto-regressive process and represents temporal effects for metabolite m of individual i at time point t. βm is a vector of dimension K that quantifies interactions between metabolite m and other omic variables encoded in the vector yit of dimension K. λmk is called the local shrinkage parameter whilst \(\sigma ^{2}_{\beta _{m}}\) is the global shrinkage parameter. St+ denotes the half Student-t distribution with τ degrees of freedom. For τ=1, this prior reduces to the horseshoe prior [51]. Intuitively, for small values of λmk the coefficient βmk is very close to 0 while for relevant variables λmk will be large. In addition, \(\sigma _{\beta _{m}}\) controls the overall shrinkage level i.e sparsity of the vector βm is more important for small values of \(\sigma _{\beta _{m}}\).

Define \(\kappa _{mk} = \frac {1}{1+ \lambda _{mk}^{2} \sigma ^{2}_{\beta _{m}}/ \tau }\) a random shrinkage coefficient such that κkm≈0 when λmk is large and κkm≈1 when λmk is small. This transformation implies the following prior distribution on κmk:

$$ {}\mathrm{p} \left(\kappa_{mk} \vert \tau, \sigma_{\beta_{m}} \right) = \frac{1}{2 \sqrt{\pi} B \left(\frac{\tau}{2}, \frac{1}{2} \right)} \frac{\sigma_{\beta_{m}}^{\tau} \kappa_{mk}^{\tau/2-1} \left(1-\kappa_{mk}\right)^{-1/2}}{\left(1- \kappa_{mk} + \kappa_{mk} \sigma_{\beta_{m}}^{2} \right)} $$

This prior density is shown in Fig. 9 for different values of \(\sigma _{\beta _{m}}\) and τ. It reduces to a Beta (τ/2,1/2) distribution if \(\sigma _{\beta _{m}}=1\) and to a Beta (1/2,1/2) which looks like a horseshoe, if in addition τ=1. When τ increases, Beta (τ/2,1/2) skews towards 1 which increases the global shrinkage power. The expectation of βm given Y,κm,τ,μtm can be expressed as:

$$\begin{array}{@{}rcl@{}} \mathbb{E} \left(\boldsymbol{\beta}_{m} \vert \boldsymbol{Y}, \boldsymbol{\kappa}_{m}, \tau, \boldsymbol{\mu}_{tm} \right) &=& \left(\sum_{t=1}^{T} \boldsymbol{Y}_{t}^{T} \Sigma_{m}^{-1} \boldsymbol{Y}_{t} + \tau \Upsilon_{m} \right)^{-1} \\ & & \times \sum_{t=1}^{T} \boldsymbol{Y}_{t}^{T} \Sigma_{m}^{-1} \boldsymbol{\mu}_{tm} \end{array} $$
Fig. 9

Shrinkage prior \(\mathrm {p} \left (\kappa _{mk} \vert \tau, \sigma _{\beta _{m}} \right)\) on κmk for different values of \(\sigma _{\beta _{m}}\) and τ. The prior distribution skews towards 1 if τ increases or \(\sigma _{\beta _{m}}\) decreases (shrinkage). It skews towards 0 if τ decreases or \(\sigma _{\beta _{m}}\) increases (no shrinkage)

where \(\Sigma _{m} = \left (\frac {\sigma _{\nu _{m}}^{2}}{1-\theta _{m}^{2}} + \sigma _{\gamma _{m}}^{2} \right) \boldsymbol {I}_{N}\) and Υm is a diagonal matrix of order K with elements 1/κmk−1. Equation (10) introduces a penalty term τΥm where Υm is a metabolite specific penalty term introduced by the horseshoe prior and τ is a global penalty term. Precisely, τ captures the overall sparsity level amongst all metabolites. The expectation of βm given Y,κm,τ,μtm is very similar to the estimate of βm under ridge regression where τΥm simply reduces to τIN.

The global sparsity level can be controlled using τ. Increasing the global sparsity level is a desired property in omic studies, as usually we deal with a large number of omic variables where only few are important. In appendix B we discuss how τ can be fixed a priori.

Fig. 10

Plates diagram of the iCARH model. Fixed variables are represented by squares, random variables by circles and observations are shaded. For clarity, all variances \(\sigma ^{2}, \sigma _{\beta _{m}}^{2}, \sigma _{\gamma _{m}}^{2}, \sigma _{\nu _{m}}^{2}\) are not represented in the diagram

Experimental design

The covariance structure might change drastically as a result of treatment if the latter affects relationships between metabolites. The model can be extended to take into account the experimental design. As specified in the previous section, αm captures the treatment effect for metabolite m, γim represents individual perturbations for metabolite m, \(\nu _{itm} \vert \nu _{i,t-1,m} \sim N\left (\theta _{m} \nu _{i,t-1,m}, \sigma _{\nu _{m}}^{2}\right)\) represents temporal effects for metabolite m of individual i at time point t in Eq. 7. In addition, we allow covariance structures C(ϕe) to be different for the control samples and the cases where e{cases, controls} designates experimental groups. This yields the overall hierarchical model:

$$\begin{array}{@{}rcl@{}} \boldsymbol{x}_{it}^{e} \vert \boldsymbol{\mu}_{it}, \boldsymbol{C}, \sigma \sim N \left(\boldsymbol{\mu}_{it}, \left(\boldsymbol{I}_{M}-\boldsymbol{C \left(\boldsymbol{\phi}^{e}\right)} \right)^{-1} \sigma^{2} \right) \\ \mu_{itm} = \alpha_{m} + \gamma_{im} + \boldsymbol{\beta}_{m} \boldsymbol{y}_{it} + \nu_{itm} \\ \beta_{mk} \vert \lambda_{mk}, \sigma_{\beta_{m}} \sim N \left(0, \lambda_{mk}^{2} \sigma^{2}_{\beta_{m}} \right) \\ \lambda_{mk} \vert \tau \sim \text{St}^{+} \left(\tau,0,1 \right) \\ \gamma_{im} \vert \sigma_{\gamma_{m}} \sim N\left(0, \sigma_{\gamma_{m}}^{2} \right) \end{array} $$
$$\begin{array}{@{}rcl@{}} \nu_{itm} \vert \theta_{m}, \sigma_{\nu_{m}} \sim N\left(\theta_{m} \nu_{i,t-1,m}, \sigma_{\nu_{m}}^{2}\right) \end{array} $$

A key point of this model is that by specifying different dependence parameters for metabolite interactions in cases and controls, the model is able to identify perturbed pathways by comparing ϕcases and ϕcontrols.

Appendix A: Global sparsity

When there is prior knowledge available, specifying τ a priori can optimize the inference and additionally, provide a more informative prior on λmk. If we fix \(\mathrm {p}\left (\sigma _{\beta _{m}}^{2}\right) \propto 1/\sigma _{\beta _{m}}^{2}\), integrating over \(\sigma _{\beta _{m}}\) gives the expected value of κmk as :

$${}{\begin{aligned} \mathbb{E} \left(\kappa_{km} \vert \tau \right) = \frac{\Gamma \left(1/2\right)^{-1}}{2 \sqrt{\pi} \Gamma \left(\tau/2\right)} G_{3,3}^{2,3} \left(\begin{array}{ll} 1, \tau/2, 0 \\ \tau/2, \tau/2-1/2, 0 \end{array} \, \middle\vert \, 1-\sigma_{\beta_{m}}^{2} \right) \end{aligned}} $$

where \( G_{\cdot,\cdot }^{\cdot,\cdot } \) is Meijer’s G-function [53]. The equation above can be used to fix τ a priori by defining the expected proportion of shrunk coefficients. In practice, different values of τ are plugged into the equation above to get the desired proportion of shrunk coefficients. However, many definite integrals can be obtained using the tables of Meijer functions in [54] for special values of parameters.

Table 1 iCARH model summary

Appendix B: Model summary

Table 1 depicts a summary of model parameters specifying parameters of interest, other inferred parameters and user specified parameters. Figure 10 shows the plates diagram of the iCARH model where fixed variables are represented by squares, random variables by circles and observations are shaded. For clarity, all variances \(\sigma ^{2}, \sigma _{\beta _{m}}^{2}, \sigma _{\gamma _{m}}^{2}, \sigma _{\nu _{m}}^{2}\) are not represented in the diagram.

The choice of the gamma distribution for \( \sigma _{\beta _{m}}^{2}, \sigma _{\gamma _{m}}^{2}, \sigma _{\nu _{m}}^{2}\) follows the same principle used in “Model” section for σ2. For each variance parameter, the gamma prior provides half pseudo-observations in addition to the available observations e.g \(\sigma _{\nu _{m}}^{2}\) has a G(T/4,T/4−1) prior such that it provides T/2 pseudo-observations in addition to T observations so that the prior is reasonably informative.

Availability of data and materials

The proposed method has been implemented in the R package iCARH that is available from CRAN. A fully worked example with publicly available data is in Additional file 1.



Conditional auto-regressive model


Dynamic probabilistic PCA


Integrative conditional auto-regressive horseshoe model


Principal component analysis


Partial least squares


Two-way orthogonal partial least squares


  1. 1

    Joyce AR, Palsson BØ. The model organism as a system: integrating’omics’ data sets. Nat Rev Mol Cell Biol. 2006; 7(3):198–210.

    CAS  PubMed  Google Scholar 

  2. 2

    Ebrahim A, Brunk E, Tan J, O’brien EJ, Kim D, Szubin R, Lerman JA, Lechner A, Sastry A, Bordbar A, et al. Multi-omic data integration enables discovery of hidden biological regularities. Nat Commun. 2016; 7.

  3. 3

    Daviss B. Growing pains for metabolomics: the newest’omic science is producing results–and more data than researchers know what to do with. The Scientist. 2005; 19(8):25–29.

    Google Scholar 

  4. 4

    Reo NV. Nmr-based metabolomics. Drug Chem Toxicol. 2002; 25(4):375–82.

    CAS  PubMed  Google Scholar 

  5. 5

    Dettmer K, Aronov PA, Hammock BD. Mass spectrometry-based metabolomics. Mass Spectrom Rev. 2007; 26(1):51–78.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. 6

    Berk M, Ebbels T, Montana G. A statistical framework for biomarker discovery in metabolomic time course data. Bioinformatics. 2011; 27(14):1979–85.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. 7

    Mei Y, Kim SB, Tsui KL. Linear-mixed effects models for feature selection in high-dimensional NMR spectra. Expert Syst Appl. 2009; 36(3 PART 1):4703–8.

    Google Scholar 

  8. 8

    Chen C, Deng L, Wei S, Nagana Gowda GA, Gu H, Chiorean EG, Abu Zaid M, Harrison ML, Pekny JF, Loehrer PJ, Zhang D, Zhang M, Raftery D. Exploring metabolic profile differences between colorectal polyp patients and controls using seemingly unrelated regression. J Proteome Res. 2015; 14(6):2492–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  9. 9

    Chen C, Nagana Gowda GA, Zhu J, Deng L, Gu H, Chiorean EG, Abu Zaid M, Harrison M, Zhang D, Zhang M, Raftery D. Altered metabolite levels and correlations in patients with colorectal cancer and polyps detected using seemingly unrelated regression analysis. Metabolomics. 2017; 13(11):125.

    PubMed  PubMed Central  Google Scholar 

  10. 10

    Pham LM, Carvalho L, Schaus S, Kolaczyk ED. Perturbation Detection Through Modeling of Gene Expression on a Latent Biological Pathway Network: A Bayesian hierarchical approach. J Am Stat Assoc. 2015; 1459(July 2016):1–61.,

    Google Scholar 

  11. 11

    Kanehisa M, Goto S. Kegg: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28(1):27–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. 12

    Wold S, Martens H, Wold H. The multivariate calibration problem in chemistry solved by the pls method. Matrix Pencils. 1983:286–93.

    Google Scholar 

  13. 13

    Trygg J, Wold S. Orthogonal projections to latent structures (o-pls). J Chemometr. 2002; 16(3):119–28.

    CAS  Google Scholar 

  14. 14

    Trygg J, Wold S. O2-pls, a two-block (x–y) latent variable regression (lvr) method with an integral osc filter. J Chemometr. 2003; 17(1):53–64.

    CAS  Google Scholar 

  15. 15

    Löfstedt T, Trygg J. Onpls—a novel multiblock method for the modelling of predictive and orthogonal variation. J Chemometr. 2011; 25(8):441–55.

    Google Scholar 

  16. 16

    Antti H, Bollard M, Ebbels T, Keun H, Lindon J, Nicholson J, Holmes E. Batch statistical processing of 1h nmr-derived urinary spectral data. J Chemometr. 2002; 16(8-10):461–8.

    CAS  Google Scholar 

  17. 17

    Kaspar MH, Ray WH. Dynamic pls modelling for process control. Chem Eng Sci. 1993; 48(20):3447–61.

    CAS  Google Scholar 

  18. 18

    Ku W, Storer RH, Georgakis C. Disturbance detection and isolation by dynamic principal component analysis. Chemom Intell Lab Syst. 1995; 30(1):179–96.

    CAS  Google Scholar 

  19. 19

    Rantalainen M, Cloarec O, Ebbels TM, Lundstedt T, Nicholson JK, Holmes E, Trygg J. Piecewise multivariate modelling of sequential metabolic profiling data. BMC Bioinformatics. 2008; 9(1):105.

    PubMed  PubMed Central  Google Scholar 

  20. 20

    Bro R. Parafac. tutorial and applications. Chemom Intell Lab Syst. 1997; 38(2):149–71.

    CAS  Google Scholar 

  21. 21

    Smilde AK, Jansen JJ, Hoefsloot HC, Lamers R-JA, Van Der Greef J, Timmerman ME. Anova-simultaneous component analysis (asca): a new tool for analyzing designed metabolomics data. Bioinformatics. 2005; 21(13):3043–8.

    CAS  PubMed  Google Scholar 

  22. 22

    Harrington P. d. B., Vieira NE, Espinoza J, Nien JK, Romero R, Yergey AL. Analysis of variance–principal component analysis: A soft tool for proteomic discovery. Analytica chimica acta. 2005; 544(1-2):118–127.

    CAS  Google Scholar 

  23. 23

    Nyamundanda G, Gormley IC, Brennan L. A dynamic probabilistic principal components model for the analysis of longitudinal metabolomics data. J R Stat Soc Ser C (Appl Stat). 2014; 63(5):763–82.

    Google Scholar 

  24. 24

    Sahra IB, Le Marchand-Brustel Y, Tanti J-F, Bost F. Metformin in cancer therapy: a new perspective for an old antidiabetic drug?Mol Cancer Ther. 2010; 9(5):1092–9.

    PubMed  Google Scholar 

  25. 25

    Group UPDSU, et al. Effect of intensive blood-glucose control with metformin on complications in overweight patients with type 2 diabetes (ukpds 34). The Lancet. 1998; 352(9131):854–65.

    Google Scholar 

  26. 26

    Anisimov VN, Berstein LM, Egormin PA, Piskunova TS, Popovich IG, Zabezhinski MA, Tyndyk ML, Yurova MV, Kovalenko IG, Poroshina TE, et al. Metformin slows down aging and extends life span of female shr mice. Cell Cycle. 2008; 7(17):2769–73.

    CAS  PubMed  Google Scholar 

  27. 27

    Rapin A, Pattaroni C, Marsland BJ, Harris NL. Microbiota analysis using an illumina miseq platform to sequence 16s rrna genes. Curr Protoc Mouse Biol. 2017:100–29.

    CAS  PubMed  Google Scholar 

  28. 28

    Kishino S, Takeuchi M, Park S-B, Hirata A, Kitamura N, Kunisawa J, Kiyono H, Iwamoto R, Isobe Y, Arita M, et al. Polyunsaturated fatty acid saturation by gut lactic acid bacteria affecting host lipid composition. Proc Natl Acad Sci. 2013; 110(44):17808–13.

    CAS  PubMed  Google Scholar 

  29. 29

    Kimura I, Ozawa K, Inoue D, Imamura T, Kimura K, Maeda T, Terasawa K, Kashihara D, Hirano K, Tani T, et al. The gut microbiota suppresses insulin-mediated fat accumulation via the short-chain fatty acid receptor gpr43. Nat Commun. 2013; 4:1829.

    PubMed  PubMed Central  Google Scholar 

  30. 30

    Chong J, Soufan O, Li C, Caraus I, Li S, Bourque G, Wishart DS, Xia J. Metaboanalyst 4.0: towards more transparent and integrative metabolomics analysis. Nucleic Acids Res. 2018; 46(W1):486–94.

    Google Scholar 

  31. 31

    Lee D, Mitchell R. Locally adaptive spatial smoothing using conditional auto-regressive models. J R Stat Soc Ser C (Appl Stat). 2013; 62(4):593–608.

    Google Scholar 

  32. 32

    Rushworth A, Lee D, Sarran C. An adaptive spatiotemporal smoothing model for estimating trends and step changes in disease risk. J R Stat Soc Ser C (Appl Stat). 2017; 66(1):141–57.,

    Google Scholar 

  33. 33

    Lee D. A comparison of conditional autoregressive models used in bayesian disease mapping. Spat Spatio-temporal Epidemiol. 2011; 2(2):79–89.

    Google Scholar 

  34. 34

    Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A, et al. The systems biology markup language (sbml): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003; 19(4):524–31.

    CAS  PubMed  Google Scholar 

  35. 35

    Bleasby A, Akrigg D, Attwood T. Owl–a non-redundant composite protein sequence database. Nucleic Acids Res. 1994; 22(17):3574.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. 36

    Moutselos K, Kanaris I, Chatziioannou A, Maglogiannis I, Kolisis FN. Keggconverter: a tool for the in-silico modelling of metabolic networks of the kegg pathways database. BMC Bioinformatics. 2009; 10(1):324.

    PubMed  PubMed Central  Google Scholar 

  37. 37

    Pilalis E, Koutsandreas T, Valavanis I, Athanasiadis E, Spyrou G, Chatziioannou A. Kenev: A web-application for the automated reconstruction and visualization of the enriched metabolic and signaling super-pathways deriving from genomic experiments. Comput Struct Biotechnol J. 2015; 13:248–255.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. 38

    Harrison LM, Green GG. A bayesian spatiotemporal model for very large data sets. NeuroImage. 2010; 50(3):1126–1141.

    CAS  PubMed  Google Scholar 

  39. 39

    Ormerod JT, You C, Müller S, et al. A variational bayes approach to variable selection. Electr J Stat. 2017; 11(2):3549–94.

    Google Scholar 

  40. 40

    Xia J, Wishart DS. Metpa: a web-based metabolomics tool for pathway analysis and visualization. Bioinformatics. 2010; 26(18):2342–4.

    CAS  PubMed  Google Scholar 

  41. 41

    Kamburov A, Cavill R, Ebbels TM, Herwig R, Keun HC. Integrated pathway-level analysis of transcriptomics and metabolomics data with impala. Bioinformatics. 2011; 27(20):2917–8.

    CAS  PubMed  Google Scholar 

  42. 42

    Chagoyen M, Pazos F. Tools for the functional interpretation of metabolomic experiments. Brief Bioinforma. 2012; 14(6):737–44.

    Google Scholar 

  43. 43

    Kankainen M, Gopalacharyulu P, Holm L, Orešič M. Mpea—metabolite pathway enrichment analysis. Bioinformatics. 2011; 27(13):1878–9.

    CAS  PubMed  Google Scholar 

  44. 44

    Gao J, Tarcea VG, Karnovsky A, Mirel BR, Weymouth TE, Beecher CW, Cavalcoli JD, Athey BD, Omenn GS, Burant CF, et al. Metscape: a cytoscape plug-in for visualizing and interpreting metabolomic data in the context of human metabolic networks. Bioinformatics. 2010; 26(7):971–3.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. 45

    Walther D, Strassburg K, Durek P, Kopka J. Metabolic pathway relationships revealed by an integrative analysis of the transcriptional and metabolic temperature stress-response dynamics in yeast. Omics J Integr Biol. 2010; 14(3):261–74.

    CAS  Google Scholar 

  46. 46

    Gipson GT, Tatsuoka KS, Sokhansanj BA, Ball RJ, Connor SC. Assignment of ms-based metabolomic datasets via compound interaction pair mapping. Metabolomics. 2008; 4(1):94–103.

    CAS  Google Scholar 

  47. 47

    Krumsiek J, Suhre K, Illig T, Adamski J, Theis FJ. Gaussian graphical modeling reconstructs pathway reactions from high-throughput metabolomics data. BMC Syst Biol. 2011; 5(1):21.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. 48

    Cressie N, Wikle CK. Statistics for Spatio-temporal Data: Wiley; 2015.

    Google Scholar 

  49. 49

    Banerjee S, Carlin BP, Gelfand AE. Hierarchical Modeling and Analysis for Spatial Data: Crc Press; 2014.

  50. 50

    Richardson S, Bottolo L, Rosenthal JS. Bayesian models for sparse regression analysis of high dimensional data. Bayesian Stat. 2010; 9:539–69.

    Google Scholar 

  51. 51

    Carvalho CM, Polson NG, Scott JG. Handling sparsity via the horseshoe. In: Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics. PMLR: Florida: 2009. p. 73–80.

    Google Scholar 

  52. 52

    Van Der Pas S, Kleijn B, Van Der Vaart A, et al. The horseshoe estimator: Posterior concentration around nearly black vectors. Electron J Stat. 2014; 8(2):2585–618.

    Google Scholar 

  53. 53

    Meijer CS. Über Whittakersche bzw. Besselsche Funktionen und deren Produkte. Nieuw Arch Wiskd, II Ser. 1936; 18(4):10–39.

    Google Scholar 

  54. 54

    Brychkov YA. Handbook of Special Functions: Derivatives, Integrals, Series and Other Formulas: CRC Press; 2008.

    Google Scholar 

  55. 55

    Brunk E, George KW, Alonso-Gutierrez J, Thompson M, Baidoo E, Wang G, Petzold CJ, McCloskey D, Monk J, Yang L, et al. Characterizing strain variation in engineered e. coli using a multi-omics-based workflow. Cell Syst. 2016; 2(5):335–46.

    CAS  PubMed  PubMed Central  Google Scholar 

Download references


Not Applicable.

Availability and requirements

Project name: iCARH.

Project home page:

Operating system(s): Platform independent.

Programming language: R, Stan.

Other requirements: R 3.5 or higher, Stan 2.18 or higher.

License: GNU GPL.

Any restrictions to use by non-academics: None.


Infrastructure support for this work was provided by the NIHR Imperial Biomedical Research Centre. TJ was supported by a Wellcome Trust ISSF Ph.D. studentship. The funding body has not played any roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information




TJ and TE jointly conducted the research and wrote the paper. All authors read and approved the manuscript.

Corresponding author

Correspondence to Takoua Jendoubi.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

Timothy M. D. Ebbels is also a member of the editorial board (Associate Editor) of BMC Bioinformatics.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Jendoubi, T., Ebbels, T.M. Integrative analysis of time course metabolic data and biomarker discovery. BMC Bioinformatics 21, 11 (2020).

Download citation


  • Metabolomics
  • Integrative analysis
  • Biomarker discovery
  • Pathways
  • Bayesian inference