3. Run ANCOM on a real cross-sectional dataset
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 (Lahti et al. 2014). The dataset is available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) 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)phyloseq-class experiment-level object
otu_table() OTU Table: [ 130 taxa and 873 samples ]
sample_data() Sample Data: [ 873 samples by 12 sample variables ]
tax_table() Taxonomy Table: [ 130 taxa by 3 taxonomic ranks ]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$res3.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