Merge pull request #12 from waldronlab/sdgamboa/scml

Add scml function for the dataset with spike-in bacteria

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MicrobiomeBenchmarkData
 Title: Datasets for benchmarking in microbiome research
-Version: 1.7.0
+Version: 1.7.1
 Description: The MicrobiomeBenchmarkData package provides functionality to 
     access microbiome datasets suitable for benchmarking. These datasets
     have some biological truth, which allows to have expected results for
@@ -58,6 +58,6 @@ URL: https://github.com/waldronlab/MicrobiomeBenchmarkData,
 BiocType: ExperimentData
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 Config/testthat/edition: 3

NAMESPACE

@@ -2,6 +2,7 @@
 
 export(getBenchmarkData)
 export(removeCache)
+export(scml)
 import(SummarizedExperiment)
 importFrom(BiocFileCache,BiocFileCache)
 importFrom(BiocFileCache,bfcadd)

NEWS.md

@@ -1,3 +1,8 @@
 # Changes in version 0.99.0 (2021-09-26)
 
 + Added a `NEWS.md` file to track changes to the package.
+
+# Changes in version  1.7.1 (2021-09-26)
+
++ Added a function (`scml`) for re-calibrating the "Stammler_2016_16S_spikein"
+dataset with SCML: spike-in-based calibration to total microbial load.

R/scml.R

@@ -0,0 +1,50 @@
+#' SCML: spike-in-based calibration to total microbial load
+#' 
+#' The \code{scml} function applies the
+#' spike-in-based calibration to total microbial load (SCML) method to
+#'
+#' @param tse A treeSummarizedExperiment from the \code{getBenchmarkData}
+#' function.
+#' @param bac A character. One of the following options:
+#' s = Salinibacter ruber (AF323500), r = Rhizobium radiobacter (AB247615),
+#' a, = Alicyclobacillus acidiphilus (AB076660)
+#' 
+#' @return A TreeSummarizedExperiment with SCML data instead of counts.
+#' @export
+#'
+#' @examples
+#' tse <- getBenchmarkData("Stammler_2016_16S_spikein", dryrun = FALSE)[[1]]
+#' tseSCML <- scml(tse, bac = "s")
+#' 
+scml <- function(tse, bac = c("s", "r", "a")) {
+    bacLetter <- match.arg(bac)
+    bacNames <- c(
+        s = "AF323500XXXX", r = "AB247615XXXX", a = "AB076660XXXX"
+    )
+    bacFullNames <- c(
+        AF323500XXXX = "Salinibacter ruber (AF323500)",
+        AB247615XXXX = "Rhizobium radiobacter (AB247615)",
+        AB076660XXXX = "Alicyclobacillus acidiphilus (AB076660)"
+
+    )
+    bacName <- bacNames[bacLetter]
+
+    lgl <- bacName %in% rownames(tse)
+    if (!lgl) {
+        stop(
+            "Feature", bacName, "not found.",
+            "Are you sure you're using the Stammler_2016_16S_spikein",
+            call. = FALSE
+        )
+    }
+    message("Re-calibrating counts with ", bacFullNames[bacName])
+    counts <- SummarizedExperiment::assay(tse, 1)
+    bacAb <- counts[bacName, ]
+    sizeFactor <- bacAb/mean(bacAb)
+    scmlData <- counts 
+    for(i in seq(ncol(scmlData))){
+        scmlData[,i] <- round(scmlData[,i] / sizeFactor[i])
+    }
+    SummarizedExperiment::assay(tse) <- scmlData
+    return(tse)
+}

tests/testthat/test-smcl.R

@@ -0,0 +1,6 @@
+test_that("scml works", {
+    tse <- suppressWarnings(getBenchmarkData("Stammler_2016_16S_spikein", dryrun = FALSE)[[1]])
+    expect_s4_class(scml(tse, bac = "s"), "TreeSummarizedExperiment")
+    expect_s4_class(scml(tse, bac = "r"), "TreeSummarizedExperiment")
+    expect_s4_class(scml(tse, bac = "a"), "TreeSummarizedExperiment")
+})

vignettes/recalibrare_spikein_data.Rmd

@@ -108,6 +108,13 @@ assay(tse) <- SCML_data
 tse
 ```
 
+## A more convenient way using the scml function included in the package:
+
+```{r}
+tse <- getBenchmarkData('Stammler_2016_16S_spikein', dryrun = FALSE)[[1]]
+tse <- scml(tse,bac = "s")
+```
+
 # Session information
 
 ```{r}