4 ANCOM-BC Implementation
4.1 Run ancombc function
out = ancombc(data = tse, assay_name = "counts",
tax_level = "Family", phyloseq = NULL,
formula = "age + region + bmi",
p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000,
group = "bmi", struc_zero = TRUE, neg_lb = TRUE, tol = 1e-5,
max_iter = 100, conserve = TRUE, alpha = 0.05, global = TRUE,
n_cl = 1, verbose = TRUE)
res = out$res
res_global = out$res_global
# ancombc also supports importing data in phyloseq format
# tse_alt = agglomerateByRank(tse, "Family")
# pseq = makePhyloseqFromTreeSummarizedExperiment(tse_alt)
# out = ancombc(data = NULL, assay_name = NULL,
# tax_level = "Family", phyloseq = pseq,
# formula = "age + region + bmi",
# p_adj_method = "holm", prv_cut = 0.10, lib_cut = 1000,
# group = "region", struc_zero = TRUE, neg_lb = TRUE, tol = 1e-5,
# max_iter = 100, conserve = TRUE, alpha = 0.05, global = TRUE,
# n_cl = 1, verbose = TRUE)4.2 ANCOMBC primary result
Result from the ANCOM-BC log-linear model to determine taxa that are differentially abundant according to the covariate of interest. It contains: 1) log fold changes; 2) standard errors; 3) test statistics; 4) p-values; 5) adjusted p-values; 6) indicators whether the taxon is differentially abundant (TRUE) or not (FALSE).
LFC
SE
Test statistic
P-values
Adjusted p-values
Differentially abundant taxa
Bias-corrected abundances
To obtain bias-corrected abundances, the following steps can be taken:
Step 1: Calculate the estimated sample-specific sampling fractions, in log scale.
Step 2: Correct the log observed abundances by subtracting the estimated sampling fraction from the log observed abundances of each sample.
It is important to note that we can only estimate sampling fractions up to an additive constant, meaning that only the difference between bias-corrected abundances is meaningful. Additionally, taxon-specific biases are not taken into account in the calculation of bias-corrected abundances, as it is assumed that these biases vary across taxa but remain constant across samples within a taxon.
samp_frac = out$samp_frac
# Replace NA with 0
samp_frac[is.na(samp_frac)] = 0
# Add pesudo-count (1) to avoid taking the log of 0
log_obs_abn = log(out$feature_table + 1)
# Adjust the log observed abundances
log_corr_abn = t(t(log_obs_abn) - samp_frac)
# Show the first 6 samples
round(log_corr_abn[, 1:6], 2) %>%
datatable(caption = "Bias-corrected log observed abundances")Visualization for age
df_lfc = data.frame(res$lfc[, -1] * res$diff_abn[, -1], check.names = FALSE) %>%
mutate(taxon_id = res$diff_abn$taxon) %>%
dplyr::select(taxon_id, everything())
df_se = data.frame(res$se[, -1] * res$diff_abn[, -1], check.names = FALSE) %>%
mutate(taxon_id = res$diff_abn$taxon) %>%
dplyr::select(taxon_id, everything())
colnames(df_se)[-1] = paste0(colnames(df_se)[-1], "SE")
df_fig_age = df_lfc %>%
dplyr::left_join(df_se, by = "taxon_id") %>%
dplyr::transmute(taxon_id, age, ageSE) %>%
dplyr::filter(age != 0) %>%
dplyr::arrange(desc(age)) %>%
dplyr::mutate(direct = ifelse(age > 0, "Positive LFC", "Negative LFC"))
df_fig_age$taxon_id = factor(df_fig_age$taxon_id, levels = df_fig_age$taxon_id)
df_fig_age$direct = factor(df_fig_age$direct,
levels = c("Positive LFC", "Negative LFC"))
p_age = ggplot(data = df_fig_age,
aes(x = taxon_id, y = age, fill = direct, color = direct)) +
geom_bar(stat = "identity", width = 0.7,
position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = age - ageSE, ymax = age + ageSE), width = 0.2,
position = position_dodge(0.05), color = "black") +
labs(x = NULL, y = "Log fold change",
title = "Log fold changes as one unit increase of age") +
scale_fill_discrete(name = NULL) +
scale_color_discrete(name = NULL) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
panel.grid.minor.y = element_blank(),
axis.text.x = element_text(angle = 60, hjust = 1))
p_age
Visualization for BMI
df_fig_bmi = df_lfc %>%
filter(bmioverweight != 0 | bmilean != 0) %>%
transmute(taxon_id,
`Overweight vs. Obese` = round(bmioverweight, 2),
`Lean vs. Obese` = round(bmilean, 2)) %>%
pivot_longer(cols = `Overweight vs. Obese`:`Lean vs. Obese`,
names_to = "group", values_to = "value") %>%
arrange(taxon_id)
lo = floor(min(df_fig_bmi$value))
up = ceiling(max(df_fig_bmi$value))
mid = (lo + up)/2
p_bmi = df_fig_bmi %>%
ggplot(aes(x = group, y = taxon_id, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
na.value = "white", midpoint = mid, limit = c(lo, up),
name = NULL) +
geom_text(aes(group, taxon_id, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Log fold changes as compared to obese subjects") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
p_bmi
4.3 ANCOMBC global test result
Result from the ANCOM-BC global test to determine taxa that are differentially abundant between at least two groups across three or more different groups. In this example, we want to identify taxa that are differentially abundant between at least two regions across CE, NE, SE, and US. The result contains: 1) test statistics; 2) p-values; 3) adjusted p-values; 4) indicators whether the taxon is differentially abundant (TRUE) or not (FALSE).
Test statistics
P-values
Adjusted p-values
Differentially abundant taxa
Visualization
sig_taxa = res_global %>%
dplyr::filter(diff_abn == TRUE) %>%
.$taxon
df_bmi = tab_lfc %>%
dplyr::select(Taxon, `Overweight - Obese`, `Lean - Obese`) %>%
filter(Taxon %in% sig_taxa)
df_heat = df_bmi %>%
pivot_longer(cols = -one_of("Taxon"),
names_to = "region", values_to = "value") %>%
mutate(value = round(value, 2))
df_heat$Taxon = factor(df_heat$Taxon, levels = sort(sig_taxa))
lo = floor(min(df_heat$value))
up = ceiling(max(df_heat$value))
mid = (lo + up)/2
p_heat = df_heat %>%
ggplot(aes(x = region, y = Taxon, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
na.value = "white", midpoint = mid, limit = c(lo, up),
name = NULL) +
geom_text(aes(region, Taxon, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL,
title = "Log fold changes for globally significant taxa") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
p_heat