CAT: a conditional association test for microbiome data using a leave-out approach

The development of next-generation sequencing techniques has enabled high-resolution profiling of the human microbiome. The challenges of analyzing microbiome data include its high dimensionality and the structural relatedness of the observed features. As a starting point in assessing the link between the microbiome and an outcome, researchers often test the global association of a phenotype with the microbiome as a whole. Global association tests addressing this question typically employ microbiome-specific metrics (also called beta diversity measures) that are more robust to extreme values than classical Euclidean distances. Popular choices include Bray-Curtis dissimilarity (Bray and Curtis, 1957), which is designed for count data, and weighted and unweighted UniFrac (Lozupone and Knight, 2005; Lozupone et al., 2007), which incorporate information on the location of features in a phylogenetic tree.

For microbiome datasets with clearly identified sample groups, a popular global association test is PERMANOVA, which utilizes permutation testing to obtain a $p$-value for the null hypothesis that there is no difference in the location of the centroids across groups (Anderson, 2001). PERMANOVA has been widely used in microbial analysis, as it is simple to apply and only requires observations of the outcome variable and the pairwise distances or dissimilarities between samples. Moreover, recent adaptations of PERMANOVA can handle nested designs and correlated outcomes (Oksanen et al., 2022).

Another popular global association test for microbiome data is MiRKAT, which is based on kernel machine regression (Zhao et al., 2015). One advantage of this method is that it can incorporate multiple candidate distance metrics to maximize power for a particular data set. In addition to continuous and binary response variables, MiRKAT has been extended to handle survival outcomes (Plantinga et al., 2017) and correlated or dependent samples (Zhan et al., 2018; Koh et al., 2019). These global association methods have been widely applied in microbiome studies. However, they cannot provide inference on specific taxa.

Thus, a practical question in following up on a significant global association test result is to identify specific microbiome features that drive the global testing result. However, existing methods for testing the effect of individual taxa often ignore the presence of related features and focus only on the marginal effect of the individual feature. In particular, popular differential abundance methods such as ALDEx2 (Fernandes et al., 2014), DESeq2 (Love et al., 2014), and ANCOM-BC (Lin and Peddada, 2020) all adopt a marginal testing framework.

An inherent problem in marginal testing of microbiome data is nested discoveries, where hits are linked within a taxonomic or phylogenetic hierarchy. Taxonomic trees reflect the traditional labeling and organization of microorganisms into groupings such as family or genus, while phylogenetic trees reflect evolutionary history, with branch points corresponding to events that gave rise to differences in the genomic sequences. Both types of trees play a key role in understanding microbiome data. Taxonomic labels serve as a basis for interpretation since they are standardized across studies. Phylogenetic trees are useful in analysis, as they encode rich information on sequence similarity, which drives phenotypic and functional similarity.

The relatedness among features can make it challenging to pinpoint which taxa play a critical role in influencing outcomes. For example, when a genus is found to be significant, the corresponding higher taxonomic units to which it belongs, such as family and order, also tend to be significant. However, the precise taxonomic level most relevant to the outcome is difficult to establish. A conditional test can provide direct quantification of the importance of a specific feature in contributing information not captured by other features in the data set, addressing a question with potentially greater biological importance than a marginal test.

Notably, the challenge of correlated predictors exists in many high-dimensional datasets, yet is particularly prominent in microbiome data. In other settings, researchers have proposed rigorous definitions of feature-outcome independence. Here, we adopt Candès et al. (2018)’s definition, where a feature is said to be “null” if and only if the outcome is independent of it conditionally on all other variables.

In this paper, we present a novel conditional test, CAT, that provides a natural next step to follow up on a significant global association test result. CAT achieves the goal of assessing the importance of individual taxa while accounting for phylogenetic structure and other features in the data set. The remainder of the paper is organized as follows: Section 2 describes the proposed CAT method in detail. Section 3 demonstrates its performance using simulated data, and Section 4 illustrates the method through applications to real datasets with binary and survival outcomes. Section 5 concludes the paper with a discussion.

We begin with an illustration highlighting the motivation behind our approach. In the left panel of Figure 1, we show a PCoA plot based on weighted UniFrac distance that depicts variation in microbiome composition between melanoma patients that responded to immunotherapy vs. those that did not. The 95% confidence regions for each group indicate there may be global differences in the microbiome profiles between the groups. In the right panel, we artificially removed counts belonging to the family Ruminococcaceae; the two clusters of patients become less separated in the PCoA plot, with a reduced distance between centroids. Removing Ruminococcaceae from the data weakened the global association between the microbiome and response, suggesting that Ruminococcaceae may play an important role in driving the global association results. In the remainder of this section, we describe our proposal for a formal testing procedure aimed at quantitatively describing this phenomenon.

In our proposed approach, we start with the finest resolution features, corresponding to the leaf nodes in the taxonomic tree. For microbiome data derived from profiling of the 16S rRNA gene, these features are typically defined as Amplicon Sequence Variants (ASVs) or operational taxonomic units (OTUs). Our method may be applied as well to features derived from whole metagenome sequencing (WGS). By comparing the representative sequence for each feature against an established reference library, the feature can be assigned a taxonomic classification. Taxonomic levels from broad to specific follow the sequence kingdom, phylum, class, order, family, genus, and species. Based on the taxonomic assignments, one can draw a taxonomic tree reflecting the relatedness of all the features in the data set.

Most existing methods for identifying individual feature associations from WGS or 16S data focus on marginal associations and do not quantify how much individual features contribute to the results from global association testing. To fill this gap, we propose CAT, which tests the association between specific features and outcomes while conditioning on the tree structure and the abundance of other features in the tree. The CAT method is rooted in the coefficient of determination $R^{2}$ for global microbiome association tests. CAT estimates the change in $R^{2}$ for a global test using the original dataset vs. a modified data set with the taxon of interest removed, and obtains a $p$-value through a bootstrap procedure, which entails sampling from the data with replacement (Efron and Tibshirani, 1994). One widely used global association test is the nonparametric PERMANOVA method (Anderson, 2001, 2017). We begin by briefly reviewing this test, which serves as a starting point for our proposed method.

We now describe how to apply CAT to test the conditional association between the outcome of interest and a specific taxon. Let $X$ denote the outcome vector for $n$ observations. Let $\mathbf{Z}$ represent the $n\times m$ matrix with the observed counts for the finest-resolution microbiome features, which correspond to the leaf nodes in a taxonomy tree $\mathcal{T}$ with $m$ leaves. We denote the set of leaf nodes for the full tree $\mathcal{L}(\mathcal{T})=\{1,\ldots,m\}$. For any internal node in the tree $t$, we let $\mathcal{L}(t)$ denote the leaf nodes corresponding to its descendants. Given these definitions, we lay out the steps of the CAT procedure as follows:

1. 1.

   Calculate the $n\times n$ sample pairwise distance matrix $\mathbf{D}$ for the original data matrix $\mathbf{Z}$.

2. 2.

   Perform PERMANOVA using $\mathbf{D}$ and $X,$ and obtain a coefficient of determination $R^{2}$ for the outcome of interest.

3. 3.

   For the specific taxon $t$ being tested by CAT, generate a new data matrix $\mathbf{Z}^{*}$ by converting all the elements of $\mathcal{L}(t)$ to have $0$ counts.

4. 4.

   Calculate a new pairwise distance matrix $\mathbf{D}^{*}$ using the modified data matrix $\mathbf{Z}^{*}$.

5. 5.

   Perform PERMANOVA using $\mathbf{D}^{*}$ and get a new coefficient of determination $R^{2*}$ for the outcome of interest.

6. 6.

   Perform bootstrap sampling, selecting the matching sample from the pairwise distance matrix $\mathbf{D}$, the modified distance matrix $\mathbf{D}^{*}$, and the outcome $X$ for $B$ bootstrap samples. For each sample, compute the coefficients of determination for the original and modified distances, $R^{2}_{(1)},R^{2}_{(2)},\dots,R^{2}_{(B)}$ and $R^{2*}_{(1)},R^{2*}_{(2)},\dots,R^{2*}_{(B)}.$

7. 7.

   Compute the differences between the original $R^{2}$ and the leave-taxon-out $R^{2*}$ for all $B$ bootstrap samples. The estimated $p$-value is the proportion of the $R^{2}$ differences that are less than zero:

   $$\hat{p}=\frac{1}{B}\sum_{i=1}^{B}I(R^{2}_{(i)}-R^{2*}_{(i)}<0).$$

In this section, we illustrate the utility of CAT as a follow-up to global testing and compare results from CAT with those from existing marginal testing approaches on simulated data. We develop realistic simulation scenarios by starting from a real microbiome data set, which will be examined more closely in Section 4. We adopt a ”spike-in” method to control the taxa driving the cross-group differences.

We now discuss the use of CAT to analyze two real data sets examining the role of the microbiome in shaping melanoma patient outcomes. First, we consider the data set from Gopalakrishnan et al. (2018), which dealt with the association of the microbiome to immunotherapy response. In Section 4.2, to illustrate the use of CAT for time-to-event outcomes, we apply CAT with MiRKAT-S to the study of Spencer et al. (2021), which characterized the role of the microbiome in shaping progression free survival.

To date, most existing methods for microbiome association testing have focused on either global or marginal testing. Here, we adopt a conditional testing framework, proposing the CAT method as a conditional association test using a leave-out approach. The leave-out idea is one of the most fundamental ideas in statistical testing, whose applications range from likelihood ratio testing to type III sum-of-squares analysis. CAT combines this classic idea with the flexibility of using various metrics and testing approaches designed for microbiome data. It is worth mentioning that though microbiome data motivate us to propose CAT, the method is widely applicable to a wide range of situations where non-Euclidean pairwise distances are used. In this paper, we only illustrate how to test the conditional association for one taxon; testing the effect of several taxa as a unit is also possible by changing leave-one-out to leave-multiple-out in Step 3 of the procedure.

Although our simulation results show that CAT may identify fewer features than marginal tests, particularly at lower levels in the tree, we do not directly address multiple testing in this paper. Commonly used multiplicity adjustments, such as the Bonferroni procedure or Benjamini-Hochberg procedure can be applied to $p$-values generated by CAT. However, the null hypotheses under testing in microbiome data sets are not independent. False discovery rate control in correlated conditional tests, particularly in the presence of a phylogenetic tree structure, is an area that we plan to address in the future.

The results from the CAT method have clear real-world relevance. There is a growing effort to develop interventions aimed at reshaping the microbiome by administering ”rationally designed” mixtures of bacterial strains (van der Lelie et al., 2021) which have been selected to confer potential benefit to the patient. Our method can identify features that are potentially influential conditioning on the existence of other bacteria. This will help clinicians identify intervention targets more efficiently.

KAD is partially supported by NIH P30CA016672, SPORE P50CA140388, CCTS TR000371 and CPRIT RP160693. CBP is partially supported by NIH R01 HL158796 and NIH/NCI CCSG P30CA016672. RRJ is partially supported by NIH R01 HL124112 and NIH R01 HL158796. YS is partially supported by NSF DMS 2310955.