4. Run ANCOM-BC2 on a real cross-sectional dataset
4.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: NULL4.2 Run ancombc2 function
set.seed(123)
# It should be noted that we have set the number of bootstrap samples (B) equal
# to 10 in the 'trend_control' function for computational expediency.
# However, it is recommended that users utilize the default value of B,
# which is 100, or larger values for optimal performance.
output = ancombc2(data = tse, assay_name = "counts", tax_level = "Family",
fix_formula = "age + region + bmi", rand_formula = NULL,
p_adj_method = "holm",
prv_cut = 0.10, lib_cut = 1000, s0_perc = 0.05,
group = "bmi", struc_zero = TRUE, neg_lb = TRUE,
alpha = 0.05, n_cl = 2, verbose = TRUE,
global = TRUE, pairwise = TRUE, dunnet = TRUE, trend = TRUE,
iter_control = list(tol = 1e-2, max_iter = 20,
verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
nrow = 2,
byrow = TRUE),
matrix(c(-1, 0, 1, -1),
nrow = 2,
byrow = TRUE)),
node = list(2, 2),
solver = "ECOS",
B = 10))4.3 Structural zeros (taxon presence/absence)
4.4 ANCOM-BC2 primary analysis
Result from the ANCOM-BC2 methodology to determine taxa that are differentially abundant according to the covariate of interest. Results contain: 1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the taxon is differentially abundant (TRUE) or not (FALSE).
Results for age
df_age = res_prim %>%
dplyr::select(taxon, ends_with("age"))
df_fig_age = df_age %>%
dplyr::filter(diff_age == 1, passed_ss_age == 1) %>%
dplyr::arrange(desc(lfc_age)) %>%
dplyr::mutate(direct = ifelse(lfc_age > 0, "Positive LFC", "Negative LFC"))
df_fig_age$taxon = factor(df_fig_age$taxon, levels = df_fig_age$taxon)
df_fig_age$direct = factor(df_fig_age$direct,
levels = c("Positive LFC", "Negative LFC"))
fig_age = df_fig_age %>%
ggplot(aes(x = taxon, y = lfc_age, fill = direct)) +
geom_bar(stat = "identity", width = 0.7, color = "black",
position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = lfc_age - se_age, ymax = lfc_age + se_age),
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))
fig_age
Results for bmi
df_bmi = res_prim %>%
dplyr::select(taxon, contains("bmi"))
df_fig_bmi = df_bmi %>%
dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) |
(diff_bmioverweight == 1 & passed_ss_bmioverweight == 1)) %>%
dplyr::mutate(lfc_overweight = ifelse(diff_bmioverweight == 1,
lfc_bmioverweight, 0),
lfc_lean = ifelse(diff_bmilean == 1,
lfc_bmilean, 0)) %>%
dplyr::transmute(taxon,
`Overweight vs. Obese` = round(lfc_overweight, 2),
`Lean vs. Obese` = round(lfc_lean, 2)) %>%
tidyr::pivot_longer(cols = `Overweight vs. Obese`:`Lean vs. Obese`,
names_to = "group", values_to = "value") %>%
dplyr::arrange(taxon)
lo = floor(min(df_fig_bmi$value))
up = ceiling(max(df_fig_bmi$value))
mid = (lo + up)/2
fig_bmi = df_fig_bmi %>%
ggplot(aes(x = group, 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(group, taxon, 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))
fig_bmi
4.5 ANCOM-BC2 global test
ANCOM-BC2 global test aims to determine taxa that are differentially abundant between at least two groups across three or more experimental groups.
In this example, we want to identify taxa that are differentially abundant between at least two groups across “lean”, “overweight”, and “obese”. The result contains: 1) test statistics, 2) p-values, 3) adjusted p-values, 4) indicators of whether the taxon is differentially abundant (TRUE) or not (FALSE).
res_global = output$res_global
df_bmi = res_prim %>%
dplyr::select(taxon, contains("bmi"))
df_fig_global = df_bmi %>%
dplyr::left_join(res_global %>%
dplyr::transmute(taxon,
diff_bmi = diff_abn,
passed_ss = passed_ss)) %>%
dplyr::filter(diff_bmi == 1, passed_ss == 1) %>%
dplyr::mutate(lfc_lean = lfc_bmilean,
lfc_overweight = lfc_bmioverweight) %>%
dplyr::transmute(taxon,
`Lean vs. Obese` = round(lfc_lean, 2),
`Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
tidyr::pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`,
names_to = "group", values_to = "value") %>%
dplyr::arrange(taxon)
lo = floor(min(df_fig_global$value))
up = ceiling(max(df_fig_global$value))
mid = (lo + up)/2
fig_global = df_fig_global %>%
ggplot(aes(x = group, 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(group, 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))
fig_global
4.6 ANCOM-BC2 multiple pairwise comparisons
ANCOM-BC2 multiple pairwise comparisons aim to determine taxa that are differentially abundant between any pair of two groups across three or more experimental groups, while controlling the mdFDR.
In this example, we want to identify taxa that are differentially abundant between any pair of two groups across “lean”, “overweight”, and “obese”. The result contains: 1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the taxon is differentially abundant (TRUE) or not (FALSE).
res_pair = output$res_pair
df_fig_pair = res_pair %>%
dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) |
(diff_bmioverweight == 1 & passed_ss_bmioverweight == 1) |
(diff_bmilean_bmioverweight == 1 & passed_ss_bmilean_bmioverweight == 1)) %>%
dplyr::mutate(lfc_lean = ifelse(diff_bmilean == 1,
lfc_bmilean, 0),
lfc_overweight = ifelse(diff_bmioverweight == 1,
lfc_bmioverweight, 0),
lfc_lean_overweight = ifelse(diff_bmilean_bmioverweight == 1,
lfc_bmilean_bmioverweight, 0)) %>%
dplyr::transmute(taxon,
`Lean vs. Obese` = round(lfc_lean, 2),
`Overweight vs. Obese` = round(lfc_overweight, 2),
`Lean vs. Overweight` = round(lfc_lean_overweight, 2)
) %>%
tidyr::pivot_longer(cols = `Lean vs. Obese`:`Lean vs. Overweight`,
names_to = "group", values_to = "value") %>%
dplyr::arrange(taxon)
df_fig_pair$group = factor(df_fig_pair$group,
levels = c("Lean vs. Obese",
"Overweight vs. Obese",
"Lean vs. Overweight"))
lo = floor(min(df_fig_pair$value))
up = ceiling(max(df_fig_pair$value))
mid = (lo + up)/2
fig_pair = df_fig_pair %>%
ggplot(aes(x = group, 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(group, taxon, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Log fold change of pairwise comparisons") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
fig_pair
4.7 ANCOM-BC2 multiple pairwise comparisons against a pre-specified group (Dunnett’s type of test)
The Dunnett’s test (Dunnett 1955; Dunnett and Tamhane 1991, 1992) is designed for making comparisons of several experimental groups with the control or the reference group. ANCOM-BC2 Dunnett’s type of test adopts the framework of Dunnett’s test while controlling the mdFDR. Of note is that the ANCOM-BC2 primary results do not control the mdFDR for the comparison of multiple groups.
In this example, we want to identify taxa that are differentially abundant between “lean”, “overweight”, and the reference group “obese”. The result contains: 1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the taxon is differentially abundant (TRUE) or not (FALSE).
res_dunn = output$res_dunn
df_fig_dunn = res_dunn %>%
dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) |
(diff_bmioverweight == 1 & passed_ss_bmioverweight == 1)) %>%
dplyr::mutate(lfc_lean = ifelse(diff_bmilean == 1, lfc_bmilean, 0),
lfc_overweight = ifelse(diff_bmioverweight == 1,
lfc_bmioverweight, 0)) %>%
dplyr::transmute(taxon,
`Lean vs. Obese` = round(lfc_lean, 2),
`Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
tidyr::pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`,
names_to = "group", values_to = "value") %>%
dplyr::arrange(taxon)
lo = floor(min(df_fig_dunn$value))
up = ceiling(max(df_fig_dunn$value))
mid = (lo + up)/2
fig_dunn = df_fig_dunn %>%
ggplot(aes(x = group, 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(group, taxon, 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))
fig_dunn
4.8 ANCOM-BC2 pattern analysis
Sometimes the experimental groups are intrinsically ordered (e.g., in a dose-response study), and we would like to see if the microbial abundances show some patterns accordingly. Examples of patterns include: monotonically increasing, monotonically decreasing, and umbrella shape.
ANCOM-BC2 is able to identify potential patterns by testing the contrast: \[Ax \ge 0\] where \(A\) is the contrast matrix and \(x\) is the vector of parameters.
For instance, in this example, we want to identify taxa that are monotonically increasing across “obese”, “overweight”, and “lean”. Note that “obese” is the reference group, and the parameters we can estimate are the differences as compared to the reference group, i.e., \[x = (\text{overweight - obese}, \text{lean - obese})^T\] To test the monotonically increasing trend: \[H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \le \text{overweight} \le \text{lean} \quad \text{with at least one strict inequality}\] We shall specify the contrast matrix \(A\) as \[A = \begin{bmatrix} 1 & 0 \\ -1 & 1 \end{bmatrix}\] Similarly, to test for the monotonically decreasing trend: \[H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \ge \text{overweight} \ge \text{lean} \quad \text{with at least one strict inequality}\] We shall specify the contrast matrix \(A\) as \[A = \begin{bmatrix} -1 & 0 \\ 1 & -1 \end{bmatrix}\] To test for monotonic trend (increasing or decreasing), one should specify the node parameter in trend_control as the last position of x. In this example, the vector of parameters \(x\) is of length 2, thus, the last position is 2. For testing umbrella shape, for instance, in this example: \[H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \le \text{overweight} \ge \text{lean} \quad \text{with at least one strict inequality}\] one should set node as the position of the turning point of x. In this example, the turning position is overweight, thus, node = 1.
We will test both the monotonically increasing and decreasing trends in this example. The result contains: 1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the taxon is differentially abundant (TRUE) or not (FALSE).
Note that the LFC and SE for the reference group (“obese” in this example) will set to be 0s.
res_trend = output$res_trend
df_fig_trend = res_trend %>%
dplyr::filter(diff_abn == 1, passed_ss == 1) %>%
dplyr::transmute(taxon,
lfc = lfc_bmioverweight,
se = se_bmioverweight,
q_val,
group = "Overweight - Obese") %>%
dplyr::bind_rows(
res_trend %>%
dplyr::filter(diff_abn, passed_ss == 1) %>%
dplyr::transmute(taxon,
lfc = lfc_bmilean,
se = se_bmilean,
q_val,
group = "Lean - Obese")
)
df_fig_trend$group = factor(df_fig_trend$group,
levels = c("Overweight - Obese", "Lean - Obese"))
fig_trend = df_fig_trend %>%
ggplot(aes(x = group, y = lfc, fill = group)) +
geom_bar(stat = "identity", position = position_dodge(), color = "black") +
geom_errorbar(aes(ymin = lfc - se, ymax = lfc + se), width = .2,
position = position_dodge(.9)) +
facet_wrap(vars(taxon), nrow = 2, scales = "free") +
labs(x = NULL, y = NULL, title = "Log fold change as compared to obese subjects") +
scale_fill_brewer(palette = "Set2", name = NULL) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "bottom")
fig_trend