In FrederickHuangLin/ANCOMBC: Microbiome differential abudance and correlation analyses with bias correction
knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA,
fig.width = 6.25, fig.height = 5)
library(ANCOMBC)
library(tidyverse)
1. Introduction
Analysis of Composition of Microbiomes (ANCOM) [@mandal2015analysis] is a
differential abundance (DA) analysis for microbial absolute abundances.
It accounts for the compositionality of microbiome data by performing
the additive log ratio (ALR) transformation. ANCOM employs a heuristic strategy
to declare taxa that are significantly differentially abundant.
For a given taxon, the output W statistic represents the number ALR transformed
models where the taxon is differentially abundant with regard to the variable
of interest. The larger the value of W, the more likely the taxon is
differentially abundant. For more details, please refer to the
ANCOM paper.
2. Installation
Download package.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ANCOMBC")
Load the package.
library(ANCOMBC)
3. Run ANCOM on a real cross-sectional dataset {.tabset}
3.1 Import example data
The HITChip Atlas dataset contains genus-level microbiota profiling with
HITChip for 1006 western adults with no reported health complications,
reported in [@lahti2014tipping]. The dataset is available via the
microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq]
format. In this tutorial, we consider the following covariates:
-
Continuous covariates: "age"
-
Categorical covariates: "region", "bmi"
-
The group variable of interest: "bmi"
-
Three groups: "lean", "overweight", "obese"
-
The reference group: "obese"
data(atlas1006, package = "microbiome")
# Subset to baseline
pseq = phyloseq::subset_samples(atlas1006, time == 0)
# Re-code the bmi group
meta_data = microbiome::meta(pseq)
meta_data$bmi = recode(meta_data$bmi_group,
obese = "obese",
severeobese = "obese",
morbidobese = "obese")
# Note that by default, levels of a categorical variable in R are sorted
# alphabetically. In this case, the reference level for `bmi` will be
# `lean`. To manually change the reference level, for instance, setting `obese`
# as the reference level, use:
meta_data$bmi = factor(meta_data$bmi, levels = c("obese", "overweight", "lean"))
# You can verify the change by checking:
# levels(meta_data$bmi)
# Create the region variable
meta_data$region = recode(as.character(meta_data$nationality),
Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE",
CentralEurope = "CE", EasternEurope = "EE",
.missing = "unknown")
phyloseq::sample_data(pseq) = meta_data
# Subset to lean, overweight, and obese subjects
pseq = phyloseq::subset_samples(pseq, bmi %in% c("lean", "overweight", "obese"))
# Discard "EE" as it contains only 1 subject
# Discard subjects with missing values of region
pseq = phyloseq::subset_samples(pseq, ! region %in% c("EE", "unknown"))
print(pseq)
3.2 Run ancom function using phyloseq data
set.seed(123)
out = ancom(data = pseq, tax_level = "Family", meta_data = NULL,
p_adj_method = "holm", prv_cut = 0.10,
lib_cut = 1000, main_var = "bmi", adj_formula = "age + region",
rand_formula = NULL, lme_control = NULL, struc_zero = TRUE,
neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE)
res = out$res
# Similarly, if the main variable of interest is continuous, such as age, the
# ancom model can be specified as
# out = ancom(data = pseq, tax_level = "Family", meta_data = NULL,
# p_adj_method = "holm", prv_cut = 0.10,
# lib_cut = 1000, main_var = "age", adj_formula = "bmi + region",
# rand_formula = NULL, lme_control = NULL, struc_zero = FALSE,
# neg_lb = FALSE, alpha = 0.05, n_cl = 2, verbose = TRUE)
3.3 Scatter plot for W statistics
q_val = out$q_data
beta_val = out$beta_data
# Only consider the effect sizes with the corresponding q-value less than alpha
beta_val = beta_val * (q_val < 0.05)
# Choose the maximum of beta's as the effect size
beta_pos = apply(abs(beta_val), 2, which.max)
beta_max = vapply(seq_along(beta_pos), function(i)
beta_val[beta_pos[i], i], FUN.VALUE = double(1))
# Number of taxa except structural zeros
n_taxa = ifelse(is.null(out$zero_ind),
nrow(tse),
sum(apply(out$zero_ind[, -1], 1, sum) == 0))
# Cutoff values for declaring differentially abundant taxa
cut_off = 0.7 * (n_taxa - 1)
df_fig_w = res %>%
dplyr::mutate(beta = beta_max,
direct = case_when(
detected_0.7 == TRUE & beta > 0 ~ "Positive",
detected_0.7 == TRUE & beta <= 0 ~ "Negative",
TRUE ~ "Not Significant"
)) %>%
dplyr::arrange(W)
df_fig_w$taxon = factor(df_fig_w$taxon, levels = df_fig_w$taxon)
df_fig_w$W = replace(df_fig_w$W, is.infinite(df_fig_w$W), n_taxa - 1)
df_fig_w$direct = factor(df_fig_w$direct,
levels = c("Negative", "Positive", "Not Significant"))
p_w = df_fig_w %>%
ggplot(aes(x = taxon, y = W, color = direct)) +
geom_point(size = 2, alpha = 0.6) +
labs(x = "Taxon", y = "W") +
scale_color_discrete(name = NULL) +
geom_hline(yintercept = cut_off, linetype = "dotted",
color = "blue", size = 1.5) +
geom_text(aes(x = 2, y = cut_off + 0.5, label = "W[0.7]"),
size = 5, vjust = -0.5, hjust = 0, color = "orange", parse = TRUE) +
theme_bw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank())
p_w
3.4 Run ancom function using tse data
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(pseq)
set.seed(123)
out = ancom(data = tse, assay_name = "counts",
tax_level = "Family", meta_data = NULL,
p_adj_method = "holm", prv_cut = 0.10,
lib_cut = 1000, main_var = "bmi", adj_formula = "age + region",
rand_formula = NULL, lme_control = NULL, struc_zero = TRUE,
neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE)
res = out$res
3.5 Run ancom function by directly providing the abundance and metadata
abundance_data = microbiome::abundances(pseq)
aggregate_data = microbiome::abundances(microbiome::aggregate_taxa(pseq, "Family"))
meta_data = microbiome::meta(pseq)
set.seed(123)
out = ancom(data = abundance_data, aggregate_data = aggregate_data,
meta_data = meta_data, p_adj_method = "holm", prv_cut = 0.10,
lib_cut = 1000, main_var = "bmi", adj_formula = "age + region",
rand_formula = NULL, lme_control = NULL, struc_zero = TRUE,
neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE)
res = out$res
4. Run ANCOM on a real longitudinal dataset {.tabset}
4.1 Import example data
A two-week diet swap study between western (USA) and traditional (rural Africa)
diets [@lahti2014tipping]. The dataset is available via the
microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq]
format. In this tutorial, we consider the following fixed effects:
-
Continuous covariates: "timepoint"
-
Categorical covariates: "nationality"
-
The group variable of interest: "group"
-
Three groups: "DI", "ED", "HE"
-
The reference group: "DI"
and the following random effects:
-
A random intercept
-
A random slope: "timepoint"
data(dietswap, package = "microbiome")
print(dietswap)
4.2 Run ancom function using phyloseq data
set.seed(123)
out = ancom(data = dietswap, tax_level = "Family",
p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000,
main_var = "group",
adj_formula = "nationality + timepoint",
rand_formula = "(timepoint | subject)",
lme_control = lme4::lmerControl(),
struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2)
res = out$res
4.3 Visualization for W statistics
q_val = out$q_data
beta_val = out$beta_data
# Only consider the effect sizes with the corresponding q-value less than alpha
beta_val = beta_val * (q_val < 0.05)
# Choose the maximum of beta's as the effect size
beta_pos = apply(abs(beta_val), 2, which.max)
beta_max = vapply(seq_along(beta_pos), function(i) beta_val[beta_pos[i], i],
FUN.VALUE = double(1))
# Number of taxa except structural zeros
n_taxa = ifelse(is.null(out$zero_ind),
nrow(tse),
sum(apply(out$zero_ind[, -1], 1, sum) == 0))
# Cutoff values for declaring differentially abundant taxa
cut_off = 0.7 * (n_taxa - 1)
df_fig_w = res %>%
dplyr::mutate(beta = beta_max,
direct = case_when(
detected_0.7 == TRUE & beta > 0 ~ "Positive",
detected_0.7 == TRUE & beta <= 0 ~ "Negative",
TRUE ~ "Not Significant"
)) %>%
dplyr::arrange(W)
df_fig_w$taxon = factor(df_fig_w$taxon, levels = df_fig_w$taxon)
df_fig_w$W = replace(df_fig_w$W, is.infinite(df_fig_w$W), n_taxa - 1)
df_fig_w$direct = factor(df_fig_w$direct,
levels = c("Negative", "Positive", "Not Significant"))
p_w = df_fig_w %>%
ggplot(aes(x = taxon, y = W, color = direct)) +
geom_point(size = 2, alpha = 0.6) +
labs(x = "Taxon", y = "W") +
scale_color_discrete(name = NULL) +
geom_hline(yintercept = cut_off, linetype = "dotted",
color = "blue", size = 1.5) +
geom_text(aes(x = 2, y = cut_off + 0.5, label = "W[0.7]"),
size = 5, vjust = -0.5, hjust = 0, color = "orange", parse = TRUE) +
theme_bw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank())
p_w
4.4 Run ancom function using tse data
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(dietswap)
set.seed(123)
out = ancom(data = tse, assay_name = "counts", tax_level = "Family",
p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000,
main_var = "group",
adj_formula = "nationality + timepoint",
rand_formula = "(timepoint | subject)",
lme_control = lme4::lmerControl(),
struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2)
res = out$res
4.5 Run ancom function by directly providing the abundance and metadata
abundance_data = microbiome::abundances(dietswap)
aggregate_data = microbiome::abundances(microbiome::aggregate_taxa(dietswap, "Family"))
meta_data = microbiome::meta(dietswap)
set.seed(123)
out = ancom(data = abundance_data, aggregate_data = aggregate_data,
meta_data = meta_data, p_adj_method = "holm",
prv_cut = 0.10, lib_cut = 1000, main_var = "group",
adj_formula = "nationality + timepoint",
rand_formula = "(timepoint | subject)",
lme_control = lme4::lmerControl(),
struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2)
res = out$res
Session information
sessionInfo()
References
FrederickHuangLin/ANCOMBC documentation built on Jan. 8, 2025, 6:20 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA, fig.width = 6.25, fig.height = 5) library(ANCOMBC) library(tidyverse)
1. Introduction
Analysis of Composition of Microbiomes (ANCOM) [@mandal2015analysis] is a differential abundance (DA) analysis for microbial absolute abundances. It accounts for the compositionality of microbiome data by performing the additive log ratio (ALR) transformation. ANCOM employs a heuristic strategy to declare taxa that are significantly differentially abundant. For a given taxon, the output W statistic represents the number ALR transformed models where the taxon is differentially abundant with regard to the variable of interest. The larger the value of W, the more likely the taxon is differentially abundant. For more details, please refer to the ANCOM paper.
2. Installation
Download package.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("ANCOMBC")
Load the package.
library(ANCOMBC)
3. Run ANCOM on a real cross-sectional dataset {.tabset}
3.1 Import example data
The HITChip Atlas dataset contains genus-level microbiota profiling with HITChip for 1006 western adults with no reported health complications, reported in [@lahti2014tipping]. The dataset is available via the microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq] format. In this tutorial, we consider the following covariates:
-
Continuous covariates: "age"
-
Categorical covariates: "region", "bmi"
-
The group variable of interest: "bmi"
-
Three groups: "lean", "overweight", "obese"
-
The reference group: "obese"
-
data(atlas1006, package = "microbiome") # Subset to baseline pseq = phyloseq::subset_samples(atlas1006, time == 0) # Re-code the bmi group meta_data = microbiome::meta(pseq) meta_data$bmi = recode(meta_data$bmi_group, obese = "obese", severeobese = "obese", morbidobese = "obese") # Note that by default, levels of a categorical variable in R are sorted # alphabetically. In this case, the reference level for `bmi` will be # `lean`. To manually change the reference level, for instance, setting `obese` # as the reference level, use: meta_data$bmi = factor(meta_data$bmi, levels = c("obese", "overweight", "lean")) # You can verify the change by checking: # levels(meta_data$bmi) # Create the region variable meta_data$region = recode(as.character(meta_data$nationality), Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE", CentralEurope = "CE", EasternEurope = "EE", .missing = "unknown") phyloseq::sample_data(pseq) = meta_data # Subset to lean, overweight, and obese subjects pseq = phyloseq::subset_samples(pseq, bmi %in% c("lean", "overweight", "obese")) # Discard "EE" as it contains only 1 subject # Discard subjects with missing values of region pseq = phyloseq::subset_samples(pseq, ! region %in% c("EE", "unknown")) print(pseq)
3.2 Run ancom function using phyloseq data
set.seed(123) out = ancom(data = pseq, tax_level = "Family", meta_data = NULL, p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "bmi", adj_formula = "age + region", rand_formula = NULL, lme_control = NULL, struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE) res = out$res # Similarly, if the main variable of interest is continuous, such as age, the # ancom model can be specified as # out = ancom(data = pseq, tax_level = "Family", meta_data = NULL, # p_adj_method = "holm", prv_cut = 0.10, # lib_cut = 1000, main_var = "age", adj_formula = "bmi + region", # rand_formula = NULL, lme_control = NULL, struc_zero = FALSE, # neg_lb = FALSE, alpha = 0.05, n_cl = 2, verbose = TRUE)
3.3 Scatter plot for W statistics
q_val = out$q_data beta_val = out$beta_data # Only consider the effect sizes with the corresponding q-value less than alpha beta_val = beta_val * (q_val < 0.05) # Choose the maximum of beta's as the effect size beta_pos = apply(abs(beta_val), 2, which.max) beta_max = vapply(seq_along(beta_pos), function(i) beta_val[beta_pos[i], i], FUN.VALUE = double(1)) # Number of taxa except structural zeros n_taxa = ifelse(is.null(out$zero_ind), nrow(tse), sum(apply(out$zero_ind[, -1], 1, sum) == 0)) # Cutoff values for declaring differentially abundant taxa cut_off = 0.7 * (n_taxa - 1) df_fig_w = res %>% dplyr::mutate(beta = beta_max, direct = case_when( detected_0.7 == TRUE & beta > 0 ~ "Positive", detected_0.7 == TRUE & beta <= 0 ~ "Negative", TRUE ~ "Not Significant" )) %>% dplyr::arrange(W) df_fig_w$taxon = factor(df_fig_w$taxon, levels = df_fig_w$taxon) df_fig_w$W = replace(df_fig_w$W, is.infinite(df_fig_w$W), n_taxa - 1) df_fig_w$direct = factor(df_fig_w$direct, levels = c("Negative", "Positive", "Not Significant")) p_w = df_fig_w %>% ggplot(aes(x = taxon, y = W, color = direct)) + geom_point(size = 2, alpha = 0.6) + labs(x = "Taxon", y = "W") + scale_color_discrete(name = NULL) + geom_hline(yintercept = cut_off, linetype = "dotted", color = "blue", size = 1.5) + geom_text(aes(x = 2, y = cut_off + 0.5, label = "W[0.7]"), size = 5, vjust = -0.5, hjust = 0, color = "orange", parse = TRUE) + theme_bw() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), panel.grid.major = element_blank()) p_w
3.4 Run ancom function using tse data
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(pseq) set.seed(123) out = ancom(data = tse, assay_name = "counts", tax_level = "Family", meta_data = NULL, p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "bmi", adj_formula = "age + region", rand_formula = NULL, lme_control = NULL, struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE) res = out$res
3.5 Run ancom function by directly providing the abundance and metadata
abundance_data = microbiome::abundances(pseq) aggregate_data = microbiome::abundances(microbiome::aggregate_taxa(pseq, "Family")) meta_data = microbiome::meta(pseq) set.seed(123) out = ancom(data = abundance_data, aggregate_data = aggregate_data, meta_data = meta_data, p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "bmi", adj_formula = "age + region", rand_formula = NULL, lme_control = NULL, struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2, verbose = TRUE) res = out$res
4. Run ANCOM on a real longitudinal dataset {.tabset}
4.1 Import example data
A two-week diet swap study between western (USA) and traditional (rural Africa) diets [@lahti2014tipping]. The dataset is available via the microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq] format. In this tutorial, we consider the following fixed effects:
-
Continuous covariates: "timepoint"
-
Categorical covariates: "nationality"
-
The group variable of interest: "group"
-
Three groups: "DI", "ED", "HE"
-
The reference group: "DI"
-
and the following random effects:
-
A random intercept
-
A random slope: "timepoint"
data(dietswap, package = "microbiome") print(dietswap)
4.2 Run ancom function using phyloseq data
set.seed(123) out = ancom(data = dietswap, tax_level = "Family", p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "group", adj_formula = "nationality + timepoint", rand_formula = "(timepoint | subject)", lme_control = lme4::lmerControl(), struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2) res = out$res
4.3 Visualization for W statistics
q_val = out$q_data beta_val = out$beta_data # Only consider the effect sizes with the corresponding q-value less than alpha beta_val = beta_val * (q_val < 0.05) # Choose the maximum of beta's as the effect size beta_pos = apply(abs(beta_val), 2, which.max) beta_max = vapply(seq_along(beta_pos), function(i) beta_val[beta_pos[i], i], FUN.VALUE = double(1)) # Number of taxa except structural zeros n_taxa = ifelse(is.null(out$zero_ind), nrow(tse), sum(apply(out$zero_ind[, -1], 1, sum) == 0)) # Cutoff values for declaring differentially abundant taxa cut_off = 0.7 * (n_taxa - 1) df_fig_w = res %>% dplyr::mutate(beta = beta_max, direct = case_when( detected_0.7 == TRUE & beta > 0 ~ "Positive", detected_0.7 == TRUE & beta <= 0 ~ "Negative", TRUE ~ "Not Significant" )) %>% dplyr::arrange(W) df_fig_w$taxon = factor(df_fig_w$taxon, levels = df_fig_w$taxon) df_fig_w$W = replace(df_fig_w$W, is.infinite(df_fig_w$W), n_taxa - 1) df_fig_w$direct = factor(df_fig_w$direct, levels = c("Negative", "Positive", "Not Significant")) p_w = df_fig_w %>% ggplot(aes(x = taxon, y = W, color = direct)) + geom_point(size = 2, alpha = 0.6) + labs(x = "Taxon", y = "W") + scale_color_discrete(name = NULL) + geom_hline(yintercept = cut_off, linetype = "dotted", color = "blue", size = 1.5) + geom_text(aes(x = 2, y = cut_off + 0.5, label = "W[0.7]"), size = 5, vjust = -0.5, hjust = 0, color = "orange", parse = TRUE) + theme_bw() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), panel.grid.major = element_blank()) p_w
4.4 Run ancom function using tse data
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(dietswap) set.seed(123) out = ancom(data = tse, assay_name = "counts", tax_level = "Family", p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "group", adj_formula = "nationality + timepoint", rand_formula = "(timepoint | subject)", lme_control = lme4::lmerControl(), struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2) res = out$res
4.5 Run ancom function by directly providing the abundance and metadata
abundance_data = microbiome::abundances(dietswap) aggregate_data = microbiome::abundances(microbiome::aggregate_taxa(dietswap, "Family")) meta_data = microbiome::meta(dietswap) set.seed(123) out = ancom(data = abundance_data, aggregate_data = aggregate_data, meta_data = meta_data, p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000, main_var = "group", adj_formula = "nationality + timepoint", rand_formula = "(timepoint | subject)", lme_control = lme4::lmerControl(), struc_zero = TRUE, neg_lb = TRUE, alpha = 0.05, n_cl = 2) res = out$res
Session information
sessionInfo()
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.