MiRKAT-S: a community-level test of association between the microbiota and survival times
- Anna Plantinga^{1},
- Xiang Zhan^{2},
- Ni Zhao^{3},
- Jun Chen^{4},
- Robert R. Jenq^{5} and
- Michael C. Wu^{1, 2}Email author
Received: 6 September 2016
Accepted: 31 January 2017
Published: 8 February 2017
Abstract
Background
Community-level analysis of the human microbiota has culminated in the discovery of relationships between overall shifts in the microbiota and a wide range of diseases and conditions. However, existing work has primarily focused on analysis of relatively simple dichotomous or quantitative outcomes, for example, disease status or biomarker levels. Recently, there is also considerable interest in the relationship between the microbiota and censored survival outcomes, such as in clinical trials. How to conduct community-level analysis with censored survival outcomes is unclear, since standard dissimilarity-based tests cannot accommodate censored survival times and no alternative methods exist.
Methods
We develop a new approach, MiRKAT-S, for community-level analysis of microbiome data with censored survival times. MiRKAT-S uses ecologically informative distance metrics, such as the UniFrac distances, to generate matrices of pairwise distances between individuals’ taxonomic profiles. The distance matrices are transformed into kernel (similarity) matrices, which are used to compare similarity in the microbiota to similarity in survival times between individuals.
Results
Simulation studies using synthetic microbial communities demonstrate correct control of type I error and adequate power. We also apply MiRKAT-S to examine the relationship between the gut microbiota and survival after allogeneic blood or bone marrow transplant.
Conclusions
We present MiRKAT-S, a method that facilitates community-level analysis of the association between the microbiota and survival outcomes and therefore provides a new approach to analysis of microbiome data arising from clinical trials.
Keywords
Survival data Kernel machine regression Distance-based analysis Community-level analysisBackground
The human microbiota, or the collection of microorganisms that inhabits the human body, plays an important role in many areas of health and disease. The development of next-generation sequencing technologies allows culture-free profiling of entire microbial communities, often by sequencing the 16S ribosomal RNA (rRNA) gene [1–3]. Similar 16S sequences are clustered into operational taxonomic units (OTUs); when clustered at the level of 97% similarity, these OTUs represent bacterial species [4]. Using these methods, the human microbiome has been studied at a variety of body sites including the gut [5], skin [6], and respiratory tract [7] and has been associated with many health conditions, such as inflammatory bowel disease, diabetes, psoriasis, and chronic obstructive pulmonary disease [8–10].
Associations between the human microbiota and health outcomes can be assessed by comparing individual OTU abundances or overall diversity metrics between samples or conditions [5, 11]. However, since taxonomic profiles are sparse and high-dimensional—hundreds to thousands of unique OTUs may be identified, many of which are present in only a subset of samples—comparisons on the level of individual OTUs may have low power. An alternative to OTU-level analysis is to compare the microbiota at the community level, i.e., to compare overall taxonomic profiles between individuals [12–14]. This class of analyses is often performed by computing pairwise distances between communities (samples), where the distance metrics are ecologically relevant and may incorporate phylogenetic structure. The matrix of pairwise distances is summarized by its top principal coordinates for visualization, and distance-based multivariate methods coupled with permutation are used to determine if dissimilarity is related to an outcome. Community-level analyses may provide power gains by utilizing phylogenetic information, avoiding the multiple testing problem, and aggregating modest effects across multiple taxa [15].
Recently, as an alternative to distance-based approaches that use permutation analysis, Zhao et al. [14] proposed the microbiome regression-based kernel association test (MiRKAT). MiRKAT uses a kernel machine framework with a variety of ecologically informative kernels to test for associations between the human microbiota and either continuous or binary outcomes. Intuitively, MiRKAT compares similarity in taxonomic profiles between communities (where similarity is measured via a kernel, which can be obtained by transforming relevant distance matrices) to similarity in outcome measures. p values are obtained analytically using a variance-component score test. MiRKAT is equivalent to distance-based analysis but has the added advantages of flexible modeling of the relationship between the microbiota and outcome measures, natural incorporation of covariates, and efficient computation of p values.
A limitation of existing community-level analysis approaches is that they cannot accommodate censored survival outcomes. However, such outcomes are of tremendous interest as microbiome profiling studies move into the clinical arena. For example, the lung microbiota has been related to progression of idiopathic pulmonary fibrosis [16] and the gut microbiota to overall survival after allogeneic blood and bone marrow transplant [5]. Additional OTU-level studies with survival outcomes have shown associations between the intestinal microbiota and development of atopic dermatitis [17] and allergic rhinitis [18] in children.
To address this critical gap in the literature, in this paper, we propose a test for association between the microbiota and censored survival outcomes (MiRKAT-S), accounting for covariates and potential confounders. We perform a distance-based analysis using the kernel machine Cox regression framework, encoding taxonomic profiles into kernel matrices via a transformation of distance metrics appropriate for microbial communities. This allows the analysis to take into account phylogenetic information and other features specific to biological communities. To formally test the association between the microbiota, as encoded in the kernel matrix, and censored survival times, we use a variance-component score test. However, when applied to microbial community profiles summarized by common kernels, the usual test statistic with p values calculated by resampling procedures is highly conservative [19, 20]. We therefore implement a small sample correction that provides proper control of type I error while maintaining adequate power, and we calculate p values analytically rather than by resampling. We demonstrate the performance of MiRKAT-S using real and simulated data summarized by a variety of kernels commonly used in microbial ecology.
This work represents the translation of existing methods in genetic studies with survival outcomes to applications in microbiome research. The first major contribution of this work is to allow survival outcomes in the kernel machine regression framework with kernels that appropriately encode microbiome data. Our small sample correction method provides proper control of type I error and improved power when using microbiota-appropriate kernels, whereas the kernel machine regression-based test as implemented for genetic studies has almost no power to detect relationships between microbial taxonomic profiles and survival. Secondly, the ability to perform the test using a variety of kernels provides robustness to the nature of the true association between the microbiota and survival. Therefore, although MiRKAT-S is technically similar to previous kernel machine regression methods, it enables microbiome analyses that are not possible using existing methods.
Methods
To associate the microbiota at the community level and censored survival times, we will relate censored survival times to taxonomic profiles using a flexible non-parametric modeling framework. We will assess significance via a variance-component score test which acknowledges the modest sample sizes of the most taxonomic profiling studies. In this section, we first describe the modeling framework, followed by the testing strategy and technical advances necessary to ensure proper control of type I error in this framework. Finally, we describe simulations encompassing a variety of true relationships between the microbiota and survival time.
Model specification
Suppose that for each of n subjects, we observe the microbial taxonomic profile, encoded by a q-vector of OTU counts Z _{ i }, and a p-vector of other covariates X _{ i }. Let T _{ i } be the survival time and C _{ i } the censoring time for the ith subject. We observe the bivariate vector (Δ _{ i }, U _{ i }), where U _{ i }=min(T _{ i },C _{ i }) is the observed time and Δ _{ i }=I(T _{ i }≤C _{ i }) is the event indicator for subject i. We wish to test whether taxonomic profiles Z are associated with survival time, adjusting for covariates X.
where λ _{0}(t) is the baseline hazard function. In the kernel machine regression framework, f(·) is generated by a positive definite kernel function K(·,·), that is, f(·) lies in the reproducing kernel Hilbert space \(\mathcal {H}_{K}\). Under the representer theorem [22], \(f(\boldsymbol {Z}_{i}) = \sum _{i^{\prime }=1}^{n} \alpha _{i^{\prime }} K(\boldsymbol {Z}_{i}, \boldsymbol {Z}_{i^{\prime }})\) for some α _{1},…,α _{ n }. Choosing different kernel functions K(·,·) allows specification of a wide variety of models. For example, the kernel function f(Z _{ i })=Z i′γ, corresponding to the linear kernel \(\phantom {\dot {i}\!}K(\boldsymbol {Z}_{i}, \boldsymbol {Z}_{i^{\prime }}) = \boldsymbol {Z}_{i} \boldsymbol {Z}_{i^{\prime }}^{\prime }\), is used to specify a linear model. Kernels are similarity matrices, so each element K _{ j,k }=K(Z _{ j },Z _{ k }) represents the pairwise similarity between samples j and k. Because we use a score test, which depends only on the null model, any kernel will result in a valid test; however, the choice of kernel does affect the power of the test.
Score test
Thus far, we have assumed that there are no tied survival times. In practice, tied survival times are fairly common due to coarse time measurements resulting from specific visit schedules or study follow-up dates. We use the Efron approximation to accommodate tied survival times [29]. This approximation performs well even with relatively small sample sizes or a high proportion of ties [30].
Small sample correction
Specifically, the fitted kernel machine Cox model (Eq. 1) is equivalent to a weighted linear model at convergence with weight matrix estimated using an iteratively reweighted least squares (IRLS) algorithm, as described in Appendix 1. We use a diagonal approximation for both the covariance of the residuals \(\hat {M}\) and the weight matrix W for IRLS. Several versions of W have been used for weighted partial least squares in the literature (e.g., [34] and [32]); we use an intermediate version that is diagonal as in [32] and [33] but whose elements are defined by the negative Hessian with respect to β as in [34].
Simulation scenarios
We carried out simulation studies in a range of settings to confirm that MiRKAT-S properly controls type I error and to assess its power using a variety of kernels. Microbiome OTU counts were generated using the same approach as [14]. Specifically, for each individual, we simulated OTU counts from a Dirichlet-multinomial distribution with dispersion parameters and proportions estimated from Charlson et al.’s real upper respiratory tract microbiome dataset, in which 856 OTUs were measured on each of 60 individuals [7]. The data for each simulated individual consists of 1000 total OTU counts distributed among the 856 OTUs of Charlson et al. We also simulated two covariates for each individual, X _{1i } and X _{2i }, from a standard normal and a Bernoulli (0.5) distribution independently of taxonomic profiles. We considered sample sizes ranging from n=25 to n=500 individuals. For all simulation scenarios, we generated datasets with approximately 25 and 50% censoring. Four simulation settings were considered, varying (1) whether OTU abundance or the presence/absence of OTUs was associated with the outcome and (2) whether phylogenetically clustered or unclustered OTUs were associated with the outcome.
where \(\bar {Z}_{(j)}\) is the average across samples of the counts for the jth OTU. This limits the ability of a single OTU to dominate the communal effect of the microbiota. Setting 2 is comparable to setting 1 when the associated cluster is common, since in both cases the abundance of common OTUs is associated with survival times, but it lacks setting 1’s close phylogenetic relationship between associated OTUs.
As in setting 1, we simulated situations where an abundant cluster, containing 19.7% of all reads, was associated with the outcome and where a rare cluster, containing 0.9% of all reads, was associated with the outcome.
Finally, in setting 4, the presence or absence of 40 randomly selected OTUs was associated with the survival time. This mimics the size of an average cluster, since the mean number of OTUs assigned to a cluster was 42.8, with cluster sizes ranging from 3 to 118 OTUs. Since the majority of OTUs are rare, the overall number of OTU reads associated with the outcome is low in this setting. The model for T _{ i } was the same as in setting 3. Setting 4 is comparable to setting 3 when the associated cluster is rare: in both cases, the presence or absence of rare OTUs is associated with survival times. However, it lacks setting 3’s close phylogenetic relationship between associated OTUs.
In all simulation settings, we considered the weighted (K _{w}) and unweighted (K _{u}) UniFrac kernels, the Bray-Curtis kernel (K _{BC}), and the generalized UniFrac kernel with α=0.5 (K _{0.5}). These kernels are expected to have high power in different simulation settings. All of the UniFrac kernels take phylogenetic information into account. The unweighted UniFrac kernel does not account for OTU abundance, whereas the weighted UniFrac kernel does; the generalized UniFrac kernel is intermediate between the weighted and unweighted. The Bray-Curtis kernel does not account for phylogenetic structure or overall abundance of an OTU but does compare both presence/absence and relative abundance between samples of each OTU. Each kernel will have highest power when its measure of distance (and therefore similarity) accurately reflects the true relationship between the microbiome and the outcome.
For each simulation setting, sample size n, and censoring proportion, and using each kernel, we applied the test described above to test for associations between OTU counts and survival time. We used 5000 simulations with γ=0 to estimate the empirical type I error rate with a nominal significance level of 0.05 and estimated empirical power across a range of γ values using 1000 simulations.
We also compare MiRKAT-S to two alternative approaches sometimes used for community-level analysis. First, we performed OTU-level tests of all OTUs. For each of the 856 OTUs in the dataset, we ran a marginal Cox regression model. The minimum p value from the 856 marginal models was compared to the null distribution to produce an overall p value for any association of the microbiota with survival times. In practice, the null distribution would be generated for an individual study using permutation; however, in the interest of computational efficiency, we generated this distribution using the minimum p values from 5000 simulations where survival times were not associated with the microbiota. Second, we performed principal coordinates analysis (PCoA) on a relevant distance matrix (see, e.g., [16]). Since it is not clear how to make PCoA plots with censored time-to-event outcomes, we followed PCoA by Cox proportional hazards regression. Specifically, we generated the UniFrac and Bray-Curtis distance matrices as above, then took the top two principal coordinates as our covariates for a Cox regression analysis. We tested the two principal coordinates jointly by using a chi-square test to compare nested models with and without the microbiota-related predictors. The covariates X _{1} and X _{2} were included in all models exactly as in the MiRKAT-S simulations.
Results
Power and type I error in simulated datasets
Empirical type I errors for n=100, 200, or 500
Number | K _{w} | K _{u} | K _{0.5} | K _{BC} |
---|---|---|---|---|
100 | 0.0544 | 0.0540 | 0.0530 | 0.0542 |
200 | 0.0494 | 0.0480 | 0.0470 | 0.0462 |
500 | 0.0506 | 0.0478 | 0.0536 | 0.0442 |
Empirical type I errors for n<100
Number | Method | K _{w} | K _{u} | K _{0.5} | K _{ BC } |
---|---|---|---|---|---|
25 | MiRKAT-S | 0.054 | 0.055 | 0.055 | 0.056 |
Permutation | 0.046 | 0.048 | 0.046 | 0.049 | |
50 | MiRKAT-S | 0.045 | 0.058 | 0.051 | 0.051 |
Permutation | 0.041 | 0.052 | 0.045 | 0.045 | |
75 | MiRKAT-S | 0.054 | 0.058 | 0.051 | 0.053 |
Permutation | 0.051 | 0.053 | 0.048 | 0.049 |
To interpret the simulation results evaluating the power of the test, recall that two aspects of the relationship between the microbiome and survival are important for understanding which kernel will provide the highest power: the relationship between associated OTUs (whether or not they cluster on a phylogenetic tree) and the importance of taxon abundance (whether OTU count or presence/absence matters). All of the UniFrac distances account for phylogeny, while the Bray-Curtis dissimilarity does not. The weighted UniFrac distance and Bray-Curtis dissimilarity both utilize taxon abundance (OTU counts), whereas the unweighted UniFrac distance only incorporates presence/absence of taxa, and the generalized UniFrac distance is a compromise between the weighted and unweighted UniFrac distances.
The power under settings 2 and 4, in which unclustered OTUs are associated with the outcome, is reported in panels e and f of Fig. 1. When the OTU counts of the ten most common OTUs are associated with survival time (panel e), the Bray-Curtis kernel has highest power, followed by the weighted UniFrac kernel and generalized UniFrac kernel with α=0.5. Since the Bray-Curtis dissimilarity metric does not incorporate phylogenetic information, this distance is designed for unclustered rather than clustered OTUs. However, since it takes abundance into account, the Bray-Curtis kernel performs better when OTU counts are associated with survival (e.g., panel a) rather than OTU presence/absence (e.g., panel d) and when the associated cluster is at least moderately abundant. When the presence or absence of a random 40 OTUs were associated with survival time (panel f), no kernel had high power, but the unweighted UniFrac kernel had higher power than any other tested kernel. The low power is likely due to the rarity of most randomly selected clusters and inability to gain power by utilizing phylogenetic information.
Kernel choice has a strong effect on the power of the test, and different kernels are optimal depending on the nature of the true relationship between the microbiota and survival time. In practice, a kernel representing relationships of particular scientific interest could be selected. For example, if a healthy microbiome at a certain body site has relatively few dominant taxa at high frequencies, and changes in the relative abundance of these taxa is hypothesized to be associated with the time to a disease outcome or death, choosing a kernel that accounts for taxon abundance will have the highest power to detect the hypothesized changes. If there is no specific hypothesized relationship, multiple kernels can be tested and then the resulting p values adjusted for multiple comparisons. Testing the four kernels discussed here is a reasonable starting point, and the limited number of tests reduces the power loss due to adjusting for multiple comparisons. In these simulations, if the analysis is performed using all four kernels included in Fig. 1 (K _{u},K _{0.5},K _{w},K _{BC}) and then the minimum p value after an FDR adjustment is used for testing, the power does not quite match the best kernel but is comparable to or better than the remaining three kernels (see Additional file 1: Figure S2).
We also compared the power of MiRKAT-S to two approaches used in current practice: performing a marginal, or OTU-level, analysis for all OTUs and including the top principal coordinates of the distance matrix as the covariates of interest in a regression model (see Additional file 1: Figure S2). We find that for most simulation settings, MiRKAT-S has substantially better power than the marginal analysis or PCoA, and in the remainder, the power is comparable between methods. In particular, the marginal analysis has power to detect an association between counts of OTUs in a cluster and survival times, but virtually no power to detect associations between presence/absence of clustered OTUs and survival or associations involving unclustered OTUs. The power of the marginal test for detecting an association with clustered OTU counts has similar power regardless of how common the cluster is. Therefore, since MiRKAT-S is more powerful for OTU counts of relatively common clusters, the marginal analysis and MiRKAT-S have similar power for rare clusters, and for very large effect sizes in this simulation setting, the marginal analysis outperforms MiRKAT-S slightly (Additional file 1: Figure S2 (panel C)). However, in all other cases, MiRKAT-S is substantially more powerful than the OTU-level analysis. PCoA regression analysis performs similarly to MiRKAT-S for each kernel when OTU counts in a common cluster are associated with survival times (Additional file 1: Figure S2 (panel A)). In most other simulation settings, PCoA matches or approaches the power of MiRKAT-S in only the case of the best kernel. That is, MiRKAT-S is more robust to kernel choice than PCoA. In addition, in settings where clustering information does not matter (Additional file 1: Figure S2 (panels e and f)), PCoA has very low power, while MiRKAT-S has moderate power provided that the associated OTUs are not too rare.
Analysis of blood and bone marrow transplant data
Acute graft-versus-host disease (aGVHD) is a leading cause of death after allogeneic blood or bone marrow transplantation. There is a suspected relationship between the intestinal microbiome and aGVHD, but previous studies in mice and humans have yielded mixed results about the presence and nature of this relationship. Therefore, Jenq et al. recently studied the association of a particular bacterial species (intestinal Blautia) and of intestinal microbiome diversity indices with time to each of aGVHD onset, aGVHD-related mortality, and adverse outcomes unrelated to aGVHD [5].
In the original study, subjects were stratified into two cohorts depending on sequencing platform. The combined dataset used here results from resequencing of the first cohort of patients using the Illumina MiSeq platform; unfortunately, four patients did not have additional DNA available for MiSeq sequencing and were excluded from our analysis. Therefore, 481 stool samples were available for 111 unique subjects, and for each sample, the Illumina MiSeq platform was used to sequence the V4–V5 regions of the 16S rRNA gene. OTUs were generated as described in [5]. Briefly, mothur version 1.34 was used to compile and process sequence data [36], and quality filters were applied as in [37]. This procedure yielded OTU counts for 2436 OTUs. As in the original paper, for each subject, we only included the sample collected closest to 12 days post-transplant in our analysis, and we excluded subjects for whom no samples were collected between 8 and 16 days post-transplant, so that 94 subjects were included. We used QIIME with default settings to align the sequences and generate a rooted phylogenetic tree. The 109 OTUs that failed to be placed on the tree were excluded from our analysis, leaving 2327 OTUs. We performed the test using the unweighted and weighted UniFrac kernels, the generalized UniFrac kernel with α=0.5, and the Bray-Curtis kernel, adjusting for age and gender. The outcomes considered were overall survival and time to stage III aGVHD.
Analysis of gut microbiome after allogeneic transplant
Outcome | Method | K _{u} | K _{0.5} | K _{w} | K _{BC} |
---|---|---|---|---|---|
Overall Survival | Uncorrected | 0.049 | 0.008 | 0.065 | 0.029 |
Corrected | 0.046 | 0.007 | 0.065 | 0.022 | |
Grade III aGVHD | Uncorrected | 0.496 | 0.514 | 0.472 | 0.849 |
Corrected | 0.560 | 0.575 | 0.518 | 0.933 |
Kaplan-Meier curves for overall survival in the two clusters are shown in Fig. 2 b. However, the simple Cox regression p value is not significant (p=0.09). That is, the similarity between individuals was measured the same way in both analyses. MiRKAT-S yielded a highly significant p value for the association of the microbiota with overall survival, whereas the analysis based on clustering individuals gave a non-significant p value. This demonstrates that MiRKAT-S has higher power to detect this association than a simple clustering analysis based on the same distance metric.
The highly significant result using MiRKAT-S with K _{0.5} may also provide information about the form of the association between the gut microbiota and survival post-transplant. The generalized UniFrac kernel incorporates phylogenetic information and represents a compromise between abundance and presence or absence of OTUs. Therefore, this kernel has highest power to detect relationships between taxonomic profiles and overall survival that occur through moderately rare clusters of OTUs or through a combination of common and rare clusters of OTUs. Accordingly, the high significance of MiRKAT-S using K _{0.5} may indicate that one of those settings holds: either moderately rare clusters of OTUs are driving the relationship between the microbiota and overall survival, or multiple clusters of OTUs, some of which are abundant and some of which are rare, are associated with overall survival. However, without further analysis, we cannot determine which OTUs or clusters are associated with survival times in aGVHD patients.
Discussion
We propose MiRKAT-S for testing the association between the human microbiota and survival outcomes. In the kernel machine Cox model framework, taxonomic profiles are modeled through a kernel function. This allows comparison of microbial community profiles using microbiome-specific distance metrics such as the UniFrac distances or Bray-Curtis dissimilarity. The kernel machine regression framework also allows linear (or, more generally, parametric) adjustment for covariates and potential confounders. We test the significance of the association between the microbiota and survival times using a variance-component score test, and we develop a small-sample correction to account for the modest sample sizes and sparse, high-dimensional data that often result from microbiome studies. In contrast to existing methods that use resampling, p values are obtained analytically using the Davies approximation.
Like other distance-based analyses, MiRKAT-S is limited to detecting the presence of an association between the microbiota and survival times. It cannot identify individual taxa that are associated with the outcome and does not provide information about relationships among taxa within a microbial community. MiRKAT-S is therefore designed to be used when the question of interest is whether an entire microbial community is associated with the outcome. Alternative ways to answer this question include testing the association of each OTU individually with the outcome of interest or using a dimension reduction technique such as PCoA and testing the top few principal coordinates. Our simulation studies show that MiRKAT-S has power at least comparable to, and often substantially greater than, either of these methods for community-level association testing. Community-level tests can be used in combination with other methods that identify taxa of interest. These include marginal tests for particular OTUs of interest, identification of OTUs with high loadings from PCA or PCoA, or penalized regression methods that account for the structure and compositional nature of the data.
Our simulation results show that MiRKAT-S correctly controls type I error. However, under conditions of extreme censoring or very small sample sizes, the analytic p values provided by MiRKAT-S may be slightly anticonservative. In these cases, obtaining p values by permutation may be preferable. Type I error is accurate regardless of the choice of kernel; that is, the test is valid even when a poor choice of kernel is made. The power of the test depends heavily on how well the selected kernel encodes the true relationship between the microbiome and the outcome of interest. For example, when the abundance of an OTU or set of OTUs is related to the outcome, a kernel that encodes abundance information, such as the weighted UniFrac or Bray-Curtis kernel, will have higher power than a kernel that encodes only taxon presence or absence. There are situations in which MiRKAT-S has low power regardless of kernel choice, but any method would have low power in those settings because the true effect size is very small. For example, if the presence or absence of a common OTU is associated with the outcome, nearly all subjects will have the OTU present in the sample, so the association will be very difficult to detect using any method.
If there is no a priori hypothesis about which kernel will best represent relationships of scientific interest, the analysis can be performed using multiple kernels and an overall p value can be obtained by permutation or adjustment for multiple comparisons. This analysis approach can provide information not only about the presence of a relationship but also about its form, depending on the distance metrics considered and their relative power for different forms of the true association. That is, if the metric with the lowest p value has the highest power to detect associations with abundance of common clusters, that may be the form of the true association. Furthermore, weighted combinations of kernels could be used to simultaneously detect different types of shifts in the microbiota. Specific combinations or kernel weights could either be selected a priori or via a grid search, again using permutation to test overall significance. As the field of microbiome analysis matures and new distance metrics are proposed, our approach will continue to increase in power.
Conclusions
We present MiRKAT-S, a method for testing the association between the microbiota, assessed on the community level, and survival (time-to-event) outcomes. Similar methods exist for binary and continuous outcomes; however, MiRKAT-S is the first community-level test for microbiome data that allows analysis of censored survival outcomes. Community-level analyses have several benefits: they often provide higher statistical power to detect associations, and they allow investigators to address additional scientific questions, such as whether the entire microbiome is collectively associated with survival time or time to development of a disease. We use the kernel machine regression framework, encoding microbiome data in ecologically relevant kernels. With judicious choice of kernels, the test can detect a wide range of true forms of association, including association of the outcome with OTU presence/absence or abundance and with either phylogenetically clustered or unclustered sets of taxa. Therefore, MiRKAT-S facilitates a robust community-level analysis of the association between the microbiota and censored survival outcomes that is not possible using existing methods.
Appendix 1
Appendix 2
so that \(\hat {M} = P_{0}^{*} M\), as claimed.
Declarations
Acknowledgements
Not applicable.
Funding
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1256082 (AP). This research was additionally supported by the National Institutes of Health (R01-HL124112, RRJ; U10-CA180819, MCW; R01-HG007508, MCW) and The Hope Foundation (MCW).
Availability of data and materials
The data used in this article are available in the European Bioinformatics Institute Nucleotide Archive under accession number PRJEB15259 (http://www.ebi.ac.uk/ena/data/view/PRJEB15259). The R package MiRKAT-S is available at http://research.fhcrc.org/wu/en/software.htmlor on GitHub at https://github.com/aplantin/MiRKATSin formats appropriate for Macintosh or Windows systems. In addition, each of these sites includes a vignette demonstrating use of the package.
Authors’ contributions
AP developed the method, analyzed the data, and wrote the manuscript. RRJ collected the data and contributed to the data analysis. JC and MCW conceived the study and critically reviewed the manuscript. XZ and NZ contributed to development and testing of the method. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable. Because this study only involved secondary analysis of an existing, de-identified dataset, it is not considered human subject research.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
- Morgan XC, Huttenhower C. Human microbiome analysis. PLoS Comput Biol. 2012; 8(12):1002808.View ArticleGoogle Scholar
- Hamady M, Knight R. Microbial community profiling for human microbiome projects: Tools, techniques, and challenges. Genome Res. 2009; 19(7):1141–1152.View ArticlePubMedPubMed CentralGoogle Scholar
- Wooley JC, Godzik A, Friedberg I. A primer on metagenomics. PLoS Comput Biol. 2010; 6(2):1000667.View ArticleGoogle Scholar
- Stackebrandt E, Goebel B. Taxonomic note: a place for DNA-DNA reassociation and 16s rRNA sequence analysis in the present species definition in bacteriology. Int J Syst Evol Microbiol. 1994; 44(4):846–9.View ArticleGoogle Scholar
- Jenq RR, Taur Y, Devlin SM, Ponce DM, Goldberg JD, Ahr KF, Littmann ER, Ling L, Gobourne AC, Miller LC, et al. Intestinal Blautia is associated with reduced death from graft-versus-host disease. Biogy Blood Marrow Transplants. 2015; 21(8):1373–1383.View ArticleGoogle Scholar
- Fitz-Gibbon S, Tomida S, Chiu BH, Nguyen L, Du C, Liu M, Elashoff D, Erfe MC, Loncaric A, Kim J, et al. Propionibacterium acnes strain populations in the human skin microbiome associated with acne. J Invest Dermatology. 2013; 133(9):2152–160.View ArticleGoogle Scholar
- Charlson ES, Chen J, Custers-Allen R, Bittinger K, Li H, Sinha R, Hwang J, Bushman FD, Collman RG. Disordered microbial communities in the upper respiratory tract of cigarette smokers. PLoS ONE. 2010; 5(12):15216.View ArticleGoogle Scholar
- Cho I, Blaser MJ. The human microbiome: at the interface of health and disease. Nature Reviews Genetics. 2012; 13(4):260–70.PubMedPubMed CentralGoogle Scholar
- Erb-Downward JR, Thompson DL, Han MK, Freeman CM, McCloskey L, Schmidt LA, Young VB, Toews GB, Curtis JL, Sundaram B, et al. Analysis of the lung microbiome in the “healthy” smoker and in COPD. PLoS ONE. 2011; 6(2):16384.View ArticleGoogle Scholar
- Giongo A, Gano KA, Crabb DB, Mukherjee N, Novelo LL, Casella G, Drew JC, Ilonen J, Knip M, Hyöty H, et al. Toward defining the autoimmune microbiome for type 1 diabetes. The ISME J. 2011; 5(1):82–91.View ArticlePubMedGoogle Scholar
- Morris A, Beck JM, Schloss PD, Campbell TB, Crothers K, Curtis JL, Flores SC, Fontenot AP, Ghedin E, Huang L, et al. Comparison of the respiratory microbiome in healthy nonsmokers and smokers. Am J Respiratory Critical Care Med. 2013; 187(10):1067–1075.View ArticleGoogle Scholar
- Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley RE, Sogin ML, Jones WJ, Roe BA, Affourtit JP, et al. A core gut microbiome in obese and lean twins. Nature. 2009; 457(7228):480–4.View ArticlePubMedGoogle Scholar
- Karlsson FH, Tremaroli V, Nookaew I, Bergström G, Behre CJ, Fagerberg B, Nielsen J, Bäckhed F. Gut metagenome in European women with normal, impaired and diabetic glucose control. Nature. 2013; 498(7452):99–103.View ArticlePubMedGoogle Scholar
- Zhao N, Chen J, Carroll IM, Ringel-Kulka T, Epstein MP, Zhou H, Zhou JJ, Ringel Y, Li H, Wu MC. Testing in microbiome-profiling studies with MiRKAT, the microbiome regression-based kernel association test. Am J Human Genet. 2015; 96(5):797–807.View ArticleGoogle Scholar
- Chen J, Li H. Kernel Methods for Regression Analysis of Microbiome Compositional Data In: Hu M, Liu Y, Lin J, editors. Topics in Applied Statistics, vol 55. New York: Springer: 2013. p. 191–201.Google Scholar
- Han MK, Zhou Y, Murray S, Tayob N, Noth I, Lama VN, Moore BB, White ES, Flaherty KR, Huffnagle GB, et al. Lung microbiome and disease progression in idiopathic pulmonary fibrosis: an analysis of the COMET study. Lancet Respir Med. 2014; 2(7):548–56.View ArticlePubMedPubMed CentralGoogle Scholar
- Penders J, Gerhold K, Stobberingh EE, Thijs C, Zimmermann K, Lau S, Hamelmann E. Establishment of the intestinal microbiota and its role for atopic dermatitis in early childhood. J Allergy Clin Immunol. 2013; 132(3):601–7.View ArticlePubMedGoogle Scholar
- Bisgaard H, Li N, Bonnelykke K, Chawes BLK, Skov T, Paludan-Müller G, Stokholm J, Smith B, Krogfelt KA. Reduced diversity of the intestinal microbiota during infancy is associated with increased risk of allergic disease at school age. J Allergy Clin Immunol. 2011; 128(3):646–52.View ArticlePubMedGoogle Scholar
- Lin X, Cai T, Wu MC, Zhou Q, Liu G, Christiani DC, Lin X. Kernel machine SNP-set analysis for censored survival outcomes in genome-wide association studies. Genet Epidemiol. 2011; 35(7):620–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Cai T, Tonini G, Lin X. Kernel machine approach to testing the significance of multiple genetic markers for risk prediction. Biometrics. 2011; 67(3):975–86.View ArticlePubMedPubMed CentralGoogle Scholar
- Cox DR. Regression models and life tables (with discussion). J Royal Stat Soc. 1972; 34:187–220.Google Scholar
- Kimeldorf GS, Wahba G. A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. Annals Math Stat. 1970; 41(2):495–502.View ArticleGoogle Scholar
- Lozupone C, Knight R. UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol. 2005; 71(12):8228–235.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen J, Bittinger K, Charlson ES, Hoffmann C, Lewis J, Wu GD, Collman RG, Bushman FD, Li H. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics. 2012; 28(16):2106–113.View ArticlePubMedPubMed CentralGoogle Scholar
- Bray JR, Curtis JT. An ordination of the upland forest communities of southern Wisconsin. Ecol Monogr. 1957; 27(4):325–49.View ArticleGoogle Scholar
- Liu D, Lin X, Ghosh D. Semiparametric regression of multidimensional genetic pathway data: least-squares kernel machines and linear mixed models. Biometrics. 2007; 63(4):1079–1088.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen H, Lumley T, Brody J, Heard-Costa NL, Fox CS, Cupples LA, Dupuis J. Sequence kernel association test for survival traits. Genet Epidemiol. 2014; 38(3):191–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Styan GP. Notes on the distribution of quadratic forms in singular normal variables. Biometrika. 1970; 57(3):567–72.View ArticleGoogle Scholar
- Efron B. The efficiency of Cox’s likelihood function for censored data. J Am Stat Assoc. 1977; 72(359):557–65.View ArticleGoogle Scholar
- Hertz-Picciotto I, Rockhill B. Validity and efficiency of approximation methods for tied survival times in Cox regression. Biometrics. 1997; 53(3):1151–1156.View ArticlePubMedGoogle Scholar
- Chen J, Chen W, Zhao N, Wu MC, Schaid DJ. Small sample kernel association tests for human genetic and microbiome association studies. Genet Epidemiol. 2016; 40(1):5–19.View ArticlePubMedGoogle Scholar
- Nygård S, Borgan Ø, Lingjærde OC, Størvold HL. Partial least squares Cox regression for genome-wide data. Lifetime Data Anal. 2008; 14(2):179–95.View ArticlePubMedGoogle Scholar
- Park PJ, Tian L, Kohane IS. Linking gene expression data with patient survival times using partial least squares. Bioinformatics. 2002; 18(suppl 1):120–7.View ArticleGoogle Scholar
- Lambert-Lacroix S, Letué F, et al. Partial least squares and Cox model with application to gene expression. Technical report. 2011. http://sites.uclouvain.be/IAP-Stat-Phase-V-VI/PhaseVI/publications_2011/TR/Lambert-LacroixLetueIAP.pdf.
- Davies RB. The distribution of a linear combination of chi-2 random variables. J R Stat Soc Ser C (Appl Stat). 1980; 29(3):323–33.Google Scholar
- Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, Lesniewski RA, Oakley BB, Parks DH, Robinson CJ, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009; 75(23):7537–541.View ArticlePubMedPubMed CentralGoogle Scholar
- Schloss PD, Gevers D, Westcott SL. Reducing the effects of pcr amplification and sequencing artifacts on 16s rRNA-based studies. PloS ONE. 2011; 6(12):27310.View ArticleGoogle Scholar
- Ward Jr JH. Hierarchical grouping to optimize an objective function. J Am Stat Assoc. 1963; 58(301):236–44.View ArticleGoogle Scholar
- Murtagh F, Legendre P. Ward’s hierarchical agglomerative clustering method: Which algorithms implement Ward’s criterion?J Classif. 2014; 31(3):274–95.View ArticleGoogle Scholar