Example microbiome dataTo ease the illustration of MiMedSurv, we employed example microbiome data to survey the mediating roles of the gut microbiome between antibiotic treatment and T1D onset11. T1D is an autoimmune disease, in which the human immune system is overly active and hence attacks even non-pathogenic normal cells. Unfortunately, T1D is increasing in incidence while decreasing in its onset age world-wide. Zhang et al.11 revealed that early-life antibiotic (known as tylosin) administration can accelerate T1D development through the gut microbiome dysbiosis. For illustration purposes, we reanalyze the data using (i) the antibiotic (tylosin) treatment as a treatment variable, (ii) the gut microbiome as a mediator, (iii) the time-to-T1D as a survival variable, and (iv) sex and elapse time after the antibiotic (tylosin) treatment as two covariate variables (Fig. 2). We expected sex and the elapse time after the antibiotic (tylosin) treatment to be potential confounders between the antibiotic (tylosin) treatment and the onset of T1D. Especially, as the elapse time increases, the effect of the antibiotic on gut microbiome can decrease and the elapse time can also be linked to the onset of T1D. Here, the time-to-T1D survival variable consists of typical right-censored survival responses with two components of (i) follow-up time (in week) and (ii) censored/event indicator with 0 for censored (no event) T1D free and 1 for event that is T1D onset.Fig. 2A graphical representation on our example microbiome causal mediation analysis with survival responses. (i) antibiotic (tylosin) treatment is a treatment variable; (ii) gut microbiome is a mediator; (iii) time-to-T1D is a survival variable; and (iv) sex and elapse time after the antibiotic (tylosin) treatment are covariate variables.The 16S ribosomal RNA gene amplicon sequence data are publicly available in the European Nucleotide Archive (ENA) (https://www.ebi.ac.uk/ena/browser/view/PRJEB14696, accession number: PRJEB14696). We processed them using a bioinformatic pipeline, QIIME 2 (https://qiime2.org)47, and a database, Greengenes (https://greengenes.secondgenome.com), to construct a feature table, taxonomic annotations, and a phylogenetic tree. We stored the final processed data as example data in the Data Input module of MiMedSurv for users to easily confirm data formatting requirements.In the following sections, we describe all the preprocessing and analytic modules of MiMedSurv step-by-step using this example microbiome data (see Application Note).Data processing: data inputUsers first need to upload their microbiome data with (i) four components: a feature (i.e., operational taxonomic unit (OTU) or amplicon sequence variants (ASV)) table, a taxonomic table, a phylogenetic tree, and metadata or (ii) three components: a feature (i.e., OTU or ASV) table, a taxonomic table, and metadata. Here, the taxonomic table must contain seven taxonomic ranks, kingdom, phylum, class, order, family and species, while the phylogenetic tree much be a rooted tree that reflects evolutionary relationships across features (i.e., OTUs or ASVs). If users upload their microbiome dataset with three components without a phylogenetic tree, only non-phylogenetic community-level (alpha- and beta-diversity) analyses will later be performed.Users can start with downloading the example microbiome data from the Example Data section. The data are stored in a widely used unified format, called phyloseq48, that can efficiently combine all essential microbiome data components. Alternatively, users can also employ them all individually. In this module, we also described all the instructions with example codes to check compatible data formats so that users can prepare microbiome data easily.Application Note: We uploaded the example microbiome data using the Browse and Upload buttons.Data processing: quality controlAs in MiMed20, users can perform quality controls with respect to (i) a microbial kingdom of interest (default: Bacteria), (ii) a minimum library size (i.e., total read count) for the study subjects to be retained (default: 3000), (iii) a minimum mean proportion for the features (i.e., OTUs or ASVs) to be retained (default: 0.02%), and (iv) errors in taxonomic names to be removed. For reference, users can download the resulting microbiome data after quality controls.Application note: We performed quality controls using the default settings and simply clicking the Run button. Then, we rescued 519 study subjects for 224 features, 7 phyla, 12 classes, 15 orders, 19 families, 25 genera, and 9 species (Fig. 3). We can see that the library sizes (i.e., total read counts) vary dramatically by study subjects, and the mean proportions are highly skewed to the left (Fig. 3).Fig. 3The status of the resulting microbiome data after quality controls: (i) sample size, the number of features, the number of phyla, the number of classes, the number of orders, the number of families, the number of genera, and the number of species after quality controls; (ii) histogram and box plot for the distribution of library sizes (i.e., total read counts) for subjects; (iii) histogram and box plot for the distribution of mean proportions for features.Data processing: data transformationAs in MiMed20, for the community-level analysis, users can compute nine alpha diversity indices (i.e., Observed, Shannon49, Simpson50, Inverse Simpson50, Fisher51, Chao152, abundance-based coverage estimator (ACE)53, incidence-based coverage estimator (ICE)54, phylogenetic diversity (PD)55) and five beta diversity indices (i.e., Jaccard dissimilarity56, Bray–Curtis dissimilarity57, Unweighted UniFrac distance58, Generalized UniFrac distance59, Weighted UniFrac distance60). For the taxonomy-level analysis, users can normalize the data using the widely used centered log-ratio (CLR) transformation method61 to relax the compositional constraint of the data, yet three other normalization methods of arcsine-root, rarefied count62 and proportion are also available. For reference, users can download all the resulting alpha and beta diversity indices and normalized taxonomic data.Application Note: We computed all the alpha and beta diversity indices and normalized taxonomic data simply clicking the Run button.Data analysis: non-mediational analysisIn this module, users can perform some baseline exploratory survival analysis to survey the disparity in survival time between medical treatments (e.g., treatment vs. placebo, new treatment vs. old treatment, and so forth) or environmental exposures (e.g., rural vs. urban, smoking vs. non-smoking, calorie restriction vs. ad libitum diet, and so forth). For this, users need to select (i) a survival time variable, (ii) a censored/event indicator variable, (iii) a treatment variable, (iv) covariate(s), and (v) an analytic method as we described in the previous section, Materials and Methods: Non-Mediational Analysis.Application note: We selected (i) ‘T1Dweek’ as a survival time variable, (ii) ‘T1D’ as a censored/event indicator variable, (iii) ‘Antibiotic’ as a treatment variable, (iv) ‘Sex’ and ‘SampleTime’ as covariates, and (v) the Cox model34 as an analytic method. Then, we found a significant disparity in survival (against T1D) rate between the normal control and antibiotic (tylosin) groups adjusting for sex and elapse time after the antibiotic (tylosin) treatment (P-value: < 0.001) (Fig. 4). We also estimated that the hazard rate (toward T1D) is higher for the antibiotic (tylosin) group than the normal control group (hazard ratio (HR): 1.758 > 1), indicating that the antibiotic (tylosin) treatment is harmful accelerating T1D onset (Fig. 4).Fig. 4The results from non-mediational analysis on the disparity in survival (against T1D) rate between the normal control and antibiotic (tylosin) groups adjusting for sex and elapse time after the antibiotic (tylosin) treatment. HR represents the hazard ratio, that is the hazard rate to T1D for the antibiotic (tylosin) group divided by the hazard rate to T1D for the normal control group.Data analysis: mediational analysisIn the previous section, Results: Data Analysis: Non-Mediational Analysis, we observed a significant disparity in survival rate between the normal control and antibiotic (tylosin) groups adjusting for sex and elapse time after the antibiotic (tylosin) treatment. In the following sections, we seek to comprehend if the gut microbiome plays a mediating role between the antibiotic (tylosin) treatment and T1D onset with respect to (i) alpha diversity, (ii) beta diversity and (iii) microbial taxa at different hierarchies (i.e., phyla, classes, orders, families, and genera).Community-level analysis: alpha diversityIn this module, users can perform the microbiome causal mediation analysis to test jointly if the treatment alters the microbial alpha diversity, and then the altered microbial alpha diversity, in turn, influences the survival responses. For this, users need to select (i) a survival time variable, (ii) a censored/event indicator variable, (iii) a treatment variable, and (iv) covariate(s) for covariate-adjusted analysis or not for univariate analysis. The available analytic method here is the Imai method35,36 coupled with the Weibull regression model37 as we described in the previous section, Materials and Methods: Mediational Analysis.Application note: We selected (i) ‘T1Dweek’ as a survival time variable, (ii) ‘T1D’ as a censored/event indicator variable, (iii) ‘Antibiotic’ as a treatment variable, (iv) ‘Sex’ and ‘SampleTime’ as covariates, and (v) the Imai method35,36 as an analytic method. Then, we found a significant mediation effect of microbial alpha diversity with respect to the Simpson index50 between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment (P-value: 0.045 < 0.05) [Fig. 5], yet all the other alpha diversity indices are not statistically significant. We also estimated that as the Simpson index50 level decreases, the T1D onset tends to be accelerated (Est. − 1.720 < 0), indicating a lower Simpson index level is harmful to T1D [Fig. 5].Fig. 5The results from mediational analysis between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment using alpha diversity indices (i.e., Observed, Shannon49, Simpson50, Inverse Simpson50, Fisher51, Chao152, ACE53, ICE54, PD55). ACME represents average causal mediation effect35,36.Community-level analysis: beta diversityIn this module, users can perform the microbiome causal mediation analysis to test jointly if the treatment alters the microbial beta diversity, and then the altered microbial beta diversity, in turn, influences the survival responses. For this, users need to select (i) a survival time variable, (ii) a censored/event indicator variable, (iii) a treatment variable, and (iv) covariate(s) for covariate-adjusted analysis or not for univariate analysis. The available analytic method here is DACT41 coupled with MiRKAT38,39 for the treatment model and MiRKAT-S39,40 for the outcome model as we described in the previous section, Materials and Methods: Mediational Analysis.Application Note: We selected (i) ‘T1Dweek’ as a survival time variable, (ii) ‘T1D’ as a censored/event indicator variable, (iii) ‘Antibiotic’ as a treatment variable, (iv) ‘Sex’ and ‘SampleTime’ as covariates, and (v) DACT41 as an analytic method. Then, we could not find any significant mediation effect of microbial beta diversity between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment [Fig. 6].Fig. 6The results from mediational analysis between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment using beta diversity indices (i.e., Jaccard dissimilarity56, Bray–Curtis dissimilarity57, Unweighted UniFrac distance58, Generalized UniFrac distance59, Weighted UniFrac distance60). Censored represents the subjects who have not had T1D during the follow-up, while Event represents the subjects who have had T1D during the follow-up.Taxonomy-level analysisIn this module, users can perform the microbiome causal mediation analysis to test jointly if the treatment alters microbial taxa, and then the altered microbial taxa, in turn, influence the survival responses. For this, users need to select (i) a data format (default: CLR61), (ii) a survival time variable, (iii) a censored/event indicator variable, (iv) a treatment variable, and (v) covariate(s) for covariate-adjusted analysis or not for univariate analysis, and (vi) the taxonomic ranks to be analyzed from phylum to genus for 16S ribosomal RNA gene amplicon sequencing1,2 or from phylum to species for shotgun metagenomics3. The available analytic method here is the Imai method35,36 coupled with the Weibull regression model37 as we described in the previous section, Materials and Methods: Mediational Analysis.Application note: We selected (i) ‘CLR’61 as a data format, (ii) ‘T1Dweek’ as a survival time variable, (iii) ‘T1D’ as a censored/event indicator variable, (iv) ‘Antibiotic’ as a treatment variable, (v) ‘Sex’ and ‘SampleTime’ as covariates, (vi) ‘Phylum-Genus’ as taxonomic ranks to be analyzed, and (vii) the Imai method35,36 as an analytic method. Then, we found a significant mediation effect of the two phyla (Antinobacteria and Verrucomicrobia), three classes (Erysipelotrichi, Actinobacteria, and Verrucomicrobiae), three orders (Erysipelotrichales, Bifidobacteriales, and Verrucomicrobiales), three families (Erysipelotrichaceae, Bifidobacteriaceae, and Verrucomicrobiaceae), and three genera (Allobaculum, Bifidobacterium, and Akkermansia) between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment (Q-value: < 0.001) (Fig. 7). We also estimated that as their relative abundance level decreases, the T1D onset tends to be accelerated (Est. < 0), indicating that they might be beneficial microbes to prevent T1D onset (Fig. 7).Fig. 7The results from mediational analysis between the antibiotic (tylosin) treatment and T1D onset adjusting for sex and elapse time after the antibiotic (tylosin) treatment using microbial taxa from phylum to genus. Q-value represents FDR-adjusted P-value43.
