One of the main characteristics of the Stammler_2016_16S_spikein dataset is the presence of spike-in bacteria with a known fixed amount of bacterial cells. These known loads of bacteria can be used to recalibrate the raw counts of the matrix and obtain recalibrated absolute counts. In this vignette, we provide an example of how to recalibrate the counts of the count matrix based on the read counts of Salinibacter ruber. This procedure is referred to as Spike-in-based calibration to total microbial load (SCML) in Sammler et al., 2016.

Import data

tse <- getBenchmarkData('Stammler_2016_16S_spikein', dryrun = FALSE)[[1]]
#> Warning: No taxonomy_tree available for Stammler_2016_16S_spikein.
#> Finished Stammler_2016_16S_spikein.
counts <- assay(tse)

Ids of the spike-in bacteria

Identifiers of the spiked-in bacteria have the suffix ‘XXXX’.

Bacteria ID Load
Salinibacter ruber AF323500XXXX 3.0 x 108
Rhizobium radiobacter AB247615XXXX 5.0 x 108
Alicyclobacillus acidiphilus AB076660XXXX 1.0 x 108

Recalibrate based on Salinibacter ruber abundance.

This recalibration is based on the original article. The only difference is that the numbers have been rounded up to obtain counts.

## AF323500XXXX is the unique OTU corresponding to S. ruber
s_ruber <- counts['AF323500XXXX', ]
size_factor <- s_ruber/mean(s_ruber)

SCML_data <- counts 
for(i in seq(ncol(SCML_data))){
    SCML_data[,i] <- round(SCML_data[,i] / size_factor[i])
}

Brief comparison of counts


no_cal <- counts |> 
    colSums() |> 
    as.data.frame() |> 
    tibble::rownames_to_column(var = 'sample_id') |> 
    magrittr::set_colnames(c('sample_id', 'colSum')) |> 
    mutate(calibrated = 'no') |> 
    as_tibble()

cal <-  SCML_data |> 
    colSums() |> 
    as.data.frame() |> 
    tibble::rownames_to_column(var = 'sample_id') |> 
    magrittr::set_colnames(c('sample_id', 'colSum')) |> 
    mutate(calibrated = 'yes') |> 
    as_tibble()

data <- bind_rows(no_cal, cal)

data |> 
    ggplot(aes(sample_id, colSum)) + 
    geom_col(aes(fill = calibrated), position = 'dodge') +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

The counts matrix can be replaced in the original tse in order to preserve the same metadata.

assay(tse) <- SCML_data
tse
#> class: TreeSummarizedExperiment 
#> dim: 4036 17 
#> metadata(0):
#> assays(1): counts
#> rownames(4036): GQ448052 EU458484 ... DQ795992 GQ492848
#> rowData names(1): taxonomy
#> colnames(17): MID26 MID27 ... MID42 MID43
#> colData names(12): dataset subject_id ... country description
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: NULL
#> rowTree: NULL
#> colLinks: NULL
#> colTree: NULL

Session information

sessionInfo()
#> R version 4.3.0 (2023-04-21)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] tidyr_1.3.0                    ggplot2_3.4.2                 
#>  [3] dplyr_1.1.2                    MicrobiomeBenchmarkData_1.3.0 
#>  [5] TreeSummarizedExperiment_2.8.0 Biostrings_2.68.0             
#>  [7] XVector_0.40.0                 SingleCellExperiment_1.22.0   
#>  [9] SummarizedExperiment_1.30.0    Biobase_2.60.0                
#> [11] GenomicRanges_1.52.0           GenomeInfoDb_1.36.0           
#> [13] IRanges_2.34.0                 S4Vectors_0.38.0              
#> [15] BiocGenerics_0.46.0            MatrixGenerics_1.12.0         
#> [17] matrixStats_0.63.0             BiocStyle_2.28.0              
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.0        farver_2.1.1            blob_1.2.4             
#>  [4] filelock_1.0.2          bitops_1.0-7            fastmap_1.1.1          
#>  [7] RCurl_1.98-1.12         lazyeval_0.2.2          BiocFileCache_2.8.0    
#> [10] digest_0.6.31           lifecycle_1.0.3         tidytree_0.4.2         
#> [13] RSQLite_2.3.1           magrittr_2.0.3          compiler_4.3.0         
#> [16] rlang_1.1.0             sass_0.4.5              tools_4.3.0            
#> [19] utf8_1.2.3              yaml_2.3.7              knitr_1.42             
#> [22] labeling_0.4.2          bit_4.0.5               curl_5.0.0             
#> [25] DelayedArray_0.25.0     BiocParallel_1.34.0     withr_2.5.0            
#> [28] purrr_1.0.1             desc_1.4.2              grid_4.3.0             
#> [31] fansi_1.0.4             colorspace_2.1-0        scales_1.2.1           
#> [34] cli_3.6.1               rmarkdown_2.21          crayon_1.5.2           
#> [37] ragg_1.2.5              treeio_1.24.0           generics_0.1.3         
#> [40] httr_1.4.5              DBI_1.1.3               ape_5.7-1              
#> [43] cachem_1.0.7            stringr_1.5.0           zlibbioc_1.46.0        
#> [46] parallel_4.3.0          BiocManager_1.30.20     vctrs_0.6.2            
#> [49] yulab.utils_0.0.6       Matrix_1.5-4            jsonlite_1.8.4         
#> [52] bookdown_0.33           bit64_4.0.5             systemfonts_1.0.4      
#> [55] jquerylib_0.1.4         glue_1.6.2              pkgdown_2.0.7          
#> [58] codetools_0.2-19        stringi_1.7.12          gtable_0.3.3           
#> [61] munsell_0.5.0           tibble_3.2.1            pillar_1.9.0           
#> [64] htmltools_0.5.5         GenomeInfoDbData_1.2.10 R6_2.5.1               
#> [67] dbplyr_2.3.2            textshaping_0.3.6       rprojroot_2.0.3        
#> [70] evaluate_0.20           lattice_0.21-8          highr_0.10             
#> [73] memoise_2.0.1           bslib_0.4.2             Rcpp_1.0.10            
#> [76] nlme_3.1-162            xfun_0.39               fs_1.6.2               
#> [79] pkgconfig_2.0.3