MMUPHin 1.0.0
adjust_batch
lm_meta
discrete_discover
continuous_discover
MMUPHin
is an R package implementing meta-analysis methods for microbial
community profiles. It has interfaces for: a) covariate-controlled batch and
study effect adjustment, b) meta-analytic differential abundance testing, and meta-analytic discovery of c) discrete (cluster-based) or d) continuous
unsupervised population structure.
Overall, MMUPHin
enables the normalization and combination of multiple
microbial community studies. It can then help in identifying microbes, genes,
or pathways that are differential with respect to combined phenotypes. Finally,
it can find clusters or gradients of sample types that reproduce consistently
among studies.
This vignette is intended to provide working examples for all four
functionalities of MMUPHin
.
library(MMUPHin)
# tidyverse packages for utilities
library(magrittr)
library(dplyr)
library(ggplot2)
MMUPHin is a Bioconductor package and can be installed via the following command.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("MMUPHin")
As input, MMUPHin
requires a properly
formatted collection of microbial community studies, with both feature
abundances and accompanying metadata. Here we use the five published
colorectal cancer (CRC)
stool metagenomic studies, incorporated in Thomas et al. (2019).
Data for the studies are already
conveniently packaged and accessible through the Biocondcutor package
curatedMetagenomicData,
though additional wranglings are needed to format input for MMUPHin
.
Importantly, MMUPHin
asks that feature abundances be provided as a
feature-by-sample matrix, and the metadata be provided as a data frame.
The two objects shoud agree on sample IDs, that is, rowname
of the feature
abundance matrix and colname
of the metadata data frame must agree.
Many popular ’omic data classes in R already enforce this standard, such as
ExpressionSet
from Biobase, or phyloseq
from
phyloseq.
To minimize users’ efforts in loading data to run the examples, we have properly formatted the five studies for easy access. The feature abundances and metadata can be loaded with the following code chunk. For the interested user, the commented out scripts were used for accessing data directly from curatedMetagenomicData and formatting. It might be worthwhile to read through these as they perform many of the common tasks for preprocessing microbial feature abundance data in R, including sample/feature subsetting, normalization, filtering, etc.
data("CRC_abd", "CRC_meta")
# CRC_abd is the feature (species) abundance matrix. Rows are features and
# columns are samples.
CRC_abd[1:5, 1, drop = FALSE]
## FengQ_2015.metaphlan_bugs_list.stool:SID31004
## s__Faecalibacterium_prausnitzii 0.11110668
## s__Streptococcus_salivarius 0.09660736
## s__Ruminococcus_sp_5_1_39BFAA 0.09115385
## s__Subdoligranulum_unclassified 0.05806767
## s__Bacteroides_stercoris 0.05685503
# CRC_meta is the metadata data frame. Columns are samples.
CRC_meta[1, 1:5]
## subjectID body_site
## FengQ_2015.metaphlan_bugs_list.stool:SID31004 SID31004 stool
## study_condition
## FengQ_2015.metaphlan_bugs_list.stool:SID31004 CRC
## disease
## FengQ_2015.metaphlan_bugs_list.stool:SID31004 CRC;fatty_liver;hypertension
## age_category
## FengQ_2015.metaphlan_bugs_list.stool:SID31004 adult
# A total of five studies are included
table(CRC_meta$studyID)
##
## FengQ_2015.metaphlan_bugs_list.stool
## 107
## HanniganGD_2017.metaphlan_bugs_list.stool
## 55
## VogtmannE_2016.metaphlan_bugs_list.stool
## 104
## YuJ_2015.metaphlan_bugs_list.stool
## 128
## ZellerG_2014.metaphlan_bugs_list.stool
## 157
# The following were used to access and format the two objects
# library(curatedMetagenomicData)
# library(phyloseq)
# library(genefilter)
# datasets <- curatedMetagenomicData(
# c("FengQ_2015.metaphlan_bugs_list.stool" ,
# "HanniganGD_2017.metaphlan_bugs_list.stool",
# "VogtmannE_2016.metaphlan_bugs_list.stool",
# "YuJ_2015.metaphlan_bugs_list.stool",
# "ZellerG_2014.metaphlan_bugs_list.stool"),
# dryrun = FALSE)
# # Construct phyloseq object from the five datasets
# physeq <-
# # Aggregate the five studies into ExpressionSet
# mergeData(datasets) %>%
# # Convert to phyloseq object
# ExpressionSet2phyloseq() %>%
# # Subset samples to only CRC and controls
# subset_samples(study_condition %in% c("CRC", "control")) %>%
# # Subset features to species
# subset_taxa(!is.na(Species) & is.na(Strain)) %>%
# # Normalize abundances to relative abundance scale
# transform_sample_counts(function(x) x / sum(x)) %>%
# # Filter features to be of at least 1e-5 relative abundance in five samples
# filter_taxa(kOverA(5, 1e-5), prune = TRUE)
# CRC_abd <- otu_table(physeq)@.Data
# CRC_meta <- data.frame(sample_data(physeq))
# CRC_meta$studyID <- factor(CRC_meta$studyID)
adjust_batch
adjust_batch
aims to correct for technical study/batch effects in microbial
feature abundances. It takes as input the feature-by-sample abundance matrix,
and performs batch effect adjustment given provided batch and optional covariate variables. It returns the batch-adjusted abundance matrix. Check
help(adjust_batch)
for additional details and options.
Here we use adjust_batch
to correct for differences in the five studies,
while controlling for the effect of CRC versus control (variable
study_condition
in CRC_meta
).
# The function call indicates for adjust_batch to correct for the effect
# of the batch variable, studyID, while controlling for the effect of the
# disease variable, study_condition. Many additional options are available
# through the control parameter, here we specify verbose=FALSE to avoid
# excessive messages, although they can often be helpful in practice!
fit_adjust_batch <- adjust_batch(feature_abd = CRC_abd,
batch = "studyID",
covariates = "study_condition",
data = CRC_meta,
control = list(verbose = FALSE))
# Note that adjust_batch returns a list of more than one components, and
# feature_abd_adj is the corrected feature abundance matrix. See
# help(adjust_batch) for the meaning of other components.
CRC_abd_adj <- fit_adjust_batch$feature_abd_adj
One way to evaluate the effect of batch adjustment is to assess the total
variability in microbial profiles attributable to study differences, before
and after adjustment. This is commonly known as a PERMANOVA test
(Tang, Chen, and Alekseyenko 2016), and can be performed with the adonis
function in
vegan.
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
# adonis requires as input sample-versus-sample dissimilarities between
# microbial profiles
D_before <- vegdist(t(CRC_abd))
D_after <- vegdist(t(CRC_abd_adj))
# fix random seed as adonis runs randomized permutations
set.seed(1)
fit_adonis_before <- adonis(D_before ~ studyID, data = CRC_meta)
fit_adonis_after <- adonis(D_after ~ studyID, data = CRC_meta)
print(fit_adonis_before)
##
## Call:
## adonis(formula = D_before ~ studyID, data = CRC_meta)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## studyID 4 13.36 3.3400 11.855 0.07991 0.001 ***
## Residuals 546 153.82 0.2817 0.92009
## Total 550 167.18 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(fit_adonis_after)
##
## Call:
## adonis(formula = D_after ~ studyID, data = CRC_meta)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## studyID 4 10.315 2.57863 9.0054 0.06189 0.001 ***
## Residuals 546 156.344 0.28634 0.93811
## Total 550 166.658 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We can see that, before study effect adjustment, study differences can expalin a total of 7.99% of the variability in microbial abundance profiles, whereas after adjustment this was reduced to 6.19%, though the effect of study is significant in both cases.
lm_meta
One of the most common meta-analysis goals is to combine association effects
across batches/studies to identify consistent overall effects. lm_meta
provides a straightforward interface to this task, by first performing
regression analysis in individual batches/studies using the well-validated
Maaslin2 packge, and then aggregating results with
established fixed/mixed effect models, realized via the
vegan package. Here, we use lm_meta
to test for
consistently differential abundant species between CRC and control samples
across the five studies, while controlling for demographic covariates (gender,
age, BMI).
# lm_meta runs regression and meta-analysis models to identify consistent
# effects of the exposure (study_condition, i.e., disease) on feature_abd
# (microbial feature abundances). Batch variable (studyID) needs to be
# specified to identify different studies. Additional covariates to include in
# the regression model can be specified via covariates (here set to gender,
# age, BMI). Check help(lm_meta) for additional parameter options.
# Note the warnings: lm_meta can tell if a covariate cannot be meaningfully fit
# within a batch and will inform the user of such cases through warnings.
fit_lm_meta <- lm_meta(feature_abd = CRC_abd,
batch = "studyID",
exposure = "study_condition",
covariates = c("gender", "age", "BMI"),
data = CRC_meta,
control = list(verbose = FALSE))
## Warning in lm_meta(feature_abd = CRC_abd, batch = "studyID", exposure = "study_condition", : Covariate gender is missing or has only one non-missing value in the following batches; will be excluded from model for these batches:
## HanniganGD_2017.metaphlan_bugs_list.stool, YuJ_2015.metaphlan_bugs_list.stool
## Warning in lm_meta(feature_abd = CRC_abd, batch = "studyID", exposure = "study_condition", : Covariate age is missing or has only one non-missing value in the following batches; will be excluded from model for these batches:
## HanniganGD_2017.metaphlan_bugs_list.stool, YuJ_2015.metaphlan_bugs_list.stool
## Warning in lm_meta(feature_abd = CRC_abd, batch = "studyID", exposure = "study_condition", : Covariate BMI is missing or has only one non-missing value in the following batches; will be excluded from model for these batches:
## HanniganGD_2017.metaphlan_bugs_list.stool, YuJ_2015.metaphlan_bugs_list.stool
# Again, lm_meta returns a list of more than one components.
# meta_fits provides the final meta-analytical testing results. See
# help(lm_meta) for the meaning of other components.
meta_fits <- fit_lm_meta$meta_fits
We can visualize the significant (FDR q < 0.05) species associated with CRC in these studies/samples. Comparing them with Figure 1b of Thomas et al. (2019), we can see that many of the significant species do agree.
meta_fits %>%
filter(qval.fdr < 0.05) %>%
arrange(coef) %>%
mutate(feature = factor(feature, levels = feature)) %>%
ggplot(aes(y = coef, x = feature)) +
geom_bar(stat = "identity") +
coord_flip()
discrete_discover
Clustering analysis of microbial profiles can help identify meaningful
discrete population subgroups (Ravel et al. 2011), but must be carried out
carefully with validations to ensure that the identified structures are
consistent (Koren et al. 2013). discrete_discover
provides the functionality to
use prediction strength (Tibshirani and Walther 2005) to evaluate discrete
clustering structures within individual microbial studies, as well as a
“generalized predicition strength” to evaluate their reproducibility in
other studies. These jointly provide meta-analytical evidence for
(or against) identifying discrete population structures in microbial profiles.
Check help(discrete_discover)
to see more details on the method and additional
options.
The gut microbiome is known to form gradients rather than discrete clusters
(Koren et al. 2013). Here we use discrete_discover
to evaluate clustering
structures among control samples in the five stool studies.
# First subset both feature abundance table and metadata to only control samples
control_meta <- subset(CRC_meta, study_condition == "control")
control_abd_adj <- CRC_abd_adj[, rownames(control_meta)]
# discrete_discover takes as input sample-by-sample dissimilarity measurements
# rather than abundance table. The former can be easily computed from the
# latter with existing R packages.
D_control <- vegdist(t(control_abd_adj))
fit_discrete <- discrete_discover(D = D_control,
batch = "studyID",
data = control_meta,
control = list(k_max = 8,
verbose = FALSE))
The internal_mean
and internal_sd
components of fit_discrete
are matrices
that provide internal evaluation statistics (prediction strength) for each
batch (column) and evaluated number of clusters (row). Similarly,
external_mean
and external_sd
provide external evaluation statistics (
generalized prediction strenght). Evidence for existence of discrete structures
would be a “peaking” of the mean statistics at a particular cluster number.
Here, for easier examination of such a pattern, we visualize the results
for the largest study, ZellerG_2014. Note that visualization for all studies
are by default automatically generated and saved to the output file
“diagnostic_discrete.pdf”.
internal <- data.frame(
# By default, fit_discrete evaluates cluster numbers 2-10
K = 2:8,
statistic =
fit_discrete$internal_mean[, "ZellerG_2014.metaphlan_bugs_list.stool"],
se =
fit_discrete$internal_se[, "ZellerG_2014.metaphlan_bugs_list.stool"],
type = "internal")
external <- data.frame(
# By default, fit_discrete evaluates cluster numbers 2-10
K = 2:8,
statistic =
fit_discrete$external_mean[, "ZellerG_2014.metaphlan_bugs_list.stool"],
se =
fit_discrete$external_se[, "ZellerG_2014.metaphlan_bugs_list.stool"],
type = "external")
rbind(internal, external) %>%
ggplot(aes(x = K, y = statistic, color = type)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_line(position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin = statistic - se, ymax = statistic + se),
position = position_dodge(width = 0.5), width = 0.5) +
ggtitle("Evaluation of discrete structure in control stool microbiome (ZellerG_2014)")
The decreasing trend for both the internal and external statistics along with number of clusters (K) suggests that discrete structures cannot be well-established. To provide a positive example, we examine the two vaginal microbiome studies provided by curatedMetagenomicData, as the vaginal microbiome is known to have distinct subtypes (Ravel et al. 2011). Again, we pre-formatted these datasets for easy access, but you can recreate them from curatedMetagenomicData with the commented out scripts.
# library(curatedMetagenomicData)
# library(phyloseq)
# datasets <- curatedMetagenomicData(
# "*metaphlan_bugs_list.vagina*",
# dryrun = FALSE)
# # Construct phyloseq object from the five datasets
# physeq <-
# # Aggregate the five studies into ExpressionSet
# mergeData(datasets) %>%
# # Convert to phyloseq object
# ExpressionSet2phyloseq() %>%
# # Subset features to species
# subset_taxa(!is.na(Species) & is.na(Strain)) %>%
# # Normalize abundances to relative abundance scale
# transform_sample_counts(function(x) x / sum(x)) %>%
# # Filter features to be of at least 1e-5 relative abundance in two samples
# filter_taxa(kOverA(2, 1e-5), prune = TRUE)
# vaginal_abd <- otu_table(physeq)@.Data
# vaginal_meta <- data.frame(sample_data(physeq))
# vaginal_meta$studyID <- factor(vaginal_meta$studyID)
data("vaginal_abd", "vaginal_meta")
D_vaginal <- vegdist(t(vaginal_abd))
fit_discrete_vag <- discrete_discover(D = D_vaginal,
batch = "studyID",
data = vaginal_meta,
control = list(verbose = FALSE,
k_max = 8))
## Warning: Removed 14 rows containing missing values (geom_errorbar).
# Examine results for the larger study, HMP_2012
data.frame(
# By default, fit_discrete evaluates cluster numbers 2-10
K = 2:8,
statistic =
fit_discrete_vag$internal_mean[, "HMP_2012.metaphlan_bugs_list.vagina"],
se =
fit_discrete_vag$internal_se[, "HMP_2012.metaphlan_bugs_list.vagina"]) %>%
ggplot(aes(x = K, y = statistic)) +
geom_point() + geom_line() +
geom_errorbar(aes(ymin = statistic - se, ymax = statistic + se),
width = 0.5) +
ggtitle("Evaluation of discrete structure in vaginal microbiome (HMP_2012)")
We can see that for the vaginal microbiome, discrete_discover
suggests the
existence of five clusters. Here we examine only the internal metrics of
HMP_2012 as the other study (FerrettiP_2018) has only
19
samples.
continuous_discover
Population structure in the microbiome can manifest as gradients rather than
discrete clusters, such as dominant phyla trade-off or disease-associated
dysbiosis. continuous_discover
provide functionality to identify such
structures as well as to validate them with meta-analysis. We again evaluate
these continuous structures in control samples of the five studies.
# Much like adjust_batch and lm_meta, continuous_discover also takes
# as input feature-by-sample abundances. control offers many tuning parameters
# and here we set one of them, var_perc_cutoff, to 0.5, which asks the method
# to include top principal components within each batch that in total explain
# at least 50% of the total variability in the batch. See
# help(continuosu_discover) for more details on the tuning parameters and
# their interpretations.
fit_continuous <- continuous_discover(feature_abd = control_abd_adj,
batch = "studyID",
data = control_meta,
control = list(var_perc_cutoff = 0.5,
verbose = FALSE))
We can visualize the identified continuous structure scores in at least two ways: first, to examine their top contributing microbial features, to get an idea of what the score is characterizing, and second, to overlay the continuous scores with an ordination visualization. Here we perform these visualizations on the first identified continuous score.
# Examine top loadings
loading <- data.frame(feature = rownames(fit_continuous$consensus_loadings),
loading1 = fit_continuous$consensus_loadings[, 1])
loading %>%
dplyr::arrange(-abs(loading1)) %>%
dplyr::slice(1:20) %>%
dplyr::arrange(loading1) %>%
dplyr::mutate(feature = factor(feature, levels = feature)) %>%
ggplot(aes(x = feature, y = loading1)) +
geom_bar(stat = "identity") +
coord_flip() +
ggtitle("Features with top loadings")
# Ordinate the samples
mds <- cmdscale(d = D_control)
colnames(mds) <- c("Axis1", "Axis2")
as.data.frame(mds) %>%
dplyr::mutate(score1 = fit_continuous$consensus_scores[, 1]) %>%
ggplot(aes(x = Axis1, y = Axis2, color = score1)) +
geom_point() +
coord_fixed()
From ordination we see that the first continuos score indeed represent strong variation across these stool samples. From the top loading features we can see that this score strongly represents a Bacteroidetes (the Bacteroides species) versus Firmicutes (the Ruminococcus species) tradeoff.
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] vegan_2.5-6 lattice_0.20-38 permute_0.9-5 ggplot2_3.2.1
## [5] dplyr_0.8.3 magrittr_1.5 MMUPHin_1.0.0 BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] mclust_5.4.5 Rcpp_1.0.2 mvtnorm_1.0-11
## [4] tidyr_1.0.0 class_7.3-15 assertthat_0.2.1
## [7] zeallot_0.1.0 digest_0.6.22 plyr_1.8.4
## [10] R6_2.4.0 backports_1.1.5 stats4_3.6.1
## [13] pcaPP_1.9-73 evaluate_0.14 pillar_1.4.2
## [16] rlang_0.4.1 diptest_0.75-7 lazyeval_0.2.2
## [19] kernlab_0.9-27 Matrix_1.2-17 hash_2.2.6.1
## [22] rmarkdown_1.16 Maaslin2_1.0.0 labeling_0.3
## [25] splines_3.6.1 stringr_1.4.0 igraph_1.2.4.1
## [28] munsell_0.5.0 compiler_3.6.1 xfun_0.10
## [31] pkgconfig_2.0.3 biglm_0.9-1 mgcv_1.8-30
## [34] htmltools_0.4.0 nnet_7.3-12 tidyselect_0.2.5
## [37] tibble_2.1.3 bookdown_0.14 logging_0.10-108
## [40] crayon_1.3.4 withr_2.1.2 MASS_7.3-51.4
## [43] grid_3.6.1 nlme_3.1-141 gtable_0.3.0
## [46] lifecycle_0.1.0 metafor_2.1-0 scales_1.0.0
## [49] stringi_1.4.3 pbapply_1.4-2 reshape2_1.4.3
## [52] flexmix_2.3-15 getopt_1.20.3 lpsymphony_1.14.0
## [55] robustbase_0.93-5 optparse_1.6.4 ellipsis_0.3.0
## [58] vctrs_0.2.0 cowplot_1.0.0 tools_3.6.1
## [61] fpc_2.2-3 glue_1.3.1 DEoptimR_1.0-8
## [64] purrr_0.3.3 parallel_3.6.1 yaml_2.2.0
## [67] colorspace_1.4-1 cluster_2.1.0 BiocManager_1.30.9
## [70] prabclus_2.3-1 knitr_1.25 modeltools_0.2-22
Koren, Omry, Dan Knights, Antonio Gonzalez, Levi Waldron, Nicola Segata, Rob Knight, Curtis Huttenhower, and Ruth E Ley. 2013. “A Guide to Enterotypes Across the Human Body: Meta-Analysis of Microbial Community Structures in Human Microbiome Datasets.” PLoS Computational Biology 9 (1). Public Library of Science:e1002863.
Ravel, Jacques, Pawel Gajer, Zaid Abdo, G Maria Schneider, Sara SK Koenig, Stacey L McCulle, Shara Karlebach, et al. 2011. “Vaginal Microbiome of Reproductive-Age Women.” Proceedings of the National Academy of Sciences 108 (Supplement 1). National Acad Sciences:4680–7.
Tang, Zheng-Zheng, Guanhua Chen, and Alexander V Alekseyenko. 2016. “PERMANOVA-S: Association Test for Microbial Community Composition That Accommodates Confounders and Multiple Distances.” Bioinformatics 32 (17). Oxford University Press:2618–25.
Thomas, Andrew Maltez, Paolo Manghi, Francesco Asnicar, Edoardo Pasolli, Federica Armanini, Moreno Zolfo, Francesco Beghini, et al. 2019. “Metagenomic Analysis of Colorectal Cancer Datasets Identifies Cross-Cohort Microbial Diagnostic Signatures and a Link with Choline Degradation.” Nature Medicine 25 (4). Nature Publishing Group:667.
Tibshirani, Robert, and Guenther Walther. 2005. “Cluster Validation by Prediction Strength.” Journal of Computational and Graphical Statistics 14 (3). Taylor & Francis:511–28.