1. Introduction

Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) (Lin and Peddada 2020) is a methodology of differential abundance (DA) analysis for microbial absolute abundances. ANCOM-BC estimates the unknown sampling fractions, corrects the bias induced by their differences through a log linear regression model including the estimated sampling fraction as an offset terms, and identifies taxa that are differentially abundant according to the variable of interest. For more details, please refer to the ANCOM-BC paper.

Analysis of Composition of Microbiomes (ANCOM) (Mandal et al. 2015) is also a 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. Larger the value of W, the more likely the taxon is differentially abundant. For more details, please refer to the ANCOM paper.

Sparse Estimation of Correlations among Microbiomes (SECOM) is a methodology which aims to detect both linear and nonlinear relationships between a pair of taxa within an ecosystem (e.g. gut) or across ecosystems (e.g. gut and tongue). SECOM corrects both sample-specific and taxon-specific biases, obtains a consistent estimator for the correlation matrix of microbial absolute abundances, while maintaining the underlying true sparsity.

2. Installation

Download package.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ANCOMBC")

Load the package.

library(ANCOMBC)

3. Example Data

Cross-sectional data

The HITChip Atlas data set (Lahti et al. 2014) is available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) format.

data(atlas1006) 
# Subset to baseline
pseq = subset_samples(atlas1006, time == 0)

# Re-code the bmi group
sample_data(pseq)$bmi_group = recode(sample_data(pseq)$bmi_group,
                                     lean = "lean",
                                     overweight = "overweight",
                                     obese = "obese",
                                     severeobese = "obese",
                                     morbidobese = "obese")
# Subset to lean and obese subjects
pseq = subset_samples(pseq, bmi_group %in% c("lean", "obese"))

# Create the region variable
sample_data(pseq)$region = recode(sample_data(pseq)$nationality,
                                  Scandinavia = "NE",
                                  UKIE = "NE",
                                  SouthEurope = "SE",
                                  CentralEurope = "CE",
                                  EasternEurope = "EE")

# Discard "EE" as it contains only 1 subject
pseq = subset_samples(pseq, region != "EE")
# Genus level data
genus_data = aggregate_taxa(pseq, "Genus")
# Family level data
family_data = aggregate_taxa(pseq, "Family")
# Phylum level data
phylum_data = aggregate_taxa(pseq, "Phylum")

print(genus_data)
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 130 taxa and 722 samples ]
sample_data() Sample Data:       [ 722 samples by 11 sample variables ]
tax_table()   Taxonomy Table:    [ 130 taxa by 3 taxonomic ranks ]
print(family_data)
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 22 taxa and 722 samples ]
sample_data() Sample Data:       [ 722 samples by 11 sample variables ]
tax_table()   Taxonomy Table:    [ 22 taxa by 3 taxonomic ranks ]
print(phylum_data)
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 8 taxa and 722 samples ]
sample_data() Sample Data:       [ 722 samples by 11 sample variables ]
tax_table()   Taxonomy Table:    [ 8 taxa by 2 taxonomic ranks ]

Longitudinal data

A two-week diet swap study between western (USA) and traditional (rural Africa) diets (Lahti et al. 2014). The data set is available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) format.

data(dietswap)

# Aggregate to family level
family_data2 = aggregate_taxa(dietswap, "Family")

print(family_data2)
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 22 taxa and 222 samples ]
sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
tax_table()   Taxonomy Table:    [ 22 taxa by 3 taxonomic ranks ]

4. ANCOMBC Implementation

4.1 Run ancombc function

out = ancombc(phyloseq = family_data, formula = "age + region + bmi_group", 
              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)

res = out$res
res_global = out$res_global

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

tab_lfc = res$lfc
col_name = c("Age", "NE - CE", "SE - CE", "US - CE", "Obese - Lean")
colnames(tab_lfc) = col_name
tab_lfc %>% 
  datatable(caption = "Log Fold Changes from the Primary Result") %>%
  formatRound(col_name, digits = 2)

SE

tab_se = res$se
colnames(tab_se) = col_name
tab_se %>% 
  datatable(caption = "SEs from the Primary Result") %>%
  formatRound(col_name, digits = 2)

Test statistic

tab_w = res$W
colnames(tab_w) = col_name
tab_w %>% 
  datatable(caption = "Test Statistics from the Primary Result") %>%
  formatRound(col_name, digits = 2)

P-values

tab_p = res$p_val
colnames(tab_p) = col_name
tab_p %>% 
  datatable(caption = "P-values from the Primary Result") %>%
  formatRound(col_name, digits = 2)

Adjusted p-values

tab_q = res$q
colnames(tab_q) = col_name
tab_q %>% 
  datatable(caption = "Adjusted p-values from the Primary Result") %>%
  formatRound(col_name, digits = 2)

Differentially abundant taxa

tab_diff = res$diff_abn
colnames(tab_diff) = col_name
tab_diff %>% 
  datatable(caption = "Differentially Abundant Taxa from the Primary Result")

Visualization for age

df_lfc = data.frame(res$lfc * res$diff_abn, check.names = FALSE) %>% 
  rownames_to_column("taxon_id")
df_se = data.frame(res$se * res$diff_abn, check.names = FALSE) %>% 
  rownames_to_column("taxon_id")
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 = "Waterfall Plot 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