ANCOM-BC Tutorial

Huang Lin\(^1\)

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

October 18, 2022

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.

2. Installation

Download package.

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

Load the package.

library(ANCOMBC)

3. 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:

data(atlas1006)

# Subset to baseline
tse = atlas1006[, atlas1006$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: NULL

4 ANCOM-BC Implementation

4.1 Run ancombc function

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

Step 1: obtain estimated sample-specific sampling fractions (in log scale).

Step 2: correct the log observed abundances by subtracting the estimated sampling fraction from log observed abundances of each sample.

Note that we are only able to estimate sampling fractions up to an additive constant. Thus, only the difference between bias-corrected abundances are meaningful.

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

Session information

sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.15-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] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] DT_0.25                        mia_1.4.0                     
 [3] MultiAssayExperiment_1.22.0    TreeSummarizedExperiment_2.4.0
 [5] Biostrings_2.64.1              XVector_0.36.0                
 [7] SingleCellExperiment_1.18.1    SummarizedExperiment_1.26.1   
 [9] Biobase_2.56.0                 GenomicRanges_1.48.0          
[11] GenomeInfoDb_1.32.4            IRanges_2.30.1                
[13] S4Vectors_0.34.0               BiocGenerics_0.42.0           
[15] MatrixGenerics_1.8.1           matrixStats_0.62.0            
[17] forcats_0.5.2                  stringr_1.4.1                 
[19] dplyr_1.0.10                   purrr_0.3.5                   
[21] readr_2.1.3                    tidyr_1.2.1                   
[23] tibble_3.1.8                   ggplot2_3.3.6                 
[25] tidyverse_1.3.2                ANCOMBC_1.6.4                 

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                  tidyselect_1.2.0           
  [3] lme4_1.1-30                 RSQLite_2.2.18             
  [5] htmlwidgets_1.5.4           grid_4.2.1                 
  [7] BiocParallel_1.30.4         gmp_0.6-6                  
  [9] munsell_0.5.0               ScaledMatrix_1.4.1         
 [11] codetools_0.2-18            interp_1.1-3               
 [13] withr_2.5.0                 colorspace_2.0-3           
 [15] energy_1.7-10               highr_0.9                  
 [17] knitr_1.40                  rstudioapi_0.14            
 [19] DescTools_0.99.46           labeling_0.4.2             
 [21] Rdpack_2.4                  emmeans_1.8.1-1            
 [23] GenomeInfoDbData_1.2.8      farver_2.1.1               
 [25] bit64_4.0.5                 coda_0.19-4                
 [27] vctrs_0.4.2                 treeio_1.20.2              
 [29] generics_0.1.3              TH.data_1.1-1              
 [31] xfun_0.33                   R6_2.5.1                   
 [33] doParallel_1.0.17           ggbeeswarm_0.6.0           
 [35] rsvd_1.0.5                  bitops_1.0-7               
 [37] cachem_1.0.6                DelayedArray_0.22.0        
 [39] assertthat_0.2.1            scales_1.2.1               
 [41] googlesheets4_1.0.1         multcomp_1.4-20            
 [43] nnet_7.3-18                 beeswarm_0.4.0             
 [45] rootSolve_1.8.2.3           gtable_0.3.1               
 [47] beachmat_2.12.0             lmom_2.9                   
 [49] sandwich_3.0-2              rlang_1.0.6                
 [51] splines_4.2.1               lazyeval_0.2.2             
 [53] gargle_1.2.1                broom_1.0.1                
 [55] checkmate_2.1.0             modelr_0.1.9               
 [57] yaml_2.3.6                  reshape2_1.4.4             
 [59] crosstalk_1.2.0             backports_1.4.1            
 [61] Hmisc_4.7-1                 tools_4.2.1                
 [63] ellipsis_0.3.2              decontam_1.16.0            
 [65] jquerylib_0.1.4             RColorBrewer_1.1-3         
 [67] proxy_0.4-27                Rcpp_1.0.9                 
 [69] plyr_1.8.7                  base64enc_0.1-3            
 [71] sparseMatrixStats_1.8.0     zlibbioc_1.42.0            
 [73] RCurl_1.98-1.9              rpart_4.1.16               
 [75] deldir_1.0-6                viridis_0.6.2              
 [77] zoo_1.8-11                  haven_2.5.1                
 [79] ggrepel_0.9.1               cluster_2.1.4              
 [81] fs_1.5.2                    DECIPHER_2.24.0            
 [83] magrittr_2.0.3              data.table_1.14.4          
 [85] lmerTest_3.1-3              reprex_2.0.2               
 [87] googledrive_2.0.0           mvtnorm_1.1-3              
 [89] gsl_2.1-7.1                 hms_1.1.2                  
 [91] evaluate_0.17               xtable_1.8-4               
 [93] jpeg_0.1-9                  readxl_1.4.1               
 [95] gridExtra_2.3               compiler_4.2.1             
 [97] scater_1.24.0               crayon_1.5.2               
 [99] minqa_1.2.4                 htmltools_0.5.3            
[101] tzdb_0.3.0                  mgcv_1.8-40                
[103] Formula_1.2-4               expm_0.999-6               
[105] Exact_3.2                   lubridate_1.8.0            
[107] DBI_1.1.3                   dbplyr_2.2.1               
[109] MASS_7.3-58.1               boot_1.3-28                
[111] Matrix_1.5-1                permute_0.9-7              
[113] cli_3.4.1                   rbibutils_2.2.9            
[115] parallel_4.2.1              pkgconfig_2.0.3            
[117] numDeriv_2016.8-1.1         foreign_0.8-83             
[119] scuttle_1.6.3               xml2_1.3.3                 
[121] foreach_1.5.2               vipor_0.4.5                
[123] bslib_0.4.0                 DirichletMultinomial_1.38.0
[125] rngtools_1.5.2              estimability_1.4.1         
[127] rvest_1.0.3                 CVXR_1.0-10                
[129] doRNG_1.8.2                 yulab.utils_0.0.5          
[131] digest_0.6.30               vegan_2.6-4                
[133] rmarkdown_2.17              cellranger_1.1.0           
[135] tidytree_0.4.1              htmlTable_2.4.1            
[137] gld_2.6.5                   DelayedMatrixStats_1.18.2  
[139] nloptr_2.0.3                lifecycle_1.0.3            
[141] nlme_3.1-160                jsonlite_1.8.2             
[143] BiocNeighbors_1.14.0        viridisLite_0.4.1          
[145] fansi_1.0.3                 pillar_1.8.1               
[147] lattice_0.20-45             fastmap_1.1.0              
[149] httr_1.4.4                  survival_3.4-0             
[151] glue_1.6.2                  png_0.1-7                  
[153] iterators_1.0.14            bit_4.0.4                  
[155] class_7.3-20                stringi_1.7.8              
[157] sass_0.4.2                  blob_1.2.3                 
[159] BiocSingular_1.12.0         latticeExtra_0.6-30        
[161] memoise_2.0.1               Rmpfr_0.8-9                
[163] irlba_2.3.5.1               e1071_1.7-11               
[165] 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.

Lin, Huang, and Shyamal Das Peddada. 2020. “Analysis of Compositions of Microbiomes with Bias Correction.” Nature Communications 11 (1): 1–11.

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.