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")
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(atlas1006)
# Subset to baseline
tse = tse[, tse$time == 0]
# Re-code the bmi group
tse$bmi = recode(tse$bmi_group,
obese = "obese",
severeobese = "obese",
morbidobese = "obese")
# Subset to lean, overweight, and obese subjects
tse = tse[, tse$bmi %in% c("lean", "overweight", "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:
tse$bmi = factor(tse$bmi, levels = c("obese", "overweight", "lean"))
# You can verify the change by checking:
# levels(sample_data(tse)$bmi)
# Create the region variable
tse$region = recode(as.character(tse$nationality),
Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE",
CentralEurope = "CE", EasternEurope = "EE",
.missing = "unknown")
# Discard "EE" as it contains only 1 subject
# Discard subjects with missing values of region
tse = tse[, ! tse$region %in% c("EE", "unknown")]
print(tse)class: TreeSummarizedExperiment
dim: 130 873
metadata(0):
assays(1): counts
rownames(130): Actinomycetaceae Aerococcus ... Xanthomonadaceae
Yersinia et rel.
rowData names(3): Phylum Family Genus
colnames(873): Sample-1 Sample-2 ... Sample-1005 Sample-1006
colData names(12): age sex ... bmi region
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
rowLinks: NULL
rowTree: NULL
colLinks: NULL
colTree: NULL3.2 Run ancom function
set.seed(123)
out = ancom(data = tse, assay_name = "counts",
tax_level = "Family", phyloseq = 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)
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 = tse, assay_name = "counts",
# tax_level = "Family", phyloseq = 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)
# ancom also supports importing data in phyloseq format
# tse_alt = agglomerateByRank(tse, "Family")
# pseq = makePhyloseqFromTreeSummarizedExperiment(tse_alt)
# out = ancom(data = NULL, assay_name = NULL,
# tax_level = "Family", phyloseq = pseq,
# 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)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
