Modeling complex measurement error in microbiome experiments to estimate relative abundances and detection effects

Next generation sequencing (NGS) has profoundly impacted the study of microbial communities (microbiomes), which reside in the host-associated environments (such as the human body) as well as natural environments (such as soils, salt- and freshwater systems, and aquifers). NGS identifies contiguous genetic sequences, which can be grouped and counted to provide a measure of their abundances. However, even when sequences can be contextualized using databases of microbial DNA, NGS measurements of microbial abundances are not proportional to the abundances of organisms in their originating community. In this paper, we consider the problem of recovering microbial relative abundances from noisy and distorted NGS abundance measurements.

Two main avenues of statistical research have considered the analysis of microbial abundance data from NGS: batch effects removal and differential abundance. Batch effects are systemic distortions in observed abundance data due either to true biological variation (e.g., cage/tank effects) or measurement error (e.g., variation in extraction kits). Tools to remove batch effects from microbiome data include both methods adapted from RNA-seq and microarray analysis (e.g., (Leek and Storey, 2007; Johnson et al., 2007; Gagnon-Bartsch and Speed, 2012; Sims et al., 2008)) as well as microbiome-specific approaches (Gibbons et al., 2018; Dai et al., 2019). A key motivation for batch effect removal is to allow for more precise comparisons of biological differences between microbial communities, and methods for “differential abundance” offer one avenue for assessing difference. Differential abundance methods typically aim to detect microbial units (such as species or gene abundances) that are present in the data at different average levels across groups. Differential abundance methods differ widely in generative models for the data (Love et al., 2014; Martin et al., 2019; Li et al., 2021) and/or their approaches to transforming the data before performing regression analyses (Fernandes et al., 2014; Mandal et al., 2015; Mallick et al., 2021).

In this paper, we consider a distinct goal from batch effects removal and differential abundance: modeling the relationship between NGS measurements and the true microbial composition of the community that was sequenced. Our methodology can estimate relative abundances of individual microbial species in specimens by correcting for the unequal detectability of microbial species, accounting for unequal depths of sequencing across samples, and removing batch-specific contamination in samples. In addition, it can estimate species detectabilities, sample intensities, and cross-sample contamination. Our model incorporates information about shared origins of samples (e.g., replicates and dilution series), shared processing of samples (e.g., sequencing or DNA extraction batches) and known information about sample composition (e.g., mock communities and reference protocols), allowing us to identify and estimate statistical parameters that are not accounted for in most models for microbial abundances (e.g., species detectabilities and contamination). While methods exist to estimate either species detectabilities (McLaren et al., 2019; Silverman et al., 2021; Zhao and Satten, 2021) or contamination (Knights et al., 2011; Shenhav et al., 2019; Davis et al., 2018), to our knowledge, no methods can simultaneously estimate species detectabilities, true microbial relative abundances, sample intensities, and contamination.

Because we consider estimation of relative abundances, and we do not assume that all species are present in all samples, our method entails estimation of a parameter that may fall on the boundary of its parameter space. We use estimating equations to construct our estimator, prove semiparametric consistency and existence of limiting laws, construct a highly stable algorithm to find a solution that may fall on the boundary, and develop uncertainty quantification that accounts for the boundary problem. Therefore, in addition to its scientific contribution, our methodology employs contemporary statistical, probabilistic, and algorithmic tools.

We begin by constructing our mean model (Section 2), and discussing identifiability (Section 3), estimation (Section 4), and asymptotic guarantees (Section 5). We demonstrate applications of our method to relative abundance estimation, contamination removal, and the comparison of multiple experimental protocols (Section 6), and evaluate its performance under simulation (Section 7). We conclude with a discussion of our approach and areas for future research (Section 8). Software implementing the methodology are implemented in a R package available at https://github.com/statdivlab/tinyvamp. Code to reproduce all simulations and data analyses are available at https://github.com/statdivlab/tinyvamp_supplementary.

Here we propose a model to connect observed microbial abundances to true relative abundances. We use the term “read” to refer to the observed measurements (which is not required to be integer-valued), “taxon” to the categorical unit under study (e.g., species, strains or cell types), “sample” for the unit of sequencing, and “specimen” for the unique source of genetic material (which may be repeatedly sampled).

Let $W_{i}=(W_{i1},\dots,W_{iJ})$ denote observed reads from taxa $1,\dots,J$ in sample $i$. A common modeling assumption for high-throughput sequencing data is

where $\mathbf{p}_{i}\in\mathbb{S}^{J-1}$ is the unknown relative abundances of taxa $1,\dots,J$ and $\text{exp}(\gamma_{i})\in\mathbb{R}^{+}$ is a sampling intensity parameter (throughout, we let $\mathbb{S}^{J-1}$ denote the closed $J-1$-dimensional simplex). Unfortunately, microbial taxa are not detected equally well (Figure 1 (left); McLaren et al. (2019) and references therein). To account for this, we begin by considering the model

where $\beta=\beta_{1},\dots,\beta_{J}$ represents the detection effects for taxa $1,\dots,J$ in a given experiment ($\circ$ indicates element-wise multiplication). As an identifiability constraint, we set $\beta_{J}=0$, and so interpret $\exp(\beta_{j}),j=1,\ldots,J$, as the degree of multiplicative over- or under-detection of taxon $j$ relative to taxon $J$. We discuss identifiability in more detail in Section 3.

We now generalize model (2) to multiple samples and one or more experimental protocols. For a study involving $n$ samples of $K$ unique specimens ($n\geq K$), we define the sample design matrix $\mathbf{Z}\in\mathbb{R}^{n\times K}$ to link samples to specimens ($\mathbf{Z}_{i\cdot}\in\mathbb{S}^{K-1}$). In most experiments, $Z_{ik}=\mathbf{1}_{\{\text{sample $i$ taken from specimen $k$}\}}$, but more complex designs are also possible (e.g., $Z_{ik}=Z_{ik^{\prime}}=\frac{1}{2}$ denotes that sample $i$ is a 1:1 mixture of specimens $k$ and $k^{\prime}$). Letting $\mathbf{p}$ be a $K\times J$ matrix such that the $k$-th row of $\mathbf{p}$ gives the true relative abundances of taxa in specimen $k$, we have that the relative abundance vector for sample $i$ is $(\mathbf{Z}\mathbf{p})_{i\cdot}$. We allow differing detections across samples to be specified via the detection design matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$. For example, if samples are processed using one of $p$ different protocols, we might specify $X_{iq}=\mathbf{1}_{\{\text{sample $i$ processed with protocol $q$}\}}$. Accordingly, we now consider a detection effect matrix $\boldsymbol{\beta}\in\mathbb{R}^{p\times J}$, with $\boldsymbol{\beta}_{\cdot J}=\mathbf{0}_{p}$ as before. Therefore, one generalization of model (2) is

where $\mathbf{D}_{\boldsymbol{\gamma}}=\text{diag}(\text{exp}(\boldsymbol{\gamma}))$, and where exponentiation is element-wise.

We now extend this model to reflect contributions of contaminant sources (Figure 1 (right)). We consider $\tilde{K}$ sources of contamination with relative abundance profiles given in the rows of $\tilde{\mathbf{p}}\in\mathbb{R}^{\tilde{K}\times J}$. To link sources of contamination to samples, we let $\tilde{\mathbf{Z}}\in\mathbb{R}^{n\times\tilde{K}}$ be a spurious read design matrix. Most commonly we expect $\tilde{Z}_{i\tilde{k}}=\mathbf{1}_{\{\text{source $\tilde{k}$ may contribute reads to sample $i$}\}}$, but we give an example of an analysis with more complex $\tilde{\mathbf{Z}}$ in Section 6.2. Then, along with contaminant read intensities $\tilde{\boldsymbol{\gamma}}=[\tilde{\gamma}_{1},\dots,\tilde{\gamma}_{\tilde{k}}]^{\text{T}}$, we propose to model

Here, $\mathbf{D}_{\boldsymbol{\gamma},\boldsymbol{\alpha}}$ is a diagonal matrix with diagonal elements equal to $\exp(\boldsymbol{\gamma}+\mathbf{V}\boldsymbol{\alpha})$ for design matrix $\mathbf{V}\in\mathbbm{R}^{n\times K_{\alpha}}$ and parameter $\boldsymbol{\alpha}\in\mathbbm{R}^{K_{\alpha}}$. As discussed in more detail in Section 6.2, $\mathbf{V}\boldsymbol{\alpha}$ allows the intensity of contamination to vary across sources, but typically $\boldsymbol{\alpha}=\mathbf{0}_{K_{\alpha}}$. Additionally, while we could incorporate a detection design matrix for contaminant reads (replacing (4) with $\mathbf{D}_{\boldsymbol{\gamma},\boldsymbol{\alpha}}\tilde{\mathbf{Z}}\left[\tilde{\mathbf{p}}\circ\text{exp}\left(\tilde{\boldsymbol{\gamma}}\mathbf{1}_{J}^{T}+\tilde{\mathbf{X}}\boldsymbol{\beta}\right)\right]$ for $\tilde{\mathbf{X}}\in\mathbb{R}^{\tilde{K}\times p}$), for most practical applications it is sufficient to identify $\tilde{\mathbf{p}}$ up to detection effects. Therefore, combining models (3) and (4), we propose the following mean model for next-generation sequencing data $\mathbf{W}\in\mathbb{R}^{n\times J}$:

Our highly flexible mean model (5) encompasses a wide variety of experimental designs and targets of estimation. As a result, without additional knowledge, the parameters in the mean model (5) will generally be unidentifiable. We provide a general condition for identifiability in SI Lemma 1, and prove identifiability the following specific cases:

1. (a)

   Estimating relative detectability under different sampling protocols using specimens of unknown composition (see SI Section 9.2 for identifiability; Section 6.1 for analysis)

2. (b)

   Estimating contamination using dilution series (see SI Section 9.3 for identifiability; Section 6.2 for analysis)

3. (c)

   Estimating contamination, detectability and composition when some samples have known composition (see SI Section 9.4 for identifiability; Section 7 for analysis)

In addition to knowledge of $\mathbf{X}$, $\mathbf{Z}$ and $\tilde{\mathbf{Z}}$, constraints on the entries of $\boldsymbol{\beta}$ and/or $\mathbf{p}$ can result in the identifiability of model parameters. For example, in (c), some rows of $\mathbf{p}$ are known. As mentioned previously, an identifiability constraint on $\beta$ (e.g., $\boldsymbol{\beta}_{\cdot J}=\mathbf{0}_{p}$) is always necessary.

We propose to estimate parameters $\boldsymbol{\theta}^{\star}:=(\boldsymbol{\theta},\boldsymbol{\gamma}):=$ $(\boldsymbol{\beta},\mathbf{p},\tilde{\mathbf{p}},\tilde{\boldsymbol{\gamma}},\boldsymbol{\gamma},\boldsymbol{\alpha})$ as M-estimators. We use likelihoods to define estimating equations but do not require or assume that the distribution of $\mathbf{W}$ lies in any particular parametric class. We show consistency and weak convergence of our estimators of $\boldsymbol{\theta}:=(\boldsymbol{\beta},\mathbf{p},\tilde{\mathbf{p}},\tilde{\boldsymbol{\gamma}},\boldsymbol{\alpha})$ under mild conditions in SI Section 10. We use $\boldsymbol{\theta}_{0}$ to denote the true value of $\boldsymbol{\theta}$. Our unweighted objective is given by a Poisson log-likelihood:

We use $M_{n}(\boldsymbol{\theta})$ to indicate the profile log-likelihood $\text{sup}_{\boldsymbol{\gamma}\in\mathbb{R}^{n}}M^{\star}_{n}(\boldsymbol{\theta},\boldsymbol{\gamma})$, considering the elements of $\boldsymbol{\gamma}\in\mathbb{R}^{n}$ as sample-specific nuisance parameters. Consistency of $\hat{\boldsymbol{\theta}}$ for $\boldsymbol{\theta}_{0}$ does not require $W_{ij}$ to follow a Poisson distribution, however, the proposed estimator will be inefficient if the relationship between $\mathbbm{E}[W_{ij}|Z_{i},X_{i},\tilde{Z}_{i},\gamma_{i},\boldsymbol{\theta}_{0}]$ and $\text{Var}[W_{ij}|Z_{i},X_{i},\tilde{Z}_{i},\gamma_{i},\boldsymbol{\theta}_{0}]$ is not linear (McCullagh, 1983). Therefore, to obtain a more efficient estimator, we also consider maximizing a reweighted Poisson log-likelihood, with weights chosen on the basis of a flexibly estimated mean-variance relationship. We motivate our specific choice of weights via the Poisson score equations (see SI Section 9). For weighting vector $\hat{\mathbf{v}}\in\mathbb{R}_{+}^{nJ}$ such that $\mathbf{1}^{T}\hat{\mathbf{v}}=nJ$, we define the weighted Poisson log-likelihood as

We define $M^{\hat{\mathbf{v}}_{n}}_{n}(\boldsymbol{\theta})$ by analogy with $M_{n}(\boldsymbol{\theta})$ above. We select $\hat{\mathbf{v}}$ via a centered isotonic regression (Oron and Flournoy, 2017) of squared residuals $\text{vec}\big{[}(\mathbf{W}-\boldsymbol{\hat{\mu}})^{2}\big{]}$ on fitted means $\text{vec}\big{[}\boldsymbol{\hat{\mu}}\big{]}$ obtained from the unweighted objective. Full details are given in SI Section 10, but briefly, we set $\hat{v}_{ij}\propto\frac{\hat{\mu}_{ij}+1}{\hat{\sigma}^{2}_{ij}+1}$, where $\hat{\sigma}^{2}_{ij}$ is the monotone regression fitted value for $\mu_{ij}$.

The estimators defined by optima of the weighted or unweighted Poisson likelihoods given above are consistent for the true value of $\boldsymbol{\theta}$ and converge weakly to well-defined limiting distributions at $\sqrt{n}$ rate, though the form of this distribution in general depends on the true value of $\boldsymbol{\theta}$. We leverage an approach from Van der Vaart (2000) to prove consistency, and we combine a bracketing argument with a directional delta method theorem of Dümbgen (1993) to show weak convergence. Proofs are given in SI Section 11.

Computing maximum (weighted) likelihood estimates of $\boldsymbol{\theta}^{\star}$ is a constrained optimization problem, as the relative abundance parameters in our model are simplex-valued, and the estimate may lie on the boundary of the simplex. Therefore, we minimize $f_{n}(\boldsymbol{\theta}^{\star})=-M^{\star}_{n}(\boldsymbol{\theta}^{\star})$ or $f_{n}(\boldsymbol{\theta}^{\star})=-M^{\star\hat{\mathbf{v}}_{n}}_{n}(\boldsymbol{\theta}^{\star})$ in two steps. In the first step, we employ the barrier method, converting our constrained optimization problem into a sequence of unconstrained optimizations, permitting solutions progressively closer to the boundary. That is, for barrier penalty parameter $t$, we update $\boldsymbol{\theta}$ as

and set $t^{(r+1)}=at^{(r)}$ where $a>1$ is a prespecified incrementing factor, iterating until $t^{(r)}>t_{\text{cutoff}}$ for large $t_{\text{cutoff}}$. In practice we find that $t^{(0)}=1$, $a=10$, and $t_{\text{cutoff}}=10^{12}$ yield good performance. We enforce the sum-to-one constraints (9) by reparametrizing $\mathbf{p}$ and $\tilde{\mathbf{p}}$ as $\boldsymbol{\rho}$ and $\tilde{\boldsymbol{\rho}}$, with $\rho_{kj}:=\text{log}\frac{p_{kj}}{p_{kJ}}$ and $\tilde{\rho}_{\tilde{k}j}:=\text{log}\frac{\tilde{p}_{\tilde{k}j}}{\tilde{p}_{\tilde{k}J}}$ for $j=1,\dots,J-1$, which are well-defined because of the logarithmic penalty terms in (8).

In the second step of our optimization procedure, we apply a constrained Newton algorithm within an augmented Lagrangian algorithm to allow elements of $\hat{\mathbf{p}}$ and $\hat{\tilde{\mathbf{p}}}$ to equal zero. Iteratively in each row $\mathbf{p}_{k}$ of $\mathbf{p}$, we approximately solve

where we write $f_{n}(\mathbf{p}_{k})$ as a function of only $\mathbf{p}_{k}$ to reflect that we fix all other parameters at values obtained in previous optimization steps. We choose update directions for $\mathbf{p}_{k}$ via an augmented Lagrangian algorithm of Bazaraa (2006) applied to

where $Q_{k}^{(t)}$ is a quadratic approximation to $f_{n}(\mathbf{p}_{k})$ at $\mathbf{p}_{k}^{(t)}$ and $\nu$ and $\mu$ are chosen using the algorithm of Bazaraa (2006). The augmented Lagrangian algorithm iteratively updates $\nu$ and $\mu$ until solutions to $\mathcal{L}_{k}$ satisfy $|\sum_{j=1}^{J}p_{kj}-1|<\epsilon$ for a small prespecified value of $\epsilon$ (we use $10^{-10}$ by default). Within each iteration of the augmented Lagrangian algorithm, we minimize $\mathcal{L}_{k}$ via fast non-negative least squares to preserve nonnegativity of $\mathbf{p}_{k}$. Through the augmented Lagrangian algorithm, we obtain a value $\mathbf{p}_{\mathcal{L}_{k}}^{(t)}$ of $\mathbf{p}_{k}$ that minimizes $\mathcal{L}_{k}^{(t)}$ (at final values of $\nu$ and $\mu$) subject to nonnegativity constraints. Our update direction for $\mathbf{p}_{k}$ is then given by $\mathbf{s}_{k}^{(t)}=\mathbf{p}^{(t)}_{\mathcal{L}_{k}}-\mathbf{p}^{(t)}_{k}$. We conduct a backtracking line search in direction $\mathbf{s}_{k}^{(t)}$ to find an update $\mathbf{p}_{k}^{(t+1)}$ that decreases $f_{n}(\mathbf{p}_{k})$.

We now address construction of confidence intervals and hypothesis tests. We focus on parameters $\boldsymbol{\beta}$ and $\mathbf{p}$, which we believe to be the most common targets for inference. To derive both marginal confidence intervals and more complex hypothesis tests, we consider a general setting in which we observe some estimate $\hat{\phi}=\phi(\mathbb{P}_{n})$ of population quantity $\phi=\phi(P)$, where $\mathbb{P}_{n}$ is the empirical distribution corresponding to a sample $\left\{\left(\mathbf{W}_{i},\mathbf{Z}_{i},\mathbf{X}_{i},\tilde{\mathbf{Z}}_{i}\right)\right\}_{i=1}^{n}$, $P$ is its population analogue, and $\phi$ is a Hadamard directionally differentiable map into the parameter space $\Theta$ or into $\mathbb{R}$. To derive marginal confidence intervals for $\mathbf{p}$ and $\boldsymbol{\beta}$, we let $\phi(\mathbb{P}_{n})=\boldsymbol{\hat{\theta}}_{n}$ with $\phi(P)=\boldsymbol{\theta}_{0}$. For hypothesis tests involving multiple parameters, we specify $\phi(\mathbb{P}_{n})$ as $2\left[\text{sup}_{\boldsymbol{\theta}\in\Theta}M_{n}(\boldsymbol{\theta})-\text{sup}_{\boldsymbol{\theta}\in\Theta_{0}}M_{n}(\boldsymbol{\theta})\right]$ with population analogue $\phi(P)=0$ under $H_{0}:\boldsymbol{\theta}\in\Theta_{0}$. In each case, we estimate the asymptotic distribution of $a(n)\left(\phi(\mathbb{P}_{n})-\phi(P)\right)$ for an appropriately chosen $a(n)\rightarrow\infty$.

As our model includes parameters that may lie on the boundary of the parameter space, the limiting distributions of our estimators and test statistics in general do not have a simple distributional form (Geyer, 1994), and the multinomial bootstrap will fail to produce asymptotically valid inference (Andrews, 2000). To address this, we employ a Bayesian subsampled bootstrap (Ishwaran et al., 2009), which consistently estimates the asymptotic distribution of our estimators when the true parameter is on the boundary. Let $\mathbbm{P}^{\boldsymbol{\xi}}_{n}$ be a weighted empirical distribution $\sum_{i=1}^{n}\xi_{i,n}\mathbf{1}_{\mathbf{W}_{i}}$ with weights $\boldsymbol{\xi}\sim G$ for $G\sim\text{Dirichlet}\left(\frac{m}{n}\mathbf{1}_{n}\right)$. Then the bootstrap estimator $a(m)\left(\phi\left(\mathbbm{P}_{n}^{\boldsymbol{\xi}}\right)-\phi\left(\mathbbm{P}_{n}\right)\right)$ converges weakly to the limiting distribution of $a(n)\left(\phi\left(\mathbbm{P}_{n}\right)-\phi\left(P\right)\right)$ if we choose $m=m(n)$ such that $\underset{n\rightarrow\infty}{\text{lim}}m=\infty$ and $\underset{n\rightarrow\infty}{\text{lim}}\frac{m}{n}=0$ (Ishwaran et al., 2009). We explore finite-sample behavior of the proposed bootstrap estimators with $m=\sqrt{n}$ in Section 7, finding good Type 1 error control.

To derive marginal confidence intervals for elements of $\boldsymbol{\theta}$, we let $a(n)=\sqrt{n}$ and $\phi(P)=\boldsymbol{\theta}(P)$. Then $\sqrt{m}(\boldsymbol{\hat{\theta}}^{\xi}_{n}-\boldsymbol{\hat{\theta}}_{n})$ has the same limiting distribution as $\sqrt{n}(\boldsymbol{\hat{\theta}}-\boldsymbol{\theta})$. Therefore, for $\hat{L}^{q}_{c}$ the $c$-th bootstrap quantile of the $q$-th element of $\sqrt{m}(\boldsymbol{\hat{\theta}}^{\xi}_{n}-\boldsymbol{\hat{\theta}}_{n})$ and $\hat{\theta}_{q}$ the $q$-th element of the maximum (weighted) likelihood estimate $\boldsymbol{\hat{\theta}}$, an asymptotically $100(1-\alpha)\%$ marginal confidence interval for $\theta_{q}$ is given by $\left(\hat{\theta}_{q}-\frac{1}{\sqrt{n}}\hat{L}^{q}_{1-\alpha/2},~{}\hat{\theta}_{q}-\frac{1}{\sqrt{n}}\hat{L}^{q}_{\alpha/2}\right)$.

As it may be of interest to test hypotheses about multiple parameters while leaving other parameters unrestricted (e.g., $\boldsymbol{\beta}=\mathbf{0}$ with unrestricted elements of $\mathbf{p}$), we also develop a procedure to test hypotheses of the form $H_{0}:\{\theta_{k}=c_{k}:k\in\mathcal{K}_{0}\}$ for a set of parameter indices $\mathcal{K}_{0}$ against alternatives with $\boldsymbol{\theta}$ unrestricted. Letting $\Theta_{0}$ indicate the parameter space under $H_{0}$ and $\Theta$ indicate the full parameter space, we conduct tests using test statistic $T_{n}:=nT(\mathbbm{P}_{n}):=2n[\text{sup}_{\boldsymbol{\theta}\in\Theta}M_{n}(\boldsymbol{\theta})-\text{sup}_{\boldsymbol{\theta}\in\Theta_{0}}M_{n}(\boldsymbol{\theta})]$. As noted above, $T_{n}$ is in general not asymptotically $\chi^{2}$ if (unknown) elements of $\mathbf{p}$ or $\tilde{\mathbf{p}}$ lie at the boundary, and so we instead approximate the null distribution of $T_{n}$ by bootstrap resampling from an empirical distribution projected onto an approximate null; this is closely related to the approach suggested by Hinkley (1988). Let $\mathring{\mu}_{ij}$ denote $\text{exp}(-\gamma_{i})$ times the expectation of $W_{ij}$ under the full model; $\mathring{\mu}_{ij}^{0}$ denote the analogous quantity under $H_{0}$; and define $W_{ij}^{0}=W_{ij}\frac{\mathring{\mu}_{ij}^{0}}{\mathring{\mu}_{ij}}$ $\text{ if }\mathring{\mu}_{ij}>0$ and otherwise set $W_{ij}^{0}=W_{ij}=0$. In practice, we do not know $\mathring{\mu}_{ij}$ or $\mathring{\mu}_{ij}^{0}$, so we replace $W_{ij}^{0}$ with $\hat{W}_{ij}^{0}=\frac{\hat{\mathring{\mu}}_{ij}^{0}}{\hat{\mathring{\mu}}_{ij}}$, where $\hat{\mathring{\mu}}_{ij}$ and $\hat{\mathring{\mu}}_{ij}^{0}$ are, up to proportionality constant $\text{exp}(\hat{\gamma}_{i})$, fitted means for $W_{ij}$ under the full and null models. After constructing $\hat{\mathbf{W}}^{0}$, we rescale its rows so row sums of $\hat{\mathbf{W}}^{0}$ and $\mathbf{W}$ are equal. We then approximate the null distribution of $T$ via bootstrap draws from $mT(\mathbbm{P}^{\boldsymbol{\xi}}_{0n})$ where $\mathbbm{P}^{\boldsymbol{\xi}}_{0n}$ is the Bayesian subsampled bootstrap distribution on $\hat{\mathbf{W}}^{0}_{n}$. We reject $H_{0}$ at level $\alpha$ if the observed likelihood ratio test statistic is larger than the $1-\alpha$ quantile of the bootstrap estimate of its null distribution, or equivalently, if $T_{n}\geq\tilde{L}^{m}_{1-\alpha}$ for $\tilde{L}^{m}_{1-\alpha}$ the $1-\alpha$ quantile of $mT(\mathbbm{P}^{\boldsymbol{\xi}}_{0n})$.

In this paper, we introduce a statistical method to model measurement error due to contamination and differential detection of taxa in microbiome experiments. Our method builds on previous work in several ways. By directly modeling the output of microbiome experiments, we do not rely on data transformations that discard information regarding measurement precision, such as ratio- or proportion-based transformations. This affords our method the key advantage of estimating relative abundances lying on the boundary of the simplex, which is typically precluded by transformation-based approaches. Accordingly, we implement inference tools appropriate to the non-standard parameter space that we consider. The advantage of estimating relative abundances on the boundary of the simplex is not purely theoretical, and we show that our interval estimates do indeed include boundary values, and demonstrate above-nominal empirical coverage in an analysis of data from Karstens et al. (2019). Furthermore, our reweighting estimator allows for flexible mean-variance relationships without the need to specify a parametric model. Our approach to parameter estimation does not assume that observations are counts, and therefore our method can be applied to a wide array of microbiome data types, including proportions, coverages and cell concentrations as well as counts. Finally, our method can accommodate complex experimental designs, including analysis of mixtures of samples, technical replicates, dilution series, detection effects that vary by experimental protocol or specimen type, and contamination impacting multiple samples. Contamination is commonly addressed via “pre-processing” data, thereby conditioning on the decontamination step. In contrast, by simultaneously estimating contamination along with all other model parameters, our approach captures holistic uncertainty in estimation.

Another advantage of our methodology is that we do not require the true composition of any specimens to be known. For example, our test of equal detection effects across protocols in Costea et al. (2017) can be performed without knowledge of specimen composition. Accordingly, our approach provides a framework for comparing experimental protocols in the absence of synthetic communities, which can be challenging to construct.

In addition, we expect that our method may have substantial utility applied to dilution series experiments, as illustrated in our analysis of data from Karstens et al. (2019). Dilution series are relatively low-cost and scalable (especially in comparison to synthetic communities) and may be especially advantageous when the impact of sample contamination on relative abundance estimates is of particular concern. With this said, we strongly recommend against testing the composite null $H_{0}:p_{kj}>0$ against the point alternative $H_{A}:p_{kj}=0$, as these hypotheses are statistically indistinguishable on the basis of any finite sample of reads. However, we are able to determine the degree to which our observations are consistent with the absence of any given taxon, and therefore we can meaningfully test $H_{0}:p_{kj}=0$ against a general alternative.

The focus of our paper was on $\mathbf{p}$ and $\boldsymbol{\beta}$ as targets of inference. Future research could investigate extensions to our model that connect relative abundances $\mathbf{p}$ to covariates of interest, allowing the comparison of average relative abundances across groups defined by covariates, for example. Our proposed bootstrap procedures may also aid the propagation of uncertainty to group-level comparisons of relative abundances in downstream analyses. In addition, while we focus applications of our model on microbiome data, our model could be applied to a broad variety of data structures obtained from high-throughput sequencing, such as single-cell RNAseq. We leave these applications to future work.

As discussed in the main text, our model is flexible enough to encompass a wide variety of experimental designs and targets of estimation, but as a result, the parameters in model (5) are unidentifiable without additional constraints. In this section, we demonstrate identifiability for the three use cases of our model that we considered as illustrative examples. For reference, we provide here the definition of identifiability: