ANCOM Tutorial

Huang Lin\(^1\)

\(^1\)NICHD, 6710B Rockledge Dr, Bethesda, MD 20817

December 18, 2022

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA,
                      fig.width = 6.25, fig.height = 5)
library(ANCOMBC)
library(tidyverse)

1. Introduction

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

2. Installation

Download package.

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

Load the package.

library(ANCOMBC)

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 also 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:

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: NULL

3.2 Run ancom function

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

4. Run ANCOM on a real longitudinal dataset

4.1 Import example data

A two-week diet swap study between western (USA) and traditional (rural Africa) diets (Lahti et al. 2014). The dataset is also available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) format. In this tutorial, we consider the following fixed effects:

and the following random effects:

class: TreeSummarizedExperiment 
dim: 130 222 
metadata(0):
assays(1): counts
rownames(130): Actinomycetaceae Aerococcus ... Xanthomonadaceae
  Yersinia et rel.
rowData names(3): Phylum Family Genus
colnames(222): Sample-1 Sample-2 ... Sample-221 Sample-222
colData names(8): subject sex ... timepoint.within.group bmi_group
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
rowLinks: NULL
rowTree: NULL
colLinks: NULL
colTree: NULL

4.2 Run ancom function

4.3 Visualization 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

Session information

sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB              LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] forcats_0.5.2   stringr_1.5.0   dplyr_1.0.10    purrr_0.3.5    
 [5] readr_2.1.3     tidyr_1.2.1     tibble_3.1.8    ggplot2_3.4.0  
 [9] tidyverse_1.3.2 ANCOMBC_2.0.2  

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                     tidyselect_1.2.0              
  [3] lme4_1.1-31                    RSQLite_2.2.19                
  [5] htmlwidgets_1.6.0              grid_4.2.2                    
  [7] BiocParallel_1.32.4            gmp_0.6-9                     
  [9] munsell_0.5.0                  ScaledMatrix_1.6.0            
 [11] codetools_0.2-18               interp_1.1-3                  
 [13] withr_2.5.0                    colorspace_2.0-3              
 [15] Biobase_2.58.0                 energy_1.7-10                 
 [17] highr_0.9                      knitr_1.41                    
 [19] rstudioapi_0.14                stats4_4.2.2                  
 [21] SingleCellExperiment_1.20.0    DescTools_0.99.47             
 [23] labeling_0.4.2                 MatrixGenerics_1.10.0         
 [25] Rdpack_2.4                     emmeans_1.8.3                 
 [27] GenomeInfoDbData_1.2.9         farver_2.1.1                  
 [29] bit64_4.0.5                    coda_0.19-4                   
 [31] vctrs_0.5.1                    treeio_1.22.0                 
 [33] generics_0.1.3                 TH.data_1.1-1                 
 [35] xfun_0.35                      timechange_0.1.1              
 [37] R6_2.5.1                       doParallel_1.0.17             
 [39] GenomeInfoDb_1.34.4            ggbeeswarm_0.7.1              
 [41] rsvd_1.0.5                     bitops_1.0-7                  
 [43] cachem_1.0.6                   DelayedArray_0.24.0           
 [45] assertthat_0.2.1               scales_1.2.1                  
 [47] multcomp_1.4-20                nnet_7.3-18                   
 [49] googlesheets4_1.0.1            beeswarm_0.4.0                
 [51] rootSolve_1.8.2.3              gtable_0.3.1                  
 [53] beachmat_2.14.0                lmom_2.9                      
 [55] sandwich_3.0-2                 rlang_1.0.6                   
 [57] splines_4.2.2                  lazyeval_0.2.2                
 [59] gargle_1.2.1                   broom_1.0.2                   
 [61] checkmate_2.1.0                modelr_0.1.10                 
 [63] yaml_2.3.6                     reshape2_1.4.4                
 [65] backports_1.4.1                Hmisc_4.7-2                   
 [67] tools_4.2.2                    ellipsis_0.3.2                
 [69] decontam_1.18.0                jquerylib_0.1.4               
 [71] RColorBrewer_1.1-3             proxy_0.4-27                  
 [73] BiocGenerics_0.44.0            MultiAssayExperiment_1.24.0   
 [75] Rcpp_1.0.9                     plyr_1.8.8                    
 [77] base64enc_0.1-3                sparseMatrixStats_1.10.0      
 [79] zlibbioc_1.44.0                RCurl_1.98-1.9                
 [81] rpart_4.1.19                   TreeSummarizedExperiment_2.6.0
 [83] deldir_1.0-6                   viridis_0.6.2                 
 [85] S4Vectors_0.36.1               zoo_1.8-11                    
 [87] SummarizedExperiment_1.28.0    haven_2.5.1                   
 [89] ggrepel_0.9.2                  cluster_2.1.4                 
 [91] fs_1.5.2                       DECIPHER_2.26.0               
 [93] magrittr_2.0.3                 data.table_1.14.6             
 [95] lmerTest_3.1-3                 reprex_2.0.2                  
 [97] googledrive_2.0.0              mvtnorm_1.1-3                 
 [99] matrixStats_0.63.0             gsl_2.1-7.1                   
[101] hms_1.1.2                      evaluate_0.19                 
[103] xtable_1.8-4                   jpeg_0.1-10                   
[105] readxl_1.4.1                   IRanges_2.32.0                
[107] gridExtra_2.3                  compiler_4.2.2                
[109] scater_1.26.1                  crayon_1.5.2                  
[111] minqa_1.2.5                    htmltools_0.5.4               
[113] tzdb_0.3.0                     mgcv_1.8-41                   
[115] Formula_1.2-4                  expm_0.999-6                  
[117] Exact_3.2                      lubridate_1.9.0               
[119] DBI_1.1.3                      dbplyr_2.2.1                  
[121] MASS_7.3-58.1                  boot_1.3-28.1                 
[123] Matrix_1.5-3                   permute_0.9-7                 
[125] cli_3.4.1                      rbibutils_2.2.11              
[127] parallel_4.2.2                 GenomicRanges_1.50.2          
[129] pkgconfig_2.0.3                numDeriv_2016.8-1.1           
[131] foreign_0.8-84                 scuttle_1.8.3                 
[133] xml2_1.3.3                     foreach_1.5.2                 
[135] vipor_0.4.5                    bslib_0.4.2                   
[137] DirichletMultinomial_1.40.0    rngtools_1.5.2                
[139] XVector_0.38.0                 estimability_1.4.1            
[141] mia_1.6.0                      rvest_1.0.3                   
[143] CVXR_1.0-11                    doRNG_1.8.2                   
[145] yulab.utils_0.0.5              digest_0.6.31                 
[147] vegan_2.6-4                    Biostrings_2.66.0             
[149] rmarkdown_2.19                 cellranger_1.1.0              
[151] tidytree_0.4.2                 htmlTable_2.4.1               
[153] gld_2.6.6                      DelayedMatrixStats_1.20.0     
[155] nloptr_2.0.3                   lifecycle_1.0.3               
[157] nlme_3.1-161                   jsonlite_1.8.4                
[159] BiocNeighbors_1.16.0           viridisLite_0.4.1             
[161] fansi_1.0.3                    pillar_1.8.1                  
[163] lattice_0.20-45                fastmap_1.1.0                 
[165] httr_1.4.4                     survival_3.4-0                
[167] glue_1.6.2                     png_0.1-8                     
[169] iterators_1.0.14               bit_4.0.5                     
[171] class_7.3-20                   stringi_1.7.8                 
[173] sass_0.4.4                     blob_1.2.3                    
[175] BiocSingular_1.14.0            latticeExtra_0.6-30           
[177] memoise_2.0.1                  Rmpfr_0.8-9                   
[179] irlba_2.3.5.1                  e1071_1.7-12                  
[181] ape_5.6-2                     

References

Lahti, Leo, Jarkko Salojärvi, Anne Salonen, Marten Scheffer, and Willem M De Vos. 2014. “Tipping Elements in the Human Intestinal Ecosystem.” Nature Communications 5 (1): 1–10.

Lahti, Leo, Sudarshan Shetty, T Blake, J Salojarvi, and others. 2017. “Tools for Microbiome Analysis in R.” Version 1: 10013.

Mandal, Siddhartha, Will Van Treuren, Richard A White, Merete Eggesbø, Rob Knight, and Shyamal D Peddada. 2015. “Analysis of Composition of Microbiomes: A Novel Method for Studying Microbial Composition.” Microbial Ecology in Health and Disease 26 (1): 27663.

McMurdie, Paul J, and Susan Holmes. 2013. “Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data.” PloS One 8 (4): e61217.