BEEM-Static: Accurate inference of ecological interactions from cross-sectional microbiome data.

CONCLUSION: BEEM-Static provides new opportunities for mining ecologically interpretable interactions and systems insights from the growing corpus of microbiome data.

1. Introduction

Microbial communities represent complex systems that impact various aspects related to human health, e.g. agriculture [1], food processing [2], disease biology [3, 4] and healthcare [5]. Interactions between members of a microbial community determine emergent phenomena such as homeostasis in the ecosystem [6, 7] and overall function of the microbiome [8]. Ecological interactions can be grouped into six major categories including mutualism (positive-positive), competition (negative-negative), antagonism (positive-negative, further includes predation and parasitism), commensalism (positive-neutral), amensalism (negative-neutral) and neutralism (neutral-neutral) [9]. Correspondingly, ecological modeling of microbiomes taking into account such interactions is a key step towards understanding community function [10, 11], forecasting dynamics [12, 13] and rationally designing interventions that alter community structure and function [14].

Advances in high-throughput sequencing and metagenomics have enabled several data-driven approaches to infer microbial interactions, bypassing limitations of experimental approaches in terms of time, resources and cultivability [9, 15]. In particular, correlation-based methods are widely used for their convenient applicability to cross-sectional datasets [9, 16], despite their inability to capture directionality of ecological interactions such as predation and parasitism [9]. Recent studies have also highlighted other pitfalls in correlational analysis, particularly the accuracy of interactions identified even when the data reflects known modes of microbial interactions [17, 18].

Predictive and dynamic modeling of microbiomes based on first-order differential equations (e.g. with generalized Lotka-Volterra models or gLVMs) has found increasing usage and provided useful insights into microbial interactions and dynamics [10, 19, 20]. Wider adoption of such techniques has been hampered by the need for large datasets (as the number of parameters grows quadratically with the number of species) and dense longitudinal sampling to adequately capture fine-grained dynamics [21]. Theoretical assumptions such as the availability of data where all species are at equilibrium (i.e. abundances of species will not change without external perturbation), and where absolute abundances are accurately known, make the determination of gLVM parameters from cross-sectional data solvable in principle [22]. In practice, microbiome data generated from high-throughput sequencing (16S rRNA gene amplicon or whole metagenomic sequencing) provides relative abundances and scaling these accurately enough for gLVM parameter estimation can be challenging [23]. Furthermore, real-world datasets often contain a mixture of perturbed and unperturbed microbiomes where the equilibrium status is unknown, and where data may even come from multiple models [24].

We have previously introduced an expectation maximization (EM) algorithm which couples gLVM parameter and scaling factor estimation (BEEM), specifically for longitudinal microbiome data [23]. Here we transform this EM algorithm to work with cross-sectional data from communities that are at or near equilibrium (BEEM-Static; ). In benchmarking comparisons with simulated communities against 10 other methods that infer microbial interactions from cross-sectional data, we noted that while all other methods only improved slightly over random predictions (AUC-ROC<0.63), BEEM-Static exhibited high accuracy similar to estimation using true scaling values (AUC-ROC>0.88). Similar observations were made with synthetic communities based on all-pair co-culture experiments, where BEEM-Static accurately recapitulated nearly all known interactions and their directionality. Based on statistical filters to identify non-model and/or non-equilibrium samples in real and simulated datasets, we show that BEEM-Static can be robust to up to 40% of data violating these assumptions. Applying BEEM-Static to a large public collection of human gut microbiome profiles (n = 4,617) identified multiple stable equilibria that appear to better reflect ecological enterotypes with distinct carrying capacities and interactions for key species (e.g. Prevotella copri) compared to prior clustering based definitions [24]. BEEM-Static thus provides new opportunities for mining ecologically interpretable interactions and systems insights from the growing corpus of microbiome data in the public domain.

Schematic overview for the BEEM-Static algorithm.

BEEM-Static takes relative abundances from a cross-sectional microbiome dataset (A), and runs an expectation-maximization algorithm (B) to estimate both biomass values (C), the interaction network (D) and carrying capacities (E).

2. Methods

2.1 Estimating biomass and gLVM parameters with cross-sectional data

The gLVM model is a set of differential equations describing the instantaneous growth rate of each species (dx(t)/dt) as a function of absolute cell densities (x(t)) of the p species in a community: where μ is the intrinsic growth rate of species i and β are interaction terms that define the strength of the influence of species j’s abundance on species i’s growth. In general, estimating gLVM parameters (μ and β) requires longitudinal data to measure dx(t)/dt. However, at the non-trivial equilibrium (dx(t)/dt = 0 and x>0): where the time parameter t now becomes implicit in the equation. Dividing both sides by -β and the biomass m, rearranging terms, and re-parameterizing the equation we get: where a = −μ/β, b = −β/β and is the relative abundance of species i at equilibrium (a is also known as the carrying capacity of the species). This equation allows us to estimate gLVM parameters (through a and b) from cross-sectional data, assuming that samples are at equilibrium, and absolute abundances are known (). To account for the fact that microbiome data provides relative abundances and biomass is typically not measured, BEEM-Static extends an EM framework to jointly estimate model parameters and biomass (; [23]):

Estimating biomass (E-step): for a sample, the equation for each species i provides an estimate for the biomass, and BEEM-Static takes the median of estimates across species as a robust estimator for the biomass of the sample:

Estimating model parameters (M-step): BEEM-Static estimates model parameters ( and ) for each species i with sparse regression (implemented with the ‘glmnet’ package in R) in iteration T:

Initialization and termination: BEEM-Static initializes biomass values based on normalization factors from cumulative sum scaling (CSS [25]), with a user defined scaling constant as the median of biomass values (kept constant through EM iterations). The EM process is then run until the maximum number of iterations specified (200 by default) or until convergence when the median of relative changes in biomass values is <10−3. Confidence values (Z-scores) for the final interaction matrix (non-zero off-diagonal entries) were calculated for each species i using forward stepwise regression (implemented in R package “selectiveInference”; version 1.2.5; Akaike information criterion as the stopping criterion).

2.2 Statistical filters to detect violations of modeling assumptions

BEEM-Static uses the following filters to identify and remove samples that violate modeling assumptions and could thus impact model inference:

Equilibrium filter: to identify samples that may not be at equilibrium, BEEM-Static first predicts the abundances of species at equilibrium (based on the current model) for all species that are present: where is the matrix form of (= estimate for ), is a column vector taking the value of (= estimate for ) and is a vector of predicted relative abundances at equilibrium for each species ( is set to 0 if species i is not present). Samples with median relative deviation above a user defined threshold (ϵ1) from these equilibrium values were then excluded as being potentially not at equilibrium ( and ϵ1 = 20% by default).

Model filter: to account for cases where some samples may come from an alternate gLVM, BEEM-Static calculates the median of squared errors for each sample k with respect to the current model parameter estimates ( and ):

This is done in the M-step for each iteration and samples with large median squared error, i.e. (e−median(e))/IQR(e)>ϵ2 where IQR is the inter-quartile range and ϵ2 is a user defined parameter (default value of 3), are then removed for the next iteration’s M-step.

2.3 Selecting shrinkage parameters for sparse regression

The shrinkage parameter λ in the sparse regression penalizes the number of parameters to avoid overfitting and is selected based on five-fold cross-validation in each iteration (selecting the value one standard error away from the best λ [26]). In the M-step of iteration T, a crude selection of is made in BEEM-Static from a large range from 10−10 to 10−1, and then refined with a fine-grained sequence from to .

2.4 Generating simulated datasets

Simulated gLVM data was generated based on previously described procedures [16, 19]. Specifically, to parameterize distributions for generating model parameters, MDSINE [19] was used to estimate the mean and standard deviations of growth rates and inter-/intra-species interaction parameters from the C. difficile infection experiment data provided with the software. Growth rates and intra-species interactions were sampled from normal distributions (forced to be positive and negative respectively to model logistic growth). The interaction network structure was generated by randomly adding edges from one species to another (with probability ranging from 0.1 to 0.5) and the magnitude of the interactions was sampled from a normal distribution (with 0 mean and standard deviation estimated from real data as noted above). Initial abundances of p (30) species were sampled from a uniform distribution (from 0.001 to the mean carrying capacity μ/β of all the p species), with each species having a probability of π to be absent from a sample (π was estimated as the average rate of absence for the top p most prevalent species in all healthy gut microbiome profiles from the database curatedMetagenomicData [27]). A dataset with n samples was generated by numerically integrating the gLVM with the same parameters until equilibrium, starting with n different initial abundance profiles. The abundances of a random time point along the numerical integration (>20% away from the abundance at equilibrium for >50% of the species) was selected as a sample not at equilibrium. Poisson noise was added to the abundance of each species to simulate experimental variability.

2.5 Generating datasets based on growth curves and co-culture experiments

Microbiome profiles for co-culture experiments were taken from a previous study [20] and the relative abundances were scaled using the corresponding biomass measurements (OD600). Six species pairs were randomly selected, and one of the three conditions were randomly picked for each pair in each sample: (1) only the first species was present, (2) only the second species was present and (3) both species were present. For the first two cases, a random timepoint (last 6 timepoints near the equilibrium) was taken from the growth curve (measured by OD600) of species present. For the last case, the scaled abundances of the two species near the steady state (randomly taken from the three replicates) were used. Abundances were re-scaled to relative abundances and the process repeated to generate a dataset with 500 samples. The interaction matrix reported in Venturelli et al [20] was treated as the ground truth (“M-PW1-PW2”).

2.6 Evaluation metrics

We computed the median of relative errors to assess the accuracy of predicted parameters as: where and θ are the estimated and true parameters (, and ) respectively. The area under the receiver operating characteristic curve (AUC-ROC) was computed for the interaction matrix. Z-scores were used to rank interactions (off-diagonal entries only) predicted by BEEM-Static. Sensitivity for predicting the signs of interaction was calculated as the fraction of interactions with correctly predicted signs in the true interaction matrix (non-zero off-diagonal entries only). Sign precision was computed as the fraction of interactions with correctly predicted signs.

2.7 Benchmarking with correlation-based methods

The following correlation-based methods () for inferring interactions from microbiome data were tested.

The -log(p-value) or edge stability were used to rank CCREPE and SPIEC-EASI correlations, respectively, while the absolute values of correlation coefficients were used for the other methods. AUC-ROC was calculated from the lower triangle of the inferred correlation matrix (an inference was considered correct if there was an interaction between the corresponding species regardless of the interaction direction). Sensitivity and precision of signs were calculated as described above, excluding positive-negative interactions as they cannot be differentiated from positive-positive and negative-negative interactions using correlation analysis.

2.8 Analysis of gut microbiome data

Healthy gut microbiome profiles from the database curatedMetagenomicData were preprocessed and used as the standard dataset for learning gLVMs by removing (1) replicate samples, (2) timepoints other than the first timepoint in longitudinal studies, (3) samples under antibiotic treatment and (4) samples from infants. In addition, we included two validation datasets to evaluate different aspects of the model learned by BEEM-Static: (1) all samples from Raymond et al [36] to validate BEEM-Static’s estimated growth and (2) samples from healthy infants (only the first timepoint for each subject) with ages below 12 months to evaluate BEEM-Static’s biomass estimation. Both datasets were also used to validate BEEM-Static’s ability to filter out samples violating model assumptions. To make the number of parameters tractable with the number of data points available, we only kept core species that were present (relative abundance >0.1%) in more than 30% of samples and subsequently removed samples where none of the core species were found, resulting in 42 core species () and 4,617 samples overall. BEEM-Static was applied with the “model filter” (ϵ2 = 0.9) to learn two models (1,995 and 1,145 samples for each model) in two iterations, in which samples violating the filter were removed (1,477 samples removed in total). BEEM-Static was then rerun without the filter on samples assigned to each model separately to re-learn parameters.

2.9 Estimating in situ growth using BEEM-Static and GRiD

With BEEM-Static, in situ growth can be estimated as the deviation from equilibrium: where and are estimated parameters. In addition, species replication rates for samples not under antibiotic treatment in Raymond et al [36] were estimated with the high-throughput mode of GRiD (v1.2.0; default parameters) [37]. Gut microbiome associated genomes provided with GRiD were used as references and read reassignment using pathoscope2 [38] was enabled (parameter “-p”) to resolve ambiguous mappings.

3. Results

3.1 Accurate inference of ecological interactions from cross-sectional microbiome data

To evaluate if ecological interactions can be inferred from cross-sectional microbiome data (16S rRNA gene amplicon or whole metagenomic sequencing) we first began by conducting extensive benchmarks on complex simulated communities where the truth is known (using gLVMs, number of species = 30, equilibrium conditions; Methods). Consistent with prior reports [17, 18], we noted that all 10 state-of-the-art correlation-based methods provided performances that were only slightly better than random predictions in identifying true interactions, even while ignoring directionality (AUC-ROC<0.63; Figs and ). The ability to correct for compositionality of microbiome data (CCREPE, SparCC, CCLasso, REBACCA), or reduce false positives from transitive correlations (MInt, SPIEC-EASI, BAnOCC and gCoda), did not notably change performance compared to naïve correlation calculations (Pearson or Spearman) when inferring ecological interactions. In contrast, BEEM-Static was able to infer the interactions with notably higher sensitivity and specificity than all correlation-based methods (AUC-ROC = 0.88; Figs and ). Although BEEM-Static only used noisy relative abundances for gLVM inference, its performance matched that of a positive control that assumes noise-free biomass values to scale relative abundances to absolute abundances (Figs and , BEEM-Static vs. true biomass), an assumption that is unlikely to be met in realistic settings [23].

Benchmarking performance for network structure and edge directionality.

Note that CCREPE has two versions with Pearson and Spearman correlations (CCREPE.P and CCREPE.S), while SPIEC-EASI using “mb” and “glasso” algorithms is represented as SE.mb and SE.glasso, respectively. (A) ROC curves for different correlation based methods, the regression method using true biomass values and BEEM-Static for one simulated dataset with 30 species and 500 samples. (B) Boxplots showing precision and recall for directionality/sign of interactions for 30 different simulated communities with 30 species and 500 samples each. (C) True interaction network for a synthetic community based on all-pair co-culture data (Ground truth), inferred correlation network by BAnOCC and inferred interaction network by BEEM-Static (numbers on edges represent gLVM parameter values or correlation coefficients). PC: Prevotella copri, BV: Bacteroides vulgatus, BU: Bacteroides uniformis, BO: Bacteroides ovatus, BT: Bacteroides thetaiotaomicron, FP: Faecalibacterium prausnitzii, BH: Blautia hydrogenotrophica, ER: Eubacterium rectale, CA: Collinsella aerofaciens, EL: Eggerthella lenta, DP: Desulfovibrio piger, CH: Clostridium hiranonis.

In addition to knowing that two species are interacting, a key aspect of ecological interactions is the directionality/sign of interactions, with correlation-based methods either exhibiting low sensitivity (REBACCA, SPIEC-EASI) or low precision (Pearson, Spearman, CCREPE, SparCC, CCLasso, MInt, BAnOCC, gCoda) despite the exclusion of predatory interactions (positive-negative) when calculating their performance (; Methods). BEEM-Static addresses this issue by providing high sensitivity (>80%) and precision (nearly 100%) using cross-sectional microbiome data (). These observations were recapitulated in a wide range of simulated datasets with varying network structure and edge sparsity, highlighting BEEM-Static’s robustness ( and ). Furthermore, BEEM-Static provides estimates for biomass values that were found to be consistently accurate (relative error <10%) and can be used to provide meaningful biological insights [23].

We extended the evaluations to experimental data, using all-pair co-culture and isolate growth curves for 12 species [20] to create synthetic communities where the interaction network is known, albeit sparse (; Ground truth). Not surprisingly, recapitulating the structure of such a simple interaction network was not difficult for most correlation-based methods, with BAnOCC having the best performance overall (AUC-ROC = 0.9, ). However, determining directionality of interactions was still a challenge, despite the simplicity of the network. For example, in the case of BAnOCC the commensalistic interaction between DP and FP was captured as a negative correlation, while the predatory interaction between BH and ER was captured as a positive correlation (; BAnOCC). BEEM-Static, on the other hand, was able to capture all interactions and their directionality correctly, except for one false positive (FP to DP) and one false negative (BT to EL) involving interactions with weak strength (; BEEM-Static). BEEM-Static’s utility in such datasets was consistently observed in comparison to correlation-based methods with AUC-ROC close to 1 ().

3.2 Statistical filters in BEEM-Static provide robustness to violations in modeling assumptions

While the simulations in the previous section account for experimental errors, they assume that all samples come from equilibrium states, an unlikely situation for most real datasets. Relaxing this assumption, we noted that with as little as 5% non-equilibrium samples, AUC-ROC decreased by >10%, and with 15% non-equilibrium samples AUC-ROC performance degraded to match that of correlation-based methods (Naïve algorithm, Figs and ; Methods). Incorporating a statistical filter in BEEM-Static that compares estimated species relative abundances at equilibrium with observed abundances (equilibrium filter; Methods) helped identify samples that were not at equilibrium with high sensitivity and specificity (). This in turn allowed BEEM-Static to be robust to having nearly half of the samples (45%) in the dataset being at a non-equilibrium state (performance reduction <5%; Figs and , BEEM-Static).

BEEM-Static robustly filters samples violating modeling assumption in simulated and real datasets.

(A-B) Performance reduction for BEEM-Static (with filters) and the naïve algorithm (without filtering of samples) as the percentage of samples at (A) equilibrium or (B) generated from the main model, decreases. Reduction is measured relative to BEEM-Static with no filters and all data from the model and at equilibrium. Points denote the means while error bars denote the standard deviation across 30 simulations each. (C-D) Principal coordinates plots (Bray-Curtis dissimilarity) representing gut microbiome taxonomic profiles from 4,617 samples. Points represent samples from individuals taking antibiotics (C) or from newborn infants (D), while crosses represent samples from adults who are not undergoing antibiotic treatment. Points that were filtered by BEEM-Static are colored blue and red otherwise.

We next investigated the impact of relaxing the “universal model” assumption i.e. that all samples have the same ecological conditions and model parameters. This assumption may not hold true in many real-world settings (e.g. gut microbiome samples from different enterotypes [24]), and as expected relaxing it had a strong impact (30% reduction in AUC-ROC) even with slight deviations (5%, Naïve algorithm; Figs and ; Methods). In addition, AUC-ROC performance continued to decrease beyond 60% even after nearly half the samples (40–45%) were derived from a different model. To address this, BEEM-Static implements a filter that identifies samples that have poor goodness-of-fit to the current model (model filter; Methods) and excludes them in subsequent iterations of model inference. This approach was found to provide robustness to up to ~40% of samples from a different model (performance reduction <20%, BEEM-Static; Figs and ).

Finally, we employed real microbiome datasets to test BEEM-Static’s robustness where a subset of samples is known to violate modeling assumptions i.e. some subjects undergoing oral antibiotic treatment or samples from newborn infants, where the majority of samples are from adults who are not undergoing antibiotic treatment. BEEM-Static was able to identify such samples with high sensitivity (>80%) using its model filter ( and ), and the filtered samples were significantly enriched for those from antibiotic-treated adults and infants (Fisher’s Exact test p-value<10−22).

3.3 BEEM-Static analysis of human gut microbiomes identifies distinct ecological configurations

To further assess BEEM-Static’s utility we evaluated the concordance of parameters learnt during the training process with orthogonal information for a large human gut microbiome dataset (N = 4,617; Methods). In particular, we noted that biomass estimates from BEEM-Static were significantly higher for adults versus newborn infants (~2×; Wilcoxon test p-value <10−15; ), consistent with our understanding of a maturing gut microbiome [39, 40]. Additionally, we used deviations from equilibrium (dx(t)/dt = 0) to estimate instantaneous growth (population increase or decrease) of each species in each sample (Methods), and assessed concordance with an in silico approach to estimate DNA replication rates [37]. Despite the fact that growth rates are also impacted by death rates, we observed that species predicted to grow based on BEEM-Static analysis were also found to have significantly higher DNA replication rates (GRiD values; Wilcoxon test p-value = 3×10−4; ).

Analysis of human gut microbiomes with BEEM-Static.

(A) Violin plots showing the significant difference in BEEM-Static biomass estimates for adults and newborn infants. (B) Boxplots showing differences in DNA replication rates (GRiD estimates) for species predicted to decrease and increase in population size, respectively, by BEEM-Static. Each point represents one species in a sample. (C) Scatter plot showing the first two components from a principal component analysis of equilibrium abundances for samples (as predicted by BEEM-Static). (D) Predicted carrying capacities (square root transformed) of species in the two models. Species with divergent carrying capacities (ratio is >2 standard deviations from 1) in the two models are marked with stars.

Using an iterative approach to train a new model on samples excluded from the first model, we observed that BEEM-Static classifies a majority of the samples to two models (1995 for model 1 and 1145 for model 2; ). Visualizing the expected microbial composition for each sample at equilibrium, we noted three distinct clusters in a principal component analysis plot (), where samples in clusters 1 and 2 largely correspond to profiles from models 1 and 2, respectively. Analysis of carrying capacities highlighted that the clusters were defined by Prevotella copri, which has higher carrying capacity in cluster 1 vs 2 and is absent in samples from cluster 3 (). Several other species were also found to have divergent carrying capacities in the two models, including Bifidobacterium adolescentis, Collinsella aerofaciens and Coprococcus comes (enriched in model 2), some of which have been shown to be associated with fiber intake [41, 42]. Overall the ecological models identified here have distinct sets of interactions () and do not appear to match earlier definitions of enterotypes [24] based on principal coordinate analysis and clustering (S).

4. Discussion

As microbiome research increasingly moves from descriptive studies to those that seek to provide a mechanistic understanding of microbial communities, the ability to infer microbial interactions from microbiome data is an important capability. In particular, the directionality and sign of interactions provide biologically interpretable information that is missed by correlation-based approaches. BEEM-Static provides an alternative avenue to infer this, with the caveats that it assumes a specific model (generalized Lotka-Volterra) for community dynamics and removes rare and low abundance species with limited number of samples in practice. In addition, as we show here, other values obtained from BEEM-Static models can have utility, including the strength of interactions, biomass estimates, deviation from equilibrium, and fit to model.

In addition to accounting for relative abundance estimates from microbiome data, the statistical filters employed by BEEM-Static make it robust to some of the violations in model assumptions that can be expected in real datasets. These features make BEEM-Static widely applicable, and also extends the use of ecological models with microbiome data. For instance, our analysis of large public microbiome datasets provides an alternate perspective to the discussion on microbial enterotypes [24] and universality of microbiome dynamics [43]. The ecological types observed here are characterized by distinct carrying capacities that might be a function of the environment (e.g. host factors or diet). Fiber rich diets are known to have a strong impact on the gut microbiome [42] and have been linked to some of the species with differential carrying capacities in our models [41]. We anticipate that the incorporation of such environmental factors into future models would be an exciting avenue to study their influence on microbial community structure in vivo. Finally, hybrid methods that learn models from both longitudinal and cross-sectional data represent another promising direction to explore for studying general and individual specific microbiome dynamics [23].

Performance comparison for determining interaction network structure.

Boxplots show AUC-ROC values for 30 different simulated communities with 30 species and 500 samples each.

(TIFF)

Click here for additional data file.

Robustness of BEEM-Static performance across different network structures.

Boxplots showing relative error for (A) parameter estimates for biomass (scaled to have the same median as the ground truth) and carrying capacity, and (B) interaction network AUC-ROC as well as sign recall and precision, based on 30 simulations. Diamonds mark median performance using true biomass values. Different network structures (except “random”) were generated with the SPIEC-EASI R package.

(TIFF)

Click here for additional data file.

Robustness of BEEM-Static performance with varying edge densities.

Boxplots showing relative error for (A) parameter estimates for biomass (scaled to have the same median as the ground truth) and carrying capacity, and (B) interaction network AUC-ROC as well as sign recall and precision, based on 30 simulations. Diamonds mark median performance using true biomass values.

(TIFF)

Click here for additional data file.

Benchmarking on synthetic communities created based on co-growth experimental data.

Boxplots show the performance of various methods based on 30 replicates that use randomly selected interacting pairs.

(TIFF)

Click here for additional data file.

Utility of BEEM-Static filters for avoiding performance reduction in the presence of model violations.

Results are shown for increasing proportion of samples that are, (A) not at equilibrium or (B) not from the main model. The naïve algorithm is without filtering of samples and performance reduction is measured relative to BEEM-Static with no model violations for the data. Points show means while error bars show standard deviation across 30 simulations.

(TIFF)

Click here for additional data file.

Percentage of the gut microbiome represented by the 42 core species used by BEEM-Static for modeling with the curatedMetagenomicData dataset.

(TIFF)

Click here for additional data file.

Ecological models from BEEM-static represent distinct configurations from enterotypes.

Principal coordinates plot (Bray-Curtis dissimilarity) based on gut microbiome taxonomic profiles. The contour lines highlight the two different regions where points aggregate. Points are colored by the BEEM-static models that they belong to (unassigned samples not shown), showing that while there appears to be some enrichment, model ids and enterotypes do not show a 1–1 correspondence.

(TIFF)

Click here for additional data file.

Parameters estimated by BEEM-Static for the two models learned from the curatedMetagenomicData dataset.

(XLSX)

Click here for additional data file.

Derivation of the expectation-maximization algorithm and training convergence.

(PDF)

Click here for additional data file.

10 May 2021

Dear Dr. Nagarajan,

Thank you very much for submitting your manuscript "BEEM-Static: Accurate inference of ecological interactions from cross-sectional metagenomic data" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Nicola Segata

Associate Editor

PLOS Computational Biology

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Li et al. developed a tool called "BEEM-Static" to infer microbial interactions from cross-sectional metagenomic data. These interactions are inferred using a generalized Lotka-Volterra model (gLVM). The authors highlight that in contrast to correlation-based methods BEEM-Static is able to infer directionality in addition to presence of the interactions. Further, they provide estimates of the sample biomass. They benchmarked their algorithm against 10 correlation-based methods using simulated datasets, simple pairwise experimental co-culture data and a large cross-sectional metagenomic dataset. BEEM-Static was robust to violations of model assumptions (equilibrium states and model differences) in these tests. While their method is interesting and the benchmarking results show improvements over other tools, I have several concerns:

Major concerns

• What is the main difference to the original BEEM algorithm? The filter functions? Please add more details on that to the introduction.

• Further details about model training should be provided. For example, additional figures showing how well the model converged by the end of training.

• The parameter inference method used in the paper is a coordinate descent algorithm, a general form of EM-algorithm. For the EM algorithm, it is necessary to write the log-likelihood function and formulate the E-step and M-step in this framework. That way, the M-step is maximizing the lower bound of the log-likelihood function which is an essential guarantee for the model to converge. In the general coordinate descent, it is necessary to prove that the function is bounded in order to show that the algorithm will converge. This is missing currently.

• In the large-scale metagenomic dataset, only 42 species passed the filtering. This number is very small for human gut microbial communities and a large fraction of the microbial community is not considered in many samples. I would suggest the authors provide a barplot showing the fraction of considered microbes in each sample.

• The authors didn't report how many samples are found invalid in the large-scale metagenomic dataset. They showed that the model still works well when up to 40% of samples are at non-equilibrium. Without the number of how many samples are predicted as "invalid" in the metagenomic dataset, it is uncertain if this criteria is met for the large metagenomic dataset.

• I would recommend a graphical overview of BEEM-static as one of the first figures, including input, output and potential conclusions that can be inferred from the results.

• I’m curious what the rationale is behind the biomass estimation. Why should it be possible to infer biomass based on relative abundances?

• Last paragraph of results section: Why do you point out fibre might play a role here? Without the context from the cohort, this seems a bit random. Please describe the real world datasets you used briefly (Healty cohort? Disease?) and why fibre is relevant.

• Why did you expect that enterotypes play a role in this or are affected by interactions?

Minor concerns:

• I don’t see why their approach is only applicable to metagenomic data. The abundance profiles from 16S data should work equally well.

• It would be interesting to verify the model's prediction on the equilibrium state of the metagenomic samples from the other time points. Does the next timepoint deviate from the current one? Does the second time point of the samples predicted as invalid deviate more than those predicted as valid? That would provide further confidence for the model's performance.

• Given a model, there will be multiple equilibrium states (the eigenvectors of the state matrix). It would be interesting to show how these equilibrium states are distributed in the large cohort.

• Some terms should be defined in the introduction, such as equilibrium. Examples for bacterial interactions of the different types would be helpful too (pos-pos, neg-neg, pos-neg).

• p. 15 (on my version): Define the abbreviation for the species involved in the interaction (DP, FP, ...).

• Fig 1C: Explain in the legend what those numbers in connection with the interactions mean.

• Fig 3C: The axis relative to the origin would make more sense to me. Also, I’m curious why the first PC explains so much variance? How were the clusters computed?

Reviewer #2: This manuscript introduces a new bioinformatics algorithm,

"BEEM-Static", to infer ecological interactions from cross-sectional

metagenomic data. The work is motivated by the current needs in

obtaining deeper ecological understanding of the systems underlying

observation in microbiome research in general, and in microbial

metagenomics in particular. Currently popular methods for

computational inference of microbial interactions are based on

observed co-abundance patterns but often lacking ecological

insight. This paper aims to demonstrate that ecological structure can

be recovered from cross-sectional metagenomic data. The method is an

EM algorithm (BEEM-Static) that utilizes a standard ecological (gLV)

interaction model. Application requires sample filtering but this is

justified based on the modeling assumptions: in particular, the

authors show that gLV parameters can estimated directly from

cross-sectional data, when it is being assumed that 1) the system is

close to equilibrium, and 2) absolute abundances are known. This has

considerable practical value because it helps to obtain ecologically

relevant insights from dynamical systems in a setting where the lack

of long and dense time series has so far limited the

investigations. Comparison with relevant alternatives demonstrates

remarkable improvements in the ability to detect interaction and its

directionality reliably (based on AUC/ROC scores). It is also shown

that the method is relatively robust to outliers and deviations from

stationarity or model assumptions, so that the first assumption seems

feasible. The second assumption is not generally available but the

authors propose a way to estimate this computationally, and advances

in measurement techniques may improve this further over time; the main

contribution in this work is the method itself. An interesting

demonstration of the model performance is provided by analysing

community stability in a well-known public human gut microbiome data

set, and the proposed model is potentially applicable to studying key

aspects of community dynamics across a vast spectrum of microbiome

data sets, and the results can be experimentally testable in many

cases if suitable experimental time series can be collected.

The manuscript is technically sound and written in a fluent and easily

understandable English. Experiments and statistical analyses have been

conducted rigorously and the conclusions are supported by the

data. Relevant literature is cited, and the relevant recent methods

have been included in the comparisons.

All data and code are fully available but I have not tried to

replicated the analysis. The associated R package seems good quality

and implemented according to prevailing standards for R package

development.

I only have the following minor notes:

* Minor

- The ecological models do not correspond to previously reported

enterotypes (in particular Fig S6). However, enterotypes - as

visualized on 2D PCoA space - are very broad community-level

clusters. It is possible that there are more refined subcommunities

that would be better visible on further PCoA dimensions or in a

focused subcommunity analysis and these might correspond better with

the previously reported enterotypes that are said to be driven by

specific bacterial groups. It would be informative to include some

more discussion on the possible connections. In particular,

Prevotella has been shown to have bimodal/bistable abundance

distribution at the population level, and if the BEEM-Static model

does not match this finding, this would be somewhat unexpected and

potentially indicate problems in the BEEM-Static model.

- The method is promoted as a new method for metagenomic data; but it

seems to me that this applies more generally to ecological community

profiling data, including 16S and potentially not even limited to

microbial ecological communities. Some discussion on the wider scope

might be warranted.

- The model is applied for core species that are abundant and

prevalent; applicability to lower abundance and rare groups is a

possible limitation and could be more clearly stated in Discussion.

- It would be informative and helpful to show how the term in M-step

inside the median (a / \\Sigma bx) can be derived for estimating

absolute abundances.

- The "EM framework" is not exactly the classical EM algorithm

although very close conceptually; reconsider naming e.g. "iterative

EM-type framework" or otherwise acknowledging the difference.

- The bullet point list in 2.7 might be better written out as text or

incorporated as a table?

- principle component analysis (common misspelling) -> should be:

principal component analysis; same for principal coordinate analysis

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Leo Lahti

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

19 Jul 2021

Submitted filename: response letter.pdf

Click here for additional data file.

11 Aug 2021

Dear Dr. Nagarajan,

We are pleased to inform you that your manuscript 'BEEM-Static: Accurate inference of ecological interactions from cross-sectional microbiome data' has been provisionally accepted for publication in PLOS Computational Biology. There are a couple of minor points still highlighted by one of the reviewer, but I trust you will address them in the final version of the manuscript.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Nicola Segata

Associate Editor

PLOS Computational Biology

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have addressed my concerns and I have no further comments.

Reviewer #2: The manuscript has improved, and the responses are adequate.

The remaining comments:

- In review responses the authors cite the Remien et al. preprint (https://www.biorxiv.org/content/10.1101/463372v1.full) discussing the estimation of absolute abundances; in fact Remien et al. show the expected result that absolute abundances cannot be fully identified from compositional data. This is in contrast with the approach in BEEM-Static, which provides an estimate of absolute abundances based on compositional data. The authors have mentioned, however, that this is only an approximation. It is important to make sure that the inability to fully estimate absolute abundances from compositional data are mentioned clearly when discussing the limitations of the method. Better estimates may become available as absolute quantification methods are developing, further increasing the potential value of this method. therefore I do not consider the inaccessibility of absolute estimates a major limitation in this case.

- Text is understandable but I suggest to pay attention to final proof-reading and figure quality

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Leo Lahti

1 Sep 2021

PCOMPBIOL-D-21-00585R1

BEEM-Static: Accurate inference of ecological interactions from cross-sectional microbiome data

Dear Dr Nagarajan,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Andrea Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol

Table 1

Correlation methods tested.

Algorithm	Version	Parameters	Note
Pearson and Spearman correlation			Directly calculated from relative abundances
CCREPE [28]	v1.2.0	1000 iterations	Correction done for both Pearson and Spearman correlations
SparCC [29]	commit id: 9a1142c	Default
CCLasso [30]	v1.0	Default
REBACCA [31]	https://tinyurl.com/ymdts4wy	nbootstrap = 50, B = 500, FWER = 0.01
MInt [32]	v1.0.1	Default
SPIEC-EASI [33]	v0.1.4	lambda.min.ratio = 0.01, nlambda = 20, rep.num = 50	Both “mb” and “glasso” algorithms
BAnOCC [34]	v1.0.1	Default
gCoda [35]	commit id: 584bd07	Default