A statistical model for microbial community profiles

SparseDOSSA is a hierarchical model for microbial count and relative abundance profiles (Fig 1), with components specifically accommodating the major distributional characteristics of such data, namely zero-inflation, compositionality (and thus sequencing depth), feature-feature non-independence, feature-environment interactions, and high-dimensionality. Briefly (Fig 1A, details in Methods), the model a) specifies zero-inflated log-normal marginal distributions for each microbial feature to allow for both biological and technical absences, b) imposes distributions on the “absolute”, i.e. pre-normalized, microbial abundances to satisfy compositionality (similar to models such as the Dirichlet [20] or gamma-MGH [18]), c) models feature-feature correlations through a multivariate Gaussian copula [21], and d) adopts a penalized fitting procedure to address high-dimensionality [22]. Conditional on feature relative abundances and total read depth, count observations are modeled with a standard multinomial sampling procedure, and per-sample read depth is modelled with a log-normal distribution. For implementation, we adopted a penalized Expectation-Maximization procedure for model fitting, and we have evaluated and provided options for cross-validated selection of the optimal penalization parameter (Methods). SparseDOSSA 2 is implemented as an R/Bioconductor package (http://huttenhower.sph.harvard.edu/sparsedossa2). Note that this differs somewhat from our previous implementation of SparseDOSSA that was already publicly available; we provide a summary of changes between the versions in Methods. SparseDOSSA can be accurately fit to a wide variety of different microbial community structures to capture both (inferred) absolute and relative count observations (Fig 1B and 1C, details below).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. A hierarchical model for microbial community feature profiles.

A) SparseDOSSA comprises a hierarchical model to capture the generation mechanism of microbial sequencing counts, including components for “hidden” absolute abundances, sequencing depth (and thus compositional relative abundances), zero inflation, and feature-feature and feature-environment interactions. Notations not defined in the figure: : cumulative density function (CDF) for the absolute abundance of feature A_j. μ_D, : mean and variance of the log normal sequencing depth distribution. B) SparseDOSSA can be fitted to varied microbial community types using cross-validation procedures by users; the software also provides pre-trained models are provided for human microbiome template datasets. This allows for C) simulation of either null or "true positive" association spiked-in synthetic datasets, to facilitate microbiome benchmarking or power analysis studies.

https://doi.org/10.1371/journal.pcbi.1008913.g001

SparseDOSSA accurately recapitulates real-world microbial community structures

We validated SparseDOSSA’s ability to accurately capture realistic microbial community feature profiles by quantifying its performance across a variety of real-world datasets (Fig 2 and S1 and S2 Tables). The studies used include: 1,2) taxonomic profiles from shotgun sequenced metagenomes of healthy human stool and posterior fornix samples from the HMP1-II, hereafter referred to as “Stool” and “Vaginal” [23], 3) shotgun sequenced stool metagenomes of inflammatory bowel disease (IBD) patients from the HMP2 Inflammatory Bowel Disease Multi-omics Database (IBDMDB, abbreviated as “IBD”) [2], and 4) 16S rRNA gene sequenced murine distal gut communities after diet perturbation [24]. By evaluating the model in different cohorts, we established its robustness under different community phenotypes, habitats (i.e. body sites), overall ecological structures, and sequencing technologies (S1 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. SparseDOSSA accurately recapitulates different microbial community structures.

We compared SparseDOSSA 2 simulated microbial counts versus those of three human microbiome training template datasets (Stool, Vaginal, and IBD). A) Bray-Curtis ordination shows global agreement between SparseDOSSA simulated microbial abundance profiles and those of their originating real-world populations. B) This was quantified by PERMANOVA R² statistics, showing that SparseDOSSA simulated samples were significantly less systematically differentiated from their targets than existing DM and metaSPARSim methods in almost all cases (Wilcoxon rank sum test p-values included in S3 Table). R² compared against randomly split original real-world data are included as baseline controls. C) Representative features from each environment are similarly distributed between real-world and SparseDOSSA simulated samples, as shown in empirical cumulative distribution functions (CDFs) of log-10 relative abundances (with pseudo value 1e-6 to visually represent zeros). D) Per-feature Kolmogorov-Smirnov summary statistics quantify that SparseDOSSA outperforms existing methods in simulating realistic feature-level relative abundance distributions. First, the similarity between the model-simulated feature abundance distribution versus that in the real-world dataset is quantified with K-S statistics. Then, the K-S statistics for SparseDOSSA and the other two models (DM and metaSPARSim) are plotted on the x- and y-axis, respectively (each point representing one feature, smaller K-S statistics represent better approximation). Lastly, the K-S statistics of SparseDOSSA versus other models are formally tested using Wilcoxon signed rank tests (p-values are significant and included in S4 Table).

https://doi.org/10.1371/journal.pcbi.1008913.g002

SparseDOSSA 2 captured community parameters and re-simulated microbial profiles with overall community structures that accurately reflected those of the original, real-world ecologies, better than alternative methods (Fig 2A and 2B), across all human datasets (murine study results reported in separate section). Overall, simulated communities yielded the same patterns of global beta-diversity as were contained within each modeled dataset (Fig 2A). This was quantitatively compared against alternative models (Dirichlet-multinomial, DM [20] and gamma-MGH, namely metaSPARSim [18]) with the PERMANOVA R² statistic [25] (Methods). We calculated ecological Bray-Curtis dissimilarities between real-world microbial profiles and those simulated by each evaluated method. We then quantified the total variability in the combined dissimilarities that could be attributed to real-world versus simulation difference, expressed as the PERMANOVA R². Smaller R²s thus indicate less deviation of the simulated community structures from the real-world target and better performance of the model.

Across almost all evaluated community types, SparseDOSSA 2 generated significantly smaller R² statistics over 25 simulation iterations than existing methods (Wilcoxon rank sum tests p < 0.05), indicating better fit to and recapitulation of the targeted communities (Fig 2B, testing results in S3 Table). Notably, this was consistent in both the human gut (Stool, IBD), where community structure forms continuous “gradients” of microbial composition [9], and the human vaginal environment (Vaginal), where communities are often characterized by a few discrete types with dominant species [26]. Only for the Stool dataset did SparseDOSSA slightly underperform when compared to metaSPARSim in terms of R² statistic, while still outperforming with respect to per-feature distributions (Fig 2D). Additionally, metaSPARSim’s simulation procedure can take as much as ~10x longer than SparseDOSSA 2 (S1 Fig), which is prohibitive for realistic data sizes (especially for benchmarking or power-analysis efforts requiring multiple simulations per parameter configuration, or for Monte-Carlo calculations). Based on the above comparative performance results, we conclude that, when evaluated for overall community structures, SparseDOSSA can capture microbial feature profiles that closely resemble those of real-world microbiomes.

The SparseDOSSA model also provided the best recapitulations of individual features’ relative abundances (Fig 2C and 2D). For representative features in each environment, the empirical cumulative distribution function (CDF) curves of samples show that SparseDOSSA simulated abundances closely resemble those of the real-world data (Fig 2C). Quantitatively, for each set of microbial features, we measured the difference of distributions between re-simulated and real-world (modeled) relative abundances with the Kolmogorov-Smirnov summary statistic (K-S, see Methods). The resulting K-S summaries provides a distance between the distribution of each feature’s relative abundances across simulated vs. modeled real-world communities. Smaller K-S statistics thus indicate better performance of the model. SparseDOSSA better approximated the targeted real-world per-feature distributions than existing methods across all evaluated datasets (Fig 2D), reaching statistical significance in each case (S4 Table, Wilcoxon signed rank tests p < 0.05). In addition to simulating existing microbial features, SparseDOSSA 2 also provides the functionality to simulate new features that resemble the targeted environment’s ecological characteristics (Methods) and was validated to generate “Stool-like”, “Vaginal-like”, or “IBD-like” new features in terms of prevalence, abundance, and variability for each of the tested datasets (S2 Fig). Thus both in overall community structure modeling and in per-feature models, SparseDOSSA was able to accurately capture and re-simulate realistic microbial observations better than alternative approaches.

SparseDOSSA captures covariation among microbes and with real or simulated “phenotypes”

Once the SparseDOSSA model is fit to a real-world microbial community profile, the “reversed” version of the model can be used not only to simulate similar, controlled ecologies, but to introduce artificial, known feature-feature or feature-covariate associations to the simulated profiles (i.e. metadata “spike-ins”). This is implemented by first capturing the “null” state of targeted real-world studies as described above and by subsequently modifying the fit model parameters to induce artificial associations. Compared to existing spiking-in paradigms [5,17,18], the model includes two important improvements (Methods). First, SparseDOSSA can model a wide variety of covariates—discrete, continuous, or any combination thereof—with multivariable linear modelling, and can thus accommodate simulations of realistic microbiome population study designs with multiple phenotypes, exposures, or confounders [1]. Second, associations with both non-zero (abundance) and zero-inflated (prevalence) components of microbial features can be captured, along with clearly defined effect sizes (fold change or odds ratio, see Methods). This enables rigorous evaluations of, for example, differential abundance testing methods for their statistical performance (e.g. power or false positive rates).

Based on models fit to the Stool and Vaginal communities, SparseDOSSA 2 accurately introduced, or “spiked-in” associations for control phenotypes in a new, simulated population (Figs 3 and S3). Specifically, for the Stool dataset, we introduced a binary covariate (similar to e.g. a case / control contrast) with non-zero effects on 16 (5% of the total 332) microbial features’ abundances (Fig 3A) and prevalences (Fig 3B). Features were selected to ensure the highest effective sample size (Methods). For simulated associations of the “phenotype” with feature abundances, log fold change of non-zero relative abundances largely agreed with the target effect sizes within 95% confidence levels (Fig 3A). Prevalence log odds ratios were also as targeted (Fig 3B), with effects in relative abundances mostly agreeing with the prescribed effect sizes. Similar abundance and prevalence results were consistently reproduced in the Vaginal environment (S3 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. SparseDOSSA can add feature-phenotype and feature-feature associations to modeled microbial community simulations.

A,B) SparseDOSSA 2 correctly simulated feature-phenotype associations targeting the prescribed non-zero relative abundance (A) and prevalence (B) effect sizes of the spiked features, while maintaining non-associations of null features. True associated (spiked) microbial features (red) are well differentiated from null features (black), through Bonferroni corrected p-values (non-significant features marked in gray; test based on linear/generalized linear regression against the spiked metadata variable, see Methods for details). The horizontal dashed lines indicate true spike-in effect sizes: red lines for the positive and negative true effect sizes, respectively, and the black line for null effect (0). C) SparseDOSSA can also prescribe feature-feature associations. Bottom right triangles are Spearman correlations in the simulated absolute abundances. As prescribed, only true association feature pairs are correlated. Top right triangles are Spearman correlations in the corresponding, simulated relative abundances. Note that in this example, Spearman correlation does not differentiate between true (“biological”) covariations versus those induced spuriously due to compositionality (as is also the case in the underlying data on which SparseDOSSA’s model is fit). As expected, both true signals and spurious correlations caused by compositionality can be observed for such data. TP: true positives.

https://doi.org/10.1371/journal.pcbi.1008913.g003

In addition to modeling associations between microbial features and external covariates, SparseDOSSA can also model community ecological interactions (i.e. correlations between microbial features, or feature-feature “spike-ins”, Figs 3C and S4). First, we note that SparseDOSSA simulated microbial communities naturally recapitulate important synergies within modeled environments, such as Firmicutes-Bacteroidetes gradients in the human gut and co-exclusion of dominant Lactobacillus species in the vaginal microbiome (S5 Fig). For introducing artificial correlation among microbial features, we extended the feature-covariate spiking process above to synthetically associate multiple features with the same hidden covariate (Methods). First, a null model fit to the Stool/Vaginal communities contains no true feature-feature associations, only those that manifest spuriously due to compositionality. Starting with this, we modified the model to induce increasingly large feature-feature “ecological” interactions. SparseDOSSA 2 produced both only and exactly the expected true feature-feature associations among absolute abundance components, and the correct induced compositional correlations after simulating the sequencing assay process (Fig 3C for subset of Stool results; full Stool/Vaginal results in S4 Fig). These results support SparseDOSSA’s ability to modify baseline, null community structures by the introduction of interactions among features or with controlled covariates, which together enable the evaluation of a wide range of statistical approaches to microbiome analysis [15,27,28].

Modeling environment-specific benchmarking and power estimation

Since most microbiome analysis methods make simplifying assumptions that may or may not be suited to particular ecologies, SparseDOSSA’s flexible model enables power and accuracy estimation in a habitat-specific manner (Fig 4). Specifically, by spiking only a limited set of known feature-phenotype associations into an otherwise guaranteed-null model, differential abundance methods can be compared directly to each other in a controlled setting (more details in [29]), enabling targeted method benchmarking (Fig 4A) or power analysis (Fig 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. SparseDOSSA enables comparative benchmarking and power analysis of microbial community statistical association tests.

For any originating community type of interest, datasets simulated based on a SparseDOSSA model fit can be spiked with known "phenotypes" and feature effect sizes to estimate methods performance (power, FPR, etc.) during (A) benchmarking as well as (B) power analysis, across controlled combinations of potential effect sizes and sample sizes. Points indicate average performance across simulation repetitions and error bars indicate standard error (Methods).

https://doi.org/10.1371/journal.pcbi.1008913.g004

To demonstrate SparseDOSSA’s use for benchmark comparison of microbial community statistical tests, we again simulated synthetic datasets based on the Stool profiles with "phenotypic" associations spiked-in for 5% of features at varying effect sizes (as in Fig 3A). Multiple replicates of the same parameter set were performed to provide performance metric mean and standard errors (Methods). Using the resulting gold standards, the performances of three different association tests–limmaVOOM [30], ANCOM [31], and MaAsLin 2 [29]—were similar for power, but false discovery rates varied strikingly (Fig 4A). Notably, the MaAsLin 2 generalized linear model showed good FDR control at small to moderate effect sizes. At higher effect sizes, non spiked-in (“null”) features are also called by MaAsLin 2 as differentially abundant. Interestingly, this is because SparseDOSSA’s spike-in effects are imposed on features’ simulated absolute abundances (Methods), and high effect size spike-ins thus also induce relative abundance change in null features due to compositionality. This highlights the important difference between true differential abundance effects corresponding to microbes’ biological variation, versus changes post normalization that are driven by other features.

In contrast, ANCOM [31] was designed to account for compositionality and draw inference about hidden absolute abundances; it successfully and maintained FDR under moderate to strong effect sizes. Arguably as a result, however, its performance suffered for small to null effects, presumably because in such cases it is difficult to distinguish between “driver” microbial features with true absolute effects versus those with changes in their relative abundances due to compositionality. Lastly, limmaVoom [30], designed primarily for RNA-Seq data, had inflated FDRs across all cases.

To demonstrate SparseDOSSA’s use for power analysis during microbial community study design, we focused targeted simulation datasets with spiked-in effects on a feature modeled on Escherichia coli, as a microbe commonly associated with dysbiosis in the human gut [2]. Using this approach, SparseDOSSA 2 can be used to estimate each association method’s expected power for similar biomarkers and populations. In this example, MaAsLin 2 has high power to detect a two-fold abundance change in "E. coli" for a sample size of at least ~500 individuals, but greatly reduced power for smaller fold-changes (Fig 4B). Since model power for differential abundance testing in sparse, compositional data is extremely difficult to determine parametrically, SparseDOSSA thus provides a way to do so by simulation tailored to any community type or feature of interest, and we provide basic guidelines for this in Methods.

SparseDOSSA reproduces an end-to-end diet-microbiome analysis

In many cases, SparseDOSSA thus captures the properties of microbial community ecologies well enough to reproduce surprisingly specific aspects of their membership and distributions, which we next demonstrated by reproducing an end-to-end example from a longitudinal, interventional diet study investigating the effects of diet on the murine gut microbiome [24] (Fig 5). Carmody et al. [24] used 16S rRNA gene sequencing to profile changes in the mouse distal gut microbiome under different dietary treatments (chow, raw and cooked tuber, and meat). To determine whether SparseDOSSA could accurately model the microbes, phenotypes, and associations observed over time in these settings, we fitted model parameters for each sample type at different time points and under different treatment assignments (Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. SparseDOSSA correctly models the effects of diet and time on the murine gut microbiome by reproducing effects from amplicon sequencing profiles.

A) SparseDOSSA 2 was fitted to subsets of samples from [24] that included up to three time points each from collections of mice fed chow, raw or cooked tubers, and meat. The resulting models were then used to simulate controlled microbial community profiles, which correctly reproduced the beta-diversity structures present in the original study (MDS ordination by Bray-Curtis dissimilarities, corresponding to Fig 1A of [24]). The SparseDOSSA model was also able to model and synthetically replicate changes in "Bacteroidetes" and "Firmicutes" phyla in response to raw vs. cooked diets, including B) overall community alpha-diversity (Shannon index), C) the resulting "Firmicutes" vs. "Bacteroidetes" ratio, and D) overall whole-community effective biomass. These correspond to [24]’s Fig 1F–1H, respectively. TRF = raw tuber (free-fd); TCF = cooked tuber (free-fed); TCR = cooked tuber (restricted ration).

https://doi.org/10.1371/journal.pcbi.1008913.g005

After re-simulating communities based on these models fits, ordinations of SparseDOSSA results closely mimic the originally observed clustering structure of dietary effects, and even the longitudinal effects of time under treatment (Fig 5A, corresponding to Fig 1A of [24]). For quantitative differential abundance effects, based on the observed difference of raw diet (TRF) samples when compared against cooked/free-fed (TCF) and cooked/restricted (TCR) within the tuber diet group [24], we additionally applied SparseDOSSA 2’s spike-in procedure to simulate a ~2x fold increase in the abundance of Bacteroidetes OTUs in TCF/TCR when compared to TCF samples, and a ~2x fold decrease in Firmicutes OTUs (Methods). Consequently, the simulated samples displayed similar differential outcomes in community diversity as measured by Shannon index, as well as Firmicutes/Bacteroidetes ratios, as seen in the original study (Fig 5B and 5C, corresponding Fig 1F–1G of [24]).

Interestingly, even though per-feature absolute abundances are theoretically unidentifiable in the SparseDOSSA model (Methods), we note the spiking-procedure recapitulated the decreased total cell counts in TRF (Fig 5D, corresponding to Fig 1H of [24]) that [24] also observed via quantitative PCR. The difference between chow versus tuber diets, on the other hand, is completely attributable to SparseDOSSA’s framework, as one can arbitrarily modify the average absolute abundances of our fit to the chow or tuber diet samples, but still yield exactly the same relative abundance profiles. These real-world application results highlight SparseDOSSA’s adaptability to community phenotypes and treatment effects, as well as demonstrate its performance for amplicon sequence datasets and microbial communities associated with non-human hosts.

The SparseDOSSA model

SparseDOSSA uses the following data generation mechanism to parameterize microbial community profiles: a) environments/samples contain microbes with absolute abundances A b) these are normalized to relative abundances X, which c) can be measured via sequenced counts C. As detailed in Fig 1A, our model specification for these components is:

• For the unobserved absolute abundances A = (A₁, A₂, …, A_p), we specify a Gaussian copula model [21] with zero-inflated log normal marginal distributions. Specifically, this involves assuming hidden multivariate Gaussian variables g = (g₁, g₂, …, g_p) for the microbial features and a mapping of these variables to the corresponding absolute abundances (A₁, A₂, …, A_p):
  1. ○. g ~ MVN(0, Ω⁻¹). That is, each g_j is a standard N(0, 1) variable and their correlation matrix is Ω⁻¹.
  2. ○. Each g_j is mapped to A_j such that A_j follows a zero-inflated log-normal distribution, parameterized by absence probability (π_j) and mean and variability of non-zero log abundances (μ_j, σ_j²):
     1. ■. A_j = 0 if g_j < Φ⁻¹(π_j)
     2. ■. if g_j ≥ Φ⁻¹(π_j)
  Where Φ is the standard normal cumulative density function and is the cumulative density function of the zero-inflated log-normal distribution, parameterized by π_j, μ_j, σ_j².

It follows from our model specification that, marginally, A_j follows the prescribed zero-inflated log-normal distribution exactly:

1. ○. With probability π_j, A_j = 0
2. ○. With probability 1 − π_j, log A_j ~ N(μ_j, σ_j²)

Jointly, correlations between microbial features’ absolute abundances are characterized through the copula parameter Ω. The benefit of adopting a copula model is to separate the parameterization and estimation of a joint distribution into its marginal and correlation components; this is illustrated in the model fitting subsection below.

• Relative abundances are directly normalized from absolute abundances: . This by definition satisfies compositionality (). Also note that because A_j’s are zero-inflated, this directly induces zero-inflation (i.e. biological absence) in X_j’s.
• For a given sample i with sequencing depth D_i, its per-feature read counts (C_i1, C_i2, …, C_ip) are assumed to follow a multinomial distribution with individual features’ probabilities given by X_ij. That is, (C_i1, C_i2, …, C_ip) ~ MultiNom(D_i, X_i1, X_i2, …, X_ip), thus also allowing technical zeros.
• Lastly, we assume the sequencing depth D_i across samples follows a log-normal distribution. That is, D_i ~ LogN(μ_D, σ_D²).

Model likelihood

It is helpful to clarify the likelihood of our model given its parameterization. First, we derive f_A, the likelihood for the unobserved absolute abundances A. The likelihood of observed data, as we show later, is an integration of f_A. For illustration purposes, we first note the special case where A_j are not zero-inflated. That is, π_j = 0 for all j’s. In this case, we have that:

Where g_j is as defined above: and ϕ(·) is the standard normal density function. The equality follows by noting that the second term (the product) is the Jacobian of the mapping . When one or more A_j’s are zero-inflated, the mapping g→A is not one-to-one, and the right hand side of the equality requires integration over g_j’s that map to zero-valued A_j’s:

To derive the likelihood for relative abundances X, we note that X, jointly with the total absolute abundance A^Σ (A^Σ≔∑_jA_j), forms a one-to-one mapping with the absolute abundances A (A = A^ΣX). Thus, the density function for X, f_X, can be obtained through integration of , which is simply f_A multiplied by the Jacobian of the transformation: (1)

Lastly, for the observed microbial count data C, the proper likelihood is:

Where f_C|X(C) is the multinomial likelihood for microbial counts given their relative abundances. In practice, to simplify computation, during model fitting we replace this likelihood with (2)

Where is the multinomial MLE for X given observed C, i.e., , const is a normalizing constant not involving the parameters. The approximation is acceptable because with modern sequencing depth [1], f_C|X(C|X) (as function of X) is highly concentrated around . The right-hand side of (2) is what we aim to maximize for estimation of our model’s parameters, π₁, μ₁, σ₁²,⋯,π_p,μ_p,σ_p², and Ω.

Identifiability

It is important to note that likelihood (1) is unidentifiable. That is, there exist different values of the parameter set (π₁, μ₁, σ₁²,⋯,π_p, μ_p, σ_p², Ω) that yield the same likelihood f_X (and consequently f_C). Intuitively, this is because X are normalized from absolute abundances A, and different A values can map to the same normalized relative abundances X—this is thus typical of any compositional setting. Regarding the identifiability of our parameters, we build on the results of [43], which is a special case of our model where π_j = 0 for all j’s. Specifically, we note that:

• π_j’s are identifiable, as X_j = 0 ⇔ A_j = 0
• μ₁,⋯,μ_p are identifiable up to a constant. That is, μ₁,⋯,μ_p and μ₁+c,⋯,μ_p+c lead to the same likelihood, for any constant c. For this reason, in our model estimation we impose the (arbitrary) constraint that ∑_jμ_j = 0.
• σ₁²,⋯,σ_p² are identifiable, given μ₁,⋯,μ_p and Ω. One can note that when π_j = 1 for all j’s, our likelihood degenerates to that in [43] with explicit analytical forms.
• Ω is not identifiable. Again, consider the special case that π_j = 1 and σ_j = 1, the form of f_X is explicit and involves , which is a multiple-to-one mapping from Ω. The issue of non-identifiable correlation matrices for microbiome abundance data has been noted and addressed in many previous works; refer to [15,27,28,43] for a partial list. We adopt the technique used in many of these previous works, namely L₁ penalization on Ω, to simultaneously address the identifiability issue as well as high-dimensionality for generic estimation of large covariance matrices [22].

We note that, importantly, good performance of SparseDOSSA does not require correct estimation of its non-identifiable parameters. Rather, so long as SparseDOSSA provides a configuration of estimated parameters that yields a joint distribution of relative abundances (X₁,⋯,X_p) that closely approximates that of the original data, its application case for capturing or simulating microbial observations is achieved. Thus, our proposed solution for non-identifiability, though not guaranteed to correctly estimate the unknown true parameters, does satisfy typical applications.

Model fitting

Given our model specification and its (non-)identifiability, we propose to minimize the following penalized negative log-likelihood function for solving the parameter set Θ = (π₁, μ₁, σ₁²,⋯,π_p, μ_p, σ_p², Ω) (the sequencing depth parameters (μ_D, σ_D²) can be fitted independently on per-sample read depths with maximum likelihood):

Subject to the constraint for μ_j as specified above: ∑_jμ_j = 0. As such, λ>0 is a penalizing tuning parameter, which we choose with K-fold cross-validation in practice (five folds by default). can be either existing relative abundance estimations or, as specified above, normalized from count observations ().

As specified in (1), the likelihood function f_X involves integration over A^Σ and is not analytically tractable. Numerically, we propose the following penalized expectation-maximization algorithm [44] for model fitting:

1. Initialize by fitting a multivariate log-normal distribution on .
2. During the r-th iteration:
   1. E-step: calculate expectation .
   2. Penalized M-step: maximize −l^(r)(X|Θ)+λ||Ω||₁ with respect to Θ to obtain . Note that do not require updates. and can be solved analytically. can be solved with standard graphical lasso [22].
3. Iterate until convergence.

Generating synthetic microbial observations and simulating new features

Given that our model is fully parametric, synthetic microbial observations, including (hidden) absolute abundances, normalized relative abundances, and sequencing counts, can be generated following the same specifications as described above. To provide model parameters, the user can adopt one of the pre-trained sets included with the software or use the SparseDOSSA 2 training procedure to estimate parameters from any microbial template dataset suited for their simulation case.

Users may also be interested in generating “new” microbial features from the same ecological environment. For this, SparseDOSSA additionally models the per-feature parameters (π_j, μ_j, σ_j²) with a three-dimensional non-parametric distribution F. That is, across features, (π_j, μ_j, σ_j²)~F. Given a set of SparseDOSSA fitted results , F can be estimated with a three-dimensional normal kernel density estimator [45]. The estimated can then be used to simulate new microbial features that follow the ecological characteristics (i.e., prevalence, abundance, and variability) of the fitted environment (S2 Fig).

Association spike-in

SparseDOSSA adopts linear and generalized linear models for flexible spiking-in in both microbial features’ non-zero abundances and prevalences, based on covariates. Let Z_i be the vector of covariate(s) for sample i and β be the targeted corresponding effect sizes (coefficients). To spike in associations between feature j’s abundance and covariates Z_i, we modify the feature’s non-zero mean log absolute abundance parameter μ_j across samples. Specifically, the post spike-in mean log abundance is modified as

For the i-th sample, A_ij can be generated with ZILogN(π_j, μ_ij, σ_j²) instead of the original ZILogN(π_j, μ_j, σ_j²). This dictates that A_ij’s are associated with Z_i in their mean non-zero log abundances. As μ_ij are specified on the log scale, β_abundance, by definition, corresponds to log fold changes.

The prevalence spike-in similarly is specified via the logistic model; we modify the presence probability parameter (1−π_j) across samples:

And generate A_ij’s correspondingly. This introduces an association between the covariates Z_i and feature j’s prevalence, with β_prevalence corresponding to log odds ratios of the feature being present. The multivariate linear modelling approach for specifying the association effects for both abundance and prevalence allows us to flexibly incorporate different variable types (e.g. binary, continuous, etc.) and study designs (e.g. existence of confounders).

We note that, importantly, our spiking-in procedure is performed on the absolute abundances, A, which induces differential effects in relative abundances X (Fig 3A and 3B). The main benefit of this approach is that both the spiked-in microbial features and the “null” (i.e. non-spiked features) are clearly defined. The alternative—specifying effects for X_j—is conceptually difficult. As X is compositional (sums to 1), prescribing enrichment effects (higher abundance or prevalence) for some microbial features must by definition lead to depletion effects for certain other features. This renders it difficult to clearly define and separate the set of "true positive" spiked-in microbial features and the set of null features. SparseDOSSA’s definition of effects for absolute abundances in its spike-in procedure align with recent efforts to rigorously characterize microbial differential abundance effects under the constraint of compositionality [31,46]. Empirically, we note that prescribed log fold changes or odds ratios for A_j often lead to similar effect sizes in the relative abundances X_j for the spiked-in feature j’s (Fig 3A and 3B).

Lastly, we note that the spiking-in procedure with metadata variables can be used to simulate association effects between pairs of microbial features (Fig 3C). Specifically, we first simulate a hidden covariate Z with standard normal distribution. For a pair of features j₁, and j₂, to enforce positive correlations between the two absolute abundances and , we simulate for them to be associated with Z in the same direction:

Where β_abundance = β_prevalence = β can be viewed here as the effect size specifying the strength of correlation between and . To spike in negative correlation between the two, we simply keep β as the effect for one of the features and use −β for the other.

Comparison between SparseDOSSA and SparseDOSSA 2

A previous implementation of SparseDOSSA has also been publicly available [47]. To help the reader understand the difference between SparseDOSSA’s old and new versions, we list the major updates below:

1. In the new version, zero-inflated log-normal modelling of microbial feature abundances is imposed on absolute abundances A_ij, instead of directly on the relative abundances X_ij, to better satisfy compositional constraints of such data types.
2. Microbial feature-feature correlations are now explicitly modelled with the Gaussian copula component; this was not present in the previous version.
3. For simulating new microbial features, which requires modelling of the triplet (π_j, μ_j, σ_j²), the previous implementation adopted a log-normal distribution on μ_j, with corresponding linear and logit link functions for σ_j and π_j. This is relaxed in the new version with the three-dimensional nonparametric Gaussian kernel density estimator to allow for additional flexibilities.
4. For spiking-in of association microbial features and metadata variables, the previous version only allows for abundance association spike-in. It also specifies effect sizes through a non-linear transformation that can be difficult for the user to interpret. The new version expanded this functionality by allowing both abundance and presence spiking-in; it also specifies effect sizes directly through linear and generalized linear models (log fold change/odds ratio), allowing more straightforward interpretations.
5. The necessary algorithm and implementation changes for the aforementioned model updates, including our new EM algorithm for model fitting, updated and expanded real-world templates for direct simulation, and software interface updates.

Evaluation with real-world datasets

For evaluation and comparison of microbiome simulation methods, we examined three real-world datasets with different host environments and disease statuses (S1 Table) [2,23]. We subset publicly available species level profiles from [23] (all healthy) to baseline time point stool (Stool) and posterior fornix (Vaginal) samples, and those from [2] to baseline time point IBD samples (IBD, including Crohn’s disease and ulcerative colitis patients). We additionally removed samples with lower than ~3,000 reads mapped to identified taxa and features present in less than 3 samples. The datasets’ dimensions (sample size, number of features), post processing and filtering, are included in S1 Table. A summary of all simulation evaluations performed in this work is included in S2 Table. Evaluation scripts used are publicly available at https://github.com/biobakery/sparsedossa_paper.

To evaluate the performance of individual methods, we randomly partitioned each dataset (Stool, Vaginal, IBD) into halves for five iterations. For each partitioning, we fit the parameterization/simulation methods (DM, metaSPARSim, and SparseDOSSA 2) on one half of the data (training). We then simulated synthetic microbial observations with the same sample size using the fitted results. Lastly, we compared the synthetic observations with the other half of the partitioned data (testing) in terms of both overall dissimilarity with PERMANOVA [25] and per-feature distribution differences (methods detailed below). Within each partitioning, this simulation was also performed five times for each method. That is, each method was used to randomly simulate five synthetic datasets for comparison with the testing half. The partitioning procedure allows us to evaluate method performance without a model overfitting effect. The DM was fitted using R package “dirmult” and metaSPARSim was fitted using the implementation referred to in its publication [18]. For metaSPARSim fitting, the percent not zero filter for features was set to 0 instead of the default 0.2. In our evaluation this led to an observed performance increase (thus a favorable assessment), likely due to the existence of many highly zero inflated microbial species.

To evaluate overall dissimilarity between the original and synthetic samples, for each partitioning we combined the testing half of the original samples with the simulated datasets (five for each fitted method). We calculated the sample-to-sample Bray-Curtis dissimilarity matrix D on the combined dataset. The univariate PERMANOVA model, D~I{Sample is simulated} was fitted. The corresponding R² statistic quantifies the percentage of variability between samples attributable to the difference between original "real-world" as compared to simulated samples. Smaller R² statistics indicate less difference, and better performance of the simulation method. For each method, a total of 25 PERMANOVA evaluations (5 original dataset partitioning × 5 simulation) was performed for each real-world dataset. Lastly, we additionally evaluated the R² between the training and testing halves of a dataset for each partitioning; this yields an estimation of minimum achievable R²’s for each dataset. R²’s between methods are compared with Wilcoxon rank sum tests (S3 Table).

To evaluate the difference between distributions of individual features in the original and synthetic datasets, we simply combined the synthetic datasets generated across all partitioning and simulation repetitions. An individual feature was compared for its relative abundance distribution between the original real-world data and combined synthetic samples. This was quantified with the Kolmogorov-Smirnov (K-S) summary statistic, which is defined as the largest absolute difference between the empirical cumulative distribution functions of the real-world and synthetic abundances. Smaller K-S statistics indicate better approximation of the targeted real-world distributions with the simulation method. To compare overall per-feature performance between models, the paired K-S summary statistics between SparseDOSSA 2 and alternative models (DM and metaSPARSim) were compared using Wilcoxon signed rank tests (matched by features, S4 Table).

Association spike-in evaluation

We simulated spiked-in associations between microbial features and a synthetic case/control variable, based on the SparseDOSSA 2 fitted results. A total of 1,000 synthetic samples were simulated (500 cases and 500 controls). For non-zero abundance spike-in (Fig 3A), the top 5% (16 total) most prevalent features were selected for spiking-in; this yields the highest effective sample size for the selected features because our abundance spiking-in targets only the non-zero component of a feature’s distribution. Half of the features were spiked for a targeted log fold change (β_abundance) of 1 in cases compared to controls, and the other half were spiked for a log fold change of -1. Actual log fold changes in the simulated relative abundances, along with 95% confidence intervals, were calculated by performing a linear regression on the log transformed non-zero relative abundances for each feature.

Similarly, for prevalence spiking-in (Fig 3B), the top 5% features with prevalence closest to 0.5 were selected; as with abundance spike-ins, this was to ensure the spiked-in features had the highest effective sample size, as the association between a binary outcome (presence/absence here) and a binary covariate (case/control) is best-powered when the sample distribution is balanced across all different outcome/covariate combinations [48]. Again, half of these features were spiked for a targeted log odds ratio (β_prevalence) of 1 in cases compared to controls, and the other half were spiked for a log odds ratio of -1. Actual log odds ratios of the simulated feature prevalence, along with 95% confidence intervals, were calculated by performing a logistic regression on the presence for each feature.

For simulation of feature-feature associations, we first set the correlation between feature pairs in SparseDOSSA 2 to zero (i.e., Ω = I where I is the identity matrix). This ensures that feature absolute abundances are independent in the “null” dataset (Fig 3C left panel, bottom right), whereas spurious correlation still exists in relative abundances due to compositionality. Two random pairs (four features) in the top ten most abundant features were selected for non-zero feature-feature association spike-in. As specified in Methods above, we simulated two independent normal synthetic hidden metadata variables, one for each feature pair to be associated. For the first feature pair, they were spiked in both abundance and prevalence with the same effect (β = β_abundance = β_prevalence) at varying sizes, for positive association. The second pair were spiked with opposing effects (β for one, −β for the other) for negative association. We used Spearman correlation to estimate the empirical association between feature pairs in the simulated absolute and relative abundances. Target association effect size was also varied (β of 0, 1, 2, and 5) to showcase the relative signals of “true” associations that exist for both absolute and relative abundances, and spurious associations that are only induced in relative abundances due to compositionality (Figs 3C and S4).

Benchmarking

Since “true” associations with prescribed effect sizes are known for SparseDOSSA synthetic datasets, they can be used for benchmarking microbiome analysis methods as well as for power analysis of microbiome study designs. For benchmarking analysis (Fig 4A), we again selected the top 5% (16 total) most prevalent features in the Stool dataset to perform abundance spike-in, such that the selected features had the highest effective sample size. A total of 200 microbial profiles were simulated to be associated with a balanced binary metadata (100 cases, 100 controls). We varied effect sizes with half spiked features at β_abundance = (0, 0.5, 1.2) and the other half with β_abundance = (0, −0.5, −1, −2), correspondingly (in the effect size 0 case no spike-in was performed and microbial profiles are generated independently of metadata). A total of 500 random simulations were performed for each parameter combination. We applied existing differential abundance analysis methods to detect the spiked-in features in each simulation dataset [29–31], with individual method configurations as reported in our previous benchmarking analysis [29]. We summarized the empirical power and FDR of a method in one simulation dataset, across the twenty random replicates for each parameter configuration, and reported the mean and standard error in Fig 4A.

Power analysis

For showcasing SparseDOSSA’s utility in a power analysis, we spiked in non-zero abundance associations with a balanced case-control variable for a simulated species parameterized by fitting Escherichia coli. This was performed at varying effect sizes (log fold change, β_abundance = (0.5, 1, 2)) and sample sizes (100 to 1000). For each parameter configuration, a total of 500 replicates were simulated. The empirical power and its standard error of using MaAsLin 2 to detect the differential abundant effect in "E. coli" was summarized across the 500 replicates and reported in Fig 4B; this was repeated for each effect size/sample size configuration.

While it remains to build a more formal power analysis using these models, we recommend the following guidelines if using SparseDOSSA directly:

• Based on the planned study design, if a pilot data set is available, the fitting functionality can be used to estimate the parameters that would generate similar data sets. Otherwise, researchers can adopt one of the pre-fitted templates that best approximate their target environment.
• Researchers can then simulate many datasets (e.g. hundreds or thousands) with the fitted parameters, across a range of spike-in effect sizes and sample sizes. We note that the SparseDOSSA 2 effect size has straightforward quantitative interpretations: log fold changes for abundance spike-in, and log odds-ratios for prevalence spike-in (details on model in Association spike-in subsection of Methods).
• This should be followed by conducting the hypothesis testing that is planned for the study on each of the simulated datasets. This allows for a measure of how often the “true” (spiked-in) associations are detected for a given sample size (i.e. statistical power). Additionally, the researcher can investigate whether other methods of testing can deliver comparable results to simulation-based estimates from SparseDOSSA 2.

Murine diet microbiome analysis

We applied SparseDOSSA 2 to the longitudinal diet dataset of the mouse gut microbiome in [24], to show that our method is capable of reproducing a complex study’s findings. To recapitulate the longitudinal diet effect as reported in [24]’s Fig 1A, we fitted SparseDOSSA 2 separately on 1) the control Chow diet samples at baseline, 2,3) Tuber diet samples at day 1 and day 5, separately, and 4,5) Meat diet samples at day 1 and day 5, separately. This approach allows SparseDOSSA 2 to independently fit subsets of the data, without assuming a priori the observed differences noted in [24]. We then used SparseDOSSA 2 fitted results to simulate synthetic observations for each diet/timepoint combination, with five times the original sample size (to reduce variability due to random sampling). Bray-Curtis MDS ordination on these synthetic data displayed a striking resemblance to that observed in [24] (Fig 5A), in that a) communities cluster according to dietary treatment, and 2) this response is consistent after one day of switching from chow to whole-food diets and is strengthened at day 5.

We next reproduced the differential gut microbial profiles observed in mice fed raw versus tuber diets as presented in [24]’s Fig 1F–1G. [24] adopted three different types of Tuber diet: the raw/free-fed (TRF), the cooked/free-fed (TCF), and the cooked/restricted (TCR). This study presented that on the phylum level, TRF induced enrichment of Bacteroidetes and depletion of Firmicutes when compared to TCF/TCR. We applied SparseDOSSA 2’s feature spike-in procedure to approximate this effect. Specifically, we generated a balanced, three category (TRF/TCF/TCR) variable. Based on our fitted model of the Tuber diet at day 5, we spiked in a two-fold (β_abundance = log 2) increase in the non-zero abundance of Bacteroidetes OTUs and a two-fold decrease in Firmicutes OTUs (β_abundance = −log 2) in TRF samples when compared to TCF/TCR samples. This roughly agrees with the presented results in [24] Fig 1D. We next simulated SparseDOSSA synthetic datasets for both the baseline Chow diet samples (sample size = 20), and the spiked-in Tuber diet samples (60 samples total, 20 each for TRF/TCF/TCR). We calculated the Shannon index and Firmicutes/Bacteroidetes ratios of these samples, and show in Fig 5B that they agreed with the corresponding findings presented in [24] Fig 1F–1G.