A zero-inflated Bayesian nonparametric approach for identifying differentially abundant taxa in multigroup microbiome data with covariates

Scientific studies conducted during the last two decades have established the central role of the microbiome in disease and health. Differential abundance analysis aims to identify microbial taxa associated with two or more sample groups defined by attributes such as disease subtype, geography, or environmental condition. The results, in turn, help clinical practitioners and researchers diagnose disease and develop new treatments more effectively. However, detecting differential abundance is uniquely challenging due to the high dimensionality, collinearity, sparsity, and compositionality of microbiome data. Further, there is a critical need for unified statistical approaches that can directly compare more than two groups and appropriately adjust for covariates. We develop a zero-inflated Bayesian nonparametric (ZIBNP) methodology that meets the multipronged challenges posed by microbiome data and identifies differentially abundant taxa in two or more groups, while also accounting for sample-specific covariates. The proposed hierarchical model flexibly adapts to unique data characteristics, casts the typically high proportion of zeros in a censoring framework, and mitigates high dimensionality and collinearity issues by utilizing the dimension reducing property of the semiparametric Chinese restaurant process. The approach relates the microbiome sampling depths to inferential precision and conforms with the compositional nature of microbiome data. In simulation studies and in the analyses of the CAnine Microbiome during Parasitism (CAMP) data on infected and uninfected dogs, and the Global Gut microbiome data on human subjects belonging to three geographical regions, we compare ZIBNP with established statistical methods for differential abundance analysis in the presence of covariates.

1. Introduction. During the last two decades, numerous scientific studies have investigated the role of microbiome in disease diagnostics and therapeutics to establish the influence of microbiome variations on health and disease (Wade, 2013;Shreiner, Kao and Young, 2015;Rautava, 2016;Tang, Li and Hazen, 2019). Giant strides in next-generation sequencing (NGS) technologies have enabled the simultaneous identification and quantification of relative abundances of all microbial taxa in study samples, offering a glimpse into the microbiome in our bodies (Tang, Li and Hazen, 2019).
Most study samples are sequenced using either amplicon or random shotgun sequencing. The former sequencing technique amplifies the highly variable regions of the 16S rRNA gene (Kembel et al., 2012). Shotgun metagenomic sequencing covers all genetic variation in the DNA compared to the specific regions targeted by amplicon sequencing. Consequently, shotgun sequencing achieves better species-level resolution and ability to detect novel viruses (Bibby, 2013;Relman, 2013), but requires much larger volumes (Sharpton, 2014). We refer the reader to Bharti and Grimm (2019) and Xia, Sun and Chen (2018a) for a comprehensive overview, along with the pros and cons, of different microbiome sequencing technologies. In general, raw microbiome sequence data are analyzed using a bioinformatics pipeline followed are motivated by sparsity and overdispersion issues. However, few statistical methods perform reliably for different microbiome datasets (Weiss et al., 2017). This paper proposes a novel unified zero-inflated Bayesian nonparametric (ZIBNP) approach that effectively meets the multipronged challenges of high-dimensionality, collinearity, sparsity, and covariate confounding to accurately detect the DA taxa across multiple groups of study samples. As its name suggests, the hierarchical model comprises a twocomponent mixture model; one mixture component is a point mass at zero, and the other component is a flexible Bayesian nonparametric model whose posterior adapts to the specific characteristics of the microbiome dataset. Technical and sampling zeros are cast in a missing data framework and stochastically inferred as a by-product of the posterior inferences. High dimensionality and collinearity (specifically, small n, large p issues) are mitigated by allocating the large number of taxa to fewer latent clusters defined by shared relative abundance patterns and induced by the Chinese restaurant process Lijoi and Prünster, 2010). The approach accounts for the available sample-specific covariates in a manner that relates the varying sampling depths to inferential precision while also conforming with the compositional interpretation of microbiome data.
We apply the ZIBNP technique to analyze two microbiome datasets and infer the set of DA taxa while adjusting for covariates: (i) The CAnine Microbiome during Parasitism (CAMP) study (Oliveira et al., 2017) explores the impact of natural parasitic infection on the gut microbiome of domesticated dogs. In addition to microbiome abundances, the data comprises animal-specific covariates and two groups of dogs determined by infection status (infected or uninfected); and (ii) Global Gut Microbiome data (Yatsunenko et al., 2012) on human subjects belonging to three geographical regions (groups).
The paper is organized as follows. Section 2 provides an overview of our proposed ZIBNP model. Section 2.1 develops the nonparametric aspects responsible for achieving dimension reduction, incorporating covariate effects, and detecting the DA taxa. Section 2.3 introduces the missing data framework for stochastically imputing technical zeros. Section 3 outlines the posterior inference and MCMC procedure for detecting DA taxa. The simulation study of Section 4, demonstrates the high accuracy achieved by the proposed method relative to some well-established analytical approaches capable of adjusting for covariates. Section 5 analyzes the motivating microbiome datasets using the ZIBNP method. Section 6 concludes with a brief discussion.

A Novel Bayesian Hierarchical
Approach. For n samples and p taxa, the microbiome data are arranged in an n × p matrix, with each cell representing the observed count or abundance of a taxon (column) in a particular study sample or subject (row). Without loss of generality, we assume that taxa corresponding to entire columns of n zeros were previously eliminated during preprocessing to give p columns of data with at least non-zero entry each. Denote the taxa abundance matrix by Z = ((Z ij )), where Z ij is the number of reads of taxon j in sample i, and Z i = (Z i1 , . . . , Z ip ) is the vector of taxa abundances for subject i.
For K ≥ 2, there are K groups determined by a categorical attribute partitioning the samples, e.g., disease status, disease subtype, or geographical region. The group membership of sample i is denoted by k i and n k samples belong to the kth group. In addition to the attribute defining the groups, there are T covariates, arranged in a matrix X = ((X it )). With L i representing sequencing depth or library size, the relative abundance vector for the ith subject, denoted byq i = (q i1 , . . . ,q ip ) , is computed asq i = Z i /L i .
We wish to identify the set of DA microbial taxa. To accomplish this end-goal, we foster a statistical framework capable of quantifying the complex relationships between microbial communities and the grouping factor of interest. In addition to the challenges posed by high dimensionality and sparsity, the association is confounded by sample-specific characteristics or covariates that produce spurious microbial associations (Vujkovic-Cvijin et al., 2020). We define a differential status variableh j , withh j = 2 (h j = 1) signifying that taxon j is (not) DA. The set of DA taxa is then D = j :h j = 2, j = 1, . . . , p , and posterior inferences about unknown parametersh 1 , . . . ,h p are of interest. As discussed in Section 1, biological zeros occur due to the actual absence of a taxon for a group, and manifest as blocks of taxa with zeros for an entire group of subjects. Therefore, these taxa are labeled as DA without further analysis. Technical zeros are related to issues such as batch effects (Wang and LêCao, 2019) often encountered during the processing and sequencing of samples. Sampling zeros randomly occur due to a taxon not being measured in a sample, e.g., because of small sampling depths. Analyzing the non-biological sources of sparsity is less straightforward and requires carefully designed statistical methods.
Consequently, our statistical framework primarily focuses on technical and sampling zeros. The proposed ZIBNP model relies on a two-component mixture for each element of Z, namely, a point mass at zero (representing technical zeros) and a Bayesian nonparametric (BNP) model for the counts, denoted by F , under which sampling zeros may stochastically occur. More formally, with F ij denoting the marginal distribution of Z ij under BNP model F , and I {0} representing a point mass at 0, the ZIBNP model assumes that where r ij is the probability of a technical zero in sample i and taxon j and R = ((r ij )) denotes the n × p matrix of probabilities of technical zeros. The second mixture component accounts for sampling zeros. Section 2.1 develops BNP model F . Section 2.2 details the stochastic process for technical zeros.
2.1. Bayesian nonparametric (BNP) model, F . The taxa abundances are modeled as where q i = (q i1 , . . . , q ip ) is a probability p-tuple, i.e., p j=1 q ij = 1. The columns of the unknown row-stochastic matrix Q = ((q ij )) are denoted byq j , j = 1, . . . , p. The marginal probability of a sampling zero under F is then Dimension reduction. Microbiome data usually involves far fewer samples than taxa, i.e., n << p, resulting in severe collinearity in the abundance matrix Z. Model F mitigates collinearity by detecting lower-dimensional structure in the matrix columns and reducing them to C unique motifs, where C is unknown and far smaller than p. Specifically, the model performs unsupervised clustering of the taxa into C latent clusters. For an allocation variable c j , let {c j = u} symbolize the random event that the jth taxon is allocated to the uth latent cluster. Let m u = p j=1 I(c j = u) be the number of taxa in cluster u, so that C u=1 m u = p. We associate the uth latent cluster with a common cluster-specific latent vector, q * u = (q * 1u , . . . , q * nu ) . The columns of matrix Q corresponding to the m u taxa are assumed to be identically equal to q * u for clusters u = 1, . . . , C. That is, (4)q j = q * cj , for all taxa j = 1, . . . , p. This formulation eliminates redundancy in the large number of taxa by detecting patterns in the matrix Q columns, with the latent vectors q * 1 , . . . , q * C representing the cluster-specific motifs or common across-sample patterns shared by all taxa belonging to a cluster. Then, conditional on the vector c = (c 1 , . . . , c p ), matrix Q is fully determined by Q * = ((q * iu )) of dimension n by C. In conjunction with model (2), this implies that all taxa j belonging to a cluster have similar, but not necessarily identical, relative abundancesq 1j , . . . ,q nj . Since Q is row-stochastic, each row of matrix Q * satisfies implying that cluster motifs q * 1 , . . . , q * C are dependent. Allocation vector c is given the semiparametric Chinese restaurant process (CRP) prior with precision or mass parameter α c  because of its ability to achieve dimension reduction (e.g., Medvedovic, Yeung and Bumgarner, 2004;Kim, Tadesse and Vannucci, 2006;Dunson and Park, 2008;Guha and Baladandayuthapani, 2016;Guha, Jung and Dunson, 2022). With [·] representing densities, we have where P p is the set of all partitions of p taxa into C non-empty clusters, and 1 ≤ C ≤ p is random. In general, CRPs achieve dimension reduction because the random number of clusters, C, is asymptotically equivalent to α c log(p) as p → ∞ (Lijoi and Prünster, 2010). Since matrix Q is a "less noisy" version of the relative abundance matrix, a consequence of model assumption (4) is that differential statuses are shared by all taxa belonging to a latent cluster, which is itself collectively DA or non-DA. Analogously to the taxon differential status, we envision cluster differential status variables h 1 , . . . , h C , with the DA taxa equivalently given by (7) D = j : h cj = 2, j = 1, . . . , p .
Incorporating covariates. We model the dependencies between the latent vectors q * 1 , . . . , q * C and the sample-specific covariates. First, we preselect and hold fixed a singleton cluster consisting only of a reference taxon; strategies for choosing this special taxon are discussed below. Without loss of generality, the reference cluster and reference taxon are labeled c = 1 and j = 1, respectively. Next, since microbiome abundance data are compositional, we model the log ratios x il β lkiu , σ 2 e , i = 1, . . . , n, and u > 1, where the latent elements are identifiable due to constraint (5). We assume that σ 2 e is small enough that most of the variability in the η iu 's is explained by the covariates. In the (unrealistic) situation in which the covariate effects are negligible, the latent elements q * iu of all samples i of a group are approximately equal.
The cluster-and group-specific regression parameters, β ku = (β 0ku , . . . , β T ku ) ∈ R T +1 , quantify the covariate relationships for samples belonging to group k and taxa belonging to non-reference clusters u = 2, . . . , K. Since the regression vectors β 1u , . . . , β Ku adjust for sample characteristics, the shared differential status of each non-reference cluster u and its member taxa are immediately available upon inspection of the regression vectors, as we discuss in the sequel.
Let X † = [1 n : X] be the matrix of dimension n by (T + 1), and let column vector ν1 M = (ν, . . . , ν) . With δ µ denoting a point mass at µ and Dir M (γ) denoting a Dirichlet distribution in R M for a parameter vector γ with positive elements, the (T + 1)-variate regression parameters are assigned a M -component finite mixture hierarchical prior: and τ 2 follows an inverse-gamma hyperprior with parameters a τ and b τ . The finite mixture prior is key because it allows the borrowing of strength between the latent clusters and groups through shared regression parameter vectors β ku . More specifically, membership variable v ku records the M -mixture component from which regression vector β ku is drawn. That is, Choosing the reference taxon Since they adjust for covariates, the group-specific regression coefficients in equation (8) are informative about the unknown cluster differential statuses. However, because of assumption (8), the differential status of the reference cluster is inextricably related to any criterion for calling the DA and non-DA clusters. Hence, by convention, we require that the reference taxon be non-DA. These are some options for choosing a reference taxon: • Minimum variance taxon Following Nearing et al. (2022), select the taxon j with the smallest variance of the relative abundancesq 1j , . . . ,q nj . Then reorder the taxa so that the reference taxon and its singleton cluster have the label 1. • Artificial reference taxon Augment abundance matrix Z with an artificial taxon with label 1 and unit abundance counts for all samples. Increase p by 1. Since the sampling depths are large, this additional "taxon" is likely to be non-DA and will not significantly change the analyses.
We have used the second option in the simulation study and data analyses.
Detecting differential statuses of taxa. The K(C − 1) regression coefficient vectors in expressions (8) and (9) adjust for covariate effects to reveal the DA clusters. Furthermore, as discussed, the finite mixture prior in (9) causes the regression vectors to be tied with positive probability. The common differential statuses of the non-reference clusters and their allocated taxa are then inferred in a straightforward member. Specifically, for non-reference clusters u = 2, . . . , C: In other words, a non-reference cluster is non-DA if and only if its K regression vectors are identical. By design, the reference cluster (taxon) is always non-DA. Further, the differential statuses of the taxa are identical as their parent clusters by (7). These ideas are demonstrated by a toy example in Figure 1, where we have n = 7 subjects belonging to K = 2 groups (adults and children), p = 50 taxa, and T subject-specific covariates. The taxa are allocated to C = 8 CRP clusters. Vector x † i ∈ R T +1 denotes the ith row of the n by (T + 1) matrix X † , and iu i.i.d.
∼ N (0, σ 2 e ). Cluster 1, containing only reference taxon 1, is not represented in the figure. For the seven non-reference clusters, K(C − 1) = 14 vectors β ku , each consisting of (T + 1) regression coefficients including the intercept, arise from a finite mixture model (9)   Cartoon illustration of the procedure for calling the taxa differential statuses. For n = 7 subjects, there are K = 2 groups (3 adults and 4 children), p = 50 taxa, T covariates, and C = 8 latent clusters. Cluster 1, containing only reference taxon 1, is non-DA and not shown here. The colors represent the components of a finite mixture model (FMM) from which the group-cluster regression coefficients, β ku ∈ R T +1 , for group k = 1, 2 and non-reference cluster u > 1 are drawn. The taxa allocated to clusters 2, 3, and 8 (i.e., taxa 2, 3, 6, 8, 10, 12, 23, and 50) are DA because, for all these clusters, the regression vectors of adults and children correspond to different FMM components or colors. Taxa 4, 5, 7, and 9, allocated to cluster 4, are non-DA because the cluster's regression vectors for adults and children, β 14 and β 24 respectively, are identical (blue). See the text for further details.
represent distinct components of the finite mixture model. For example, the regression vectors β 14 , β 24 , and β 28 are equal because they are drawn from the same mixture component (shown in blue). Applying expression (12), cluster 4 is non-DA, whereas clusters 2, 3, and 8 are DA. All taxa share the same differential attributes as their parent clusters, e.g., taxa 4, 5, 7, and 9 belonging to cluster 4 are all non-DA, whereas taxa 8 and 23 belonging to cluster 3 are DA.

Modeling technical zeros.
Returning to the first mixture component in expression (1), we model r ij , the probability of technical zeros in sample i = 1, . . . , n, and taxon j = 1, . . . , p. The log-sampling depth and covariates are important predictors of the proportion of technical zeros in a sample (Ma, Luo and Jiang, 2020). For random effects vector, λ iu = (λ 0iu , . . . , λ T +1,iu ) , accounting for cluster-specific effects shared by the taxa, we assume that ∼ N (0, τ 2 λ ) for all t, i, and u. In addition to the covariates, expression (13) includes log-sampling depth as a known predictor of the logit-likelihood of technical zeros (Jiang, Li and Li, 2021). Denote by Λ the collection of cluster-specific random matrices Λ 1 , . . . , Λ C , where Λ u = ((λ tiu )) is the n × (T + 2) matrix of regression coefficients for the uth cluster in equation (13). Because all m u taxa belonging to the uth cluster have the same probability of technical zeros for a given sample, we write r * icj = r ij , i = 1, . . . , n; j = 1, . . . , p. From expressions (1) and (3), and applying the CRP taxon-to-cluster allocations, the probability of a non-biological (i.e., either technical or sampling) zero in taxon j of sample i is Due to the lower-dimensional structure imposed by submodel F and equation (13), the abundance matrix columns of all taxa in a latent cluster display similar zero patterns. This feature of the ZIBNP model helps overcome the challenges of data sparsity in accurately detecting the latent clusters.
2.3. A missing data framework for technical zeros. BNP submodel F accounts for sampling zeros that are relevant to detecting DA taxa. By contrast, since technical zeros are related to sequencing process errors, they obfuscate the submodel F parameters in expression (1) and are not informative about the differential statuses of the taxa. Consequently, the ZIBNP model is cast in a missing data framework to make accurate inferences about the differential patterns of the groups.
Specifically, we interpret expression (1) as an independent mechanism that, with probability r * icj , converts the abundance count,Z ij , arising from submodel F to a technical zero. We denote the true abundance matrix corresponding to F byZ = ((Z ij )). Inferences about parameters directly relevant to the DA taxa rely on the partially latent matrixZ rather than the observed matrix Z.
Let δ ij represent the missingness indicator thatZ ij is missing, so that δ ij The observed abundance counts are then related to the true counts as If Z ij > 0, we clearly have δ ij = 1. On the other hand, Z ij = 0 may correspond to a sampling zero (δ ij = 1) or technical zero (δ ij = 0). Unlike standard missing data approaches, indicator δ ij is unknown for zero abundances, and the probability of a technical zero is the sampling zero probability given by the complementary event. Since the sampling depths are usually large, a technical zero is much more likely than a sampling zero, unless q * icj is very small. In any case, equation (15) can be applied to iteratively call the technical and sampling zeros (i.e., a posteriori generate the unknown missingness indicators) in the MCMC procedure outlined in Section 3.1. Let J i = {j : δ ij = 0, j = 1, . . . , p} be the taxa with missing counts in the ith sample. Using latent matrixZ, the observed sampling depths are then L i = j / ∈JiZ ij . The latent sampling depth relies on the missing taxa and is defined asS i = j∈JiZ ij . The true sampling depth is defined asL i = p j=1Z ij and equals (L i +S i ). With these definitions, the following result, whose proof appears in Supplementary Material, allows us to stochastically reconstruct the true abundance matrixZ in the posterior: THEOREM 2.1. Suppose the missingness indicators δ ij corresponding to the zero abundances are given. For sample i = 1, . . . , n, letq i = j∈Ji q * icj . Then 1. Latent sampling depthS i has a negative binomial distribution: 2. Let the vector of |J i | missing taxa abundances beZ Finally, we assign standard conjugate priors to the remaining hyperparameters. Figure 2 displays the directed acyclic graph (DAG) representation of the ZIBNP model.

Posterior Inferences.
3.1. MCMC algorithm. The model parameters are initialized using the naive estimation techniques described in Supplementary Material. Subsequently, all parameters are iteratively updated by MCMC procedures. All model parameters except η iu are updated by Gibbs samplers. Although the full conditional of η iu does not have a closed form, it is log-concave, and the η iu 's can therefore be generated by adaptive rejection sampling (Gilks and Wild, 1992); see Section 2 of Supplementary Material for an outline of the MCMC steps. Following burnin, the MCMC sample is stored and post-processed for posterior inferences.
3.2. Taxa differential statuses. First, we derive expressions for conditional posterior probabilities of the differential statuses of the latent clusters and their member taxa under the ZIBNP model. The following results follow from law of total probability and are omitted for brevity. THEOREM 3.1. For the ZIBNP model, let the symbol Θ − generically denote all model parameters not deterministically related to the quantities discussed in each the following statements. Thus, the set of parameters indicated by Θ − depends on the context. Given Θ − and Z, let P * [·] denote the conditional posterior probability, P · | Θ − , Z . Then 1. For clusters u = 1, . . . , C, the conditional posterior probability that the uth cluster is non-DA is P * [h u = 1] = M m=1 K k=1 P * v ku = m , where membership variable v ku is defined in (11). 2. For the non-reference taxa j = 2, . . . , p, the probability that the jth taxon is non-DA is , where the expression for P * [h u = 1] on the right hand side is given in Part 1 of the theorem.
The above theoretical expressions suggest MCMC-based empirical average estimators of the (unrestricted) posterior probabilities of the taxa differential statuses. For example, letP (l) [·] denote conditional posterior probability P · | Θ − (l) , Z , evaluated using the lth MCMC sample model parameters, and L be the MCMC sample size. Then, applying Part 2 of Theorem 3.1, an MCMC-based empirical average estimator of posterior probability The procedure is applied to infer the posterior probabilities of the differential statuses of the p taxa along with uncertainty estimates. For a nominal FDR (say, 5%), we apply the Bayesian FDR or direct posterior probability approach (Newton et al., 2004) to select an appropriate posterior probability threshold, κ, for calling the DA taxa. An estimate of differential status variableh j is where I(·) represents the indicator function. The set of inferred DA taxa is thenD = j : h j = 2, j = 1, . . . , p .
4. Simulation Studies. Using simulated microbiome datasets with different sparsity levels and effect sizes, we investigated the effectiveness of the proposed ZIBNP method, making comparisons with some existing statistical methods with available R packages and capable of detecting DA taxa in the presence of covariates.
Data generation. For n = 100 subjects belonging to K = 2 groups, we generated abundance matrix Z for p = 1, 002 taxa, with taxon 1 as the artificial reference taxon whose abundance counts, Z i1 , . . . , Z 100,1 , were all equal to 1. The following generation steps, which differ in several important aspects from the ZIBNP model, were implemented to generate the artificial datasets; for example, we relied on the SparseDOSSA  model, which assumes a log-normal distribution for the marginal feature abundances and a multinomial distribution for the counts.
Thirty datasets were generated in each of nine scenarios corresponding to the combinations of (i) three sparsity levels (15%, 25%, and 58%) matching typical microbiome data and regulated by simulation parameter λ 0 in Step 6 below, and (ii) three effect sizes quantified by the fold change, ω, equal to 2, 3 and 4 in the SparseDOSSA model of Step 4 below. The specific details are as follows: 1. Covariates We used the covariates from a publicly available oral microbiome dataset (Burcham, et al. and Taste Lab Citizen Scientists, 2020) consisting of adults and children as the subject groups. We selected T = 4 covariates, namely, three binary covariates (antibiotics taken during the last 6 months, brush teeth daily, and sex) and one continuous covariate (BMI). Randomly sampling the covariates of 50 subjects from this real dataset, we repeated the same set of covariates in adults and children to obtain the covariate matrix, X, of dimension 100 by 4. More formally, if X (0) denotes the 50 × 4 covariate matrix sampled from the real oral microbiome dataset, then where the first 50 subjects represent adults, and the remaining 50 subjects represent children. 2. Allocation vector We generated the true number of clusters, C, from a discrete uniform distribution on {8, 9, . . . , 20}. Applying the rpartitions package in R, we generated a C-partition of p objects using the function rand_partitions. The true taxon-to cluster allocation variables, c 1 , . . . , c p , were set equal to this randomly generated partition with C clusters, and therefore relied on an entirely different stochastic mechanism than the CRP of the ZIBNP model. 3. Cluster differential statuses With h u = 2 signifying a DA uth cluster, the differential statuses were generated as h u i.i.d.
∼ 1 + Bernoulli(π 0 ), u = 1, . . . , C, setting π 0 = 0.3. 4. Regression coefficients Using the T = 4 covariates of Step 1, the cluster-specific regression coefficients excluding the intercepts were generated as . . , C. The regression intercepts, ψ ku , of the 2C group-cluster combinations determined the effect sizes, and were evaluated as follows: (17) ψ ku = 6 + log(ω)I k = 2, h u = 2 , k = 1, 2, to obtain the regression vectors of the group-cluster combinations as β ku = (ψ ku , β − u ) T . Construct (17) ensures that the observed counts generated in the subsequent steps tend to be systematically different in the two groups for the DA clusters. 5. True taxa abundances Using the sparseDOSSA model, the non-missing absolute abundances were generated using log-normal distribution: . . , n, j = 1, . . . , p, with σ 2 chosen so that R 2 = 0.97 in the linear regression analysis of {log A ij }. To mimic the range of sampling depths observed in actual microbiome datasets, we generated the observed sampling depths,L i i.i.d.
Varying simulation parameter λ 0 in the set {−0.1, −.059, .023} produced the required range of sparsity levels in the simulation scenarios; see Table 1. In the simulation study, comparative benchmarking of the competing statistical methods in the nine scenarios corresponding to the combinations of sparsity level 15%, 25% and 58%, and fold change 2, 3, and 4.
For the 30 artificial datasets in each of the nine scenarios, and disregarding knowledge of the generation mechanism and its true parameter values, the Section 3 procedure was applied to analyze each artificial dataset and make posterior inferences for the proposed ZIBNP model. Applying the Section 3.2 strategy, we post-processed the MCMC sample to estimate the DA posterior probabilities, P h j = 2 | Z . As described in Section 3.2, assuming a nominal FDR of 5%, we applied the Bayesian FDR procedure to call the DA taxa, and thereby, evaluated the achieved false discovery rate (FDR), Matthews correlation coefficient (MCC), sensitivity, and specificity of ZIBNP for each simulated dataset. Unlike other performance metrics, MCC delivers a holistic evaluation of a method's accuracy by producing a high score near 1 only if a method's predictions are accurate with respect to all four confusion matrix metrics; namely, true positives, false negatives, true negatives, and false positives (Chicco and Jurman, 2020).
We compared our technique with some well-established methods for differential abundance analysis with covariates, focusing only on methods implemented in publicly available R packages: 1. ANCOMBC (Lin and Peddada, 2020): We used the ancombc function of R package ANCOMBC (1.0.5) with default settings. The function fits compositional data by making a log transformation and using a sample-specific offset term for bias correction. 2. DESeq2 (Love, Huber and Anders, 2014): We used the DESeq2 function of package DE-Seq2 (1.30.1) with argument sizefactor set to "poscounts". The R function performs a likelihood ratio test for a reduced model containing all covariates except the grouping factor. 3. Metagenomeseq (Paulson et al., 2013): R function fitZIG of the metagenomeSeq (1.32.0) package was applied to fit a zero-inflated Gaussian version of the technique. The workflow involved normalizing the data using cumulative sum scaling (Paulson et al., 2013) with default settings. 4. Maaslin2 : The Maaslin2 function of the Maaslin2 (1.4.0) R package was used with total sum scaling normalization and covariate fixed effects. This technique fits a generalized linear model to each feature abundance with respect to the covariates, checking for significance using the Wald test.
The competing statistical methods were corrected for multiple hypothesis testing using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) to obtain taxa-specific adjusted p-values (Storey, 2003). Using the adjusted p-values, for a target FDR of 5%, the DA taxa were called using the Benjamini-Hochberg procedure to evaluate each method's achieved FDR, MCC, sensitivity, and specificity in each dataset.
Summarizing over the 30 datasets of each simulation scenario, Figure 3 displays boxplots of different performance metrics of the competing methods in detecting the DA taxa. We find that fitZIG displays high sensitivity, but has highly inflated FDR, low MCC, and low specificity in all the scenarios. MaasLin2 exhibits inflated FDR and low MCC. It is interesting that the FDR of MaasLin2 increases with increase in fold change, an anomaly also reported by other studies (e.g., . ANCOMBC successfully controls FDR at lower sparsity levels, but is unable to reliably detect DA taxa, as evidenced by its low sensitivity in all nine scenarios. Although DESeq2 performs reliably for higher levels of fold change, its performance significantly deteriorates with smaller effect sizes. From Figure 3, we observe that ZIBNP is generally reliable with impressive sensitivity, MCC, and specificity. Its FDR is well controlled for lower sparsity levels. Although FDR is somewhat affected for the highest sparsity level (58% zeros), ZIBNP nonetheless performs well even in these scenarios, exhibiting high MCC, sensitivity, and specificity. In summary, ZIBNP appreciably outperforms the methods DESeq2, ANCOMBC, fitZIG, and MaAsLin2 in inferential accuracy and is reliable for analyzing even sparse microbiome datasets.
5. Data Analysis. We applied the ZIBNP technique to analyze publicly available microbiome data from the CAMP and Global Gut studies, consisting of K = 2 and K = 3 groups, respectively. The results were caffected with other differential analysis techniques.

CAMP study canine data.
The dataset was downloaded from the MicrobiomeDB resource (Oliveira et al., 2017). To investigate the association between eukaryotic parasite infection on the composition of the gut microbiome, the data was obtained by sequencing the V4 region of the 16S rRNA gene of the fecal samples of 155 infected (case) and 115 uninfected (control) dogs. Animal-specific attributes such as sterilization, pet ownership, age, and sex are available, allowing statistical methods for differential analysis to adjust for these covariates while detecting DA taxa. Table 2 presents the DA taxa detected by ZIBNP between the case and control groups . Some of the detected taxa have been reported by previous studies. For example, bacteria of the genus Bacteroides were detected as DA with low relative abundance in infected dogs, and bacteria of the genera Megamonas and Prevotella were significantly associated with the infection status of dogs in studies of the effect of parasite-induced infections on gut microbiome composition (Berry et al., 2020;Pilla and Suchodolski, 2020). The genus Blautia has been previously reported to have a strong dysbiosis in dogs with acute diarrhea (Pilla and Suchodolski, 2020). In addition, microbes from the genera Streptococcus (Pilla and Suchodolski, 2020), Ruminococcus gnavus (Hall et al., 2017;Boucard et al., 2021), and Alloprevotella (Boucard et al., 2021) have known associations with inflammatory bowel disease in humans and dogs.
The methods ANCOMBC, fitZIG, Maaslin2, DESeq2, and ZIBNP detected 9, 72, 7, 42, and 10 DA taxa, respectively. Table 3 displays the pairwise similarity measures using the Jaccard index (Jaccard, 1901) between the DA taxa detected by the different statistical methods, all of which attempt to adjust for covariates. The Jaccard index (J) is a measure of similarity between two sets that ranges between 0% to 100%. It is defined as the size of the intersection divided by the size of the union between two sets; J = 0% implies that the two sets have no overlap, whereas J = 100% implies the sets completely overlap. Although all the methods in Table 3 display low overlap, there is moderate overlap between ANCOMBC and Maaslin2, possibly because both methods fit generalized linear models to the data. The proposed ZIBNP method has some shared DA taxa with ANCOMBC and Maaslin2. However, ZIBNP has a low overlap with fitZIG and DESeq2, respectively relying on parametric zero-inflated Gaussian and negative binomial models. This reveals that the findings of ZIBNP are substantially different from existing statistical methods. The genera Blautia and Lachnoclostridium were the common findings of the methods ZIBNP, ANCOMBC, and Maaslin2. Additionally, the genus Ruminococcus gnavus was detected by ZIBNP as well as ANCOMBC.

5.2.
Global Gut Microbiome study. We applied the proposed ZIBNP method to analyze the motivating Global Gut Microbiome data (Yatsunenko et al., 2012). The study examines the differences in the microbial abundance between samples collected from individuals residing in Malawia, Venezuela, and the USA. More specifically, the data consists of 100 US individuals, and 83 individuals each from Malawi and Venezuela, in addition to age and sex as covariates. The abundance matrix of 266 subjects × 1, 270 taxa was comprised of 30% zeros. Along with the proposed ZIBNP approach, we analyzed the data using the methods ANCOMBC, fitZIG and DESeq2, all of which perform differential analysis with covariates. Since MaAsLin2 in its current form does not have a global significance test for more than two groups, we did not use this method to analyze the Global Gut data.
The results are shown in Table 4 and Figure 5. The four methods detected a common set of 37 DA taxa, as seen in Figure 5. The low overlap of detected DA taxa between ZIBNP and other competing methods reveals the benefits of the proposed nonparametric Bayesian approach in finding otherwise undetected taxa. ANCOMBC and fitZIG shared a large proportion of DA taxa and had a Jaccard index of 82%. FitZIG (MetagenomeSeq) detected 1, 035 (out of 1, 270) taxa as DA, and in fact, detected the largest number of DA taxa in both motivating datasets. However, as noted by several other studies (Thorsen et al., 2016;Hawinkel et al., 2017) and Section 4, the methods fitZIG and ANCOMBC often exhibit inflated FDR, and one should proceed with caution while interpreting their results. Table S1 of Supplementary Material provides the list of 201 DA taxa detected by ZIBNP. Among all DA taxa detected by ZIBNP, 5.5% and 5% taxa belonged to the genus Prevotella and Bacteroides, respectively. The original article (Yatsunenko et al., 2012) reporting the Global Gut Microbiome results focused on pairwise comparisons, such US versus non-US  individuals and Venezuela versus Malawi individuals. They are not, therefore, directly comparable to the multigroup comparisons of ZIBNP. However, Yatsunenko et al. (2012) reported several taxa belonging to the genus Prevotella as differential for the (pairwise) geographical regions. Furthermore, they stated that taxa belonging to the genus Bacteroides were significantly more abundant in US individuals compared to non-US individuals.
6. Discussion. Differential abundance analyses between multiple groups of study samples help identify novel therapeutic targets for disease treatment. However, very few existing methods for differential analysis with covariates perform consistently well for different datasets (Weiss et al., 2017). Motivated by the challenges of high dimensionality, sparsity, and compositionality inherent to microbiome data, we propose a unified methodology based on a two-component mixture model that adjusts to distinctive data characteristics. A key contributor to our strategy's success is a model-based missing data framework that relies on the stochastic relationship between latent sampling depths and the large number of zero counts to improve the accuracy of differential abundance analysis.
Latent clusters for the taxa induced by the semiparametric Chinese restaurant process alleviate collinearity issues caused by the small n, large p problem. The compositional aspects of microbiome data are modeled by a regression framework for the log-ratios of the dependent cluster-specific probabilities. The taxa differential statuses are immediately available from the group-cluster regression parameters, whose common multivariate prior has a finite number of atoms. Through simulation experiments and data analyses, we demonstrate the potential of the proposed ZIBNP method as a reliable technique for detecting DA taxa in a wide variety of microbiome datasets.
As suggested by several studies (Sankaran and Holmes, 2014;Xiao, Cao and Chen, 2017), the accuracy of statistical methods for differential analysis can be further improved by incorporating the phylogenetic distances between the taxa. Going forward, our research will focus on utilizing this potentially valuable information. We are also developing an R package for ZIBNP that will be made available on GitHub for differential analysis of multiple groups with covariates.
Funding. The third author (SG) was supported by the National Science Foundation under Award DMS-1854003.

SUPPLEMENTARY MATERIAL
MCMC implementation Section 1.1 describes the MCMC initialization for the ZIBNP model parameters. Section 1.2 outlines the MCMC updates.
Proof of Theorem 2.1 Section 2 of the supplementary material describes the proofs of Parts 1 and 2 of Theorem 2.1.

Table S1
Excel file listing the DA taxa detected by ZIBNP for the Global Gut Microbiome study.