Add full SummarizedExperiment support

.github/workflows/R-CMD-check.yaml

@@ -64,7 +64,7 @@ jobs:
         run: |
           remotes::install_deps(dependencies = TRUE)
           remotes::install_cran("rcmdcheck")
-          if (!requireNamespace("BiocManager", quietly = TRUE)){install.packages("BiocManager")}; BiocManager::install(c("phyloseq", "limma"), ask = FALSE)
+          if (!requireNamespace("BiocManager", quietly = TRUE)){install.packages("BiocManager")}; BiocManager::install(c("phyloseq", "limma", "SummarizedExperiment"), ask = FALSE)
         shell: Rscript {0}
 
       - name: Check

.github/workflows/test-coverage.yaml

@@ -37,7 +37,7 @@ jobs:
       - name: Install dependencies
         run: |
           install.packages(c("remotes"))
-          if (!requireNamespace("BiocManager", quietly = TRUE)){install.packages("BiocManager")}; BiocManager::install(c("phyloseq", "limma"), ask = FALSE)
+          if (!requireNamespace("BiocManager", quietly = TRUE)){install.packages("BiocManager")}; BiocManager::install(c("phyloseq", "limma", "SummarizedExperiment"), ask = FALSE)
           remotes::install_deps(dependencies = TRUE)
           remotes::install_cran("covr")
         shell: Rscript {0}

DESCRIPTION

@@ -17,7 +17,7 @@ License: GPL (>= 2)
 Imports: stats, utils, VGAM, numDeriv, ggplot2, trust, dplyr, magrittr, detectseparation, scales, rlang
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: knitr,
     rmarkdown,
     testthat,
@@ -26,5 +26,6 @@ Suggests: knitr,
     slam,
     R.rsp,
     optimx,
-    phyloseq
+    phyloseq,
+    SummarizedExperiment
 VignetteBuilder: knitr

NAMESPACE

@@ -12,8 +12,10 @@ export(HDIbetabinom)
 export(bbdml)
 export(checkNested)
 export(clean_taxa_names)
+export(clean_taxa_names_se)
 export(contrastsTest)
 export(convert_phylo)
+export(convert_sumexp)
 export(coth)
 export(differentialTest)
 export(fishZ)
@@ -24,6 +26,7 @@ export(invlogit)
 export(logit)
 export(lrtest)
 export(otu_to_taxonomy)
+export(otu_to_taxonomy_se)
 export(pbLRT)
 export(pbRao)
 export(pbWald)
@@ -35,6 +38,7 @@ export(score)
 export(waldchisq)
 export(waldt)
 export(warn_phyloseq)
+export(warn_sumexp)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
 importFrom(stats,simulate)

R/bbdml.R

@@ -2,7 +2,7 @@
 #'
 #' @param formula an object of class \code{formula}: a symbolic description of the model to be fitted to the abundance
 #' @param phi.formula an object of class \code{formula} without the response: a symbolic description of the model to be fitted to the dispersion
-#' @param data a data frame or \code{phyloseq} object containing the variables in the models
+#' @param data a data frame, \code{phyloseq}, or \code{SummarizedExperiment} object containing the variables in the models
 #' @param link link function for abundance covariates, defaults to \code{"logit"}
 #' @param phi.link link function for dispersion covariates, defaults to \code{"logit"}
 #' @param method optimization method, defaults to \code{"trust"}, or see \code{\link[optimx]{optimr}} for other options
@@ -60,6 +60,14 @@ bbdml <- function(formula, phi.formula, data,
     formula <- stats::update(formula, cbind(W, M - W) ~ .)
   }
 
+  # Convert SummarizedExperiment objects
+  if (inherits(data, "SummarizedExperiment")) {
+    selection <- all.vars(formula)[1]
+    data <- convert_sumexp(data, select = selection)
+    # Update formula to match convert_phylo specification
+    formula <- stats::update(formula, cbind(W, M - W) ~ .)
+  }
+
   # Record call
   call <- match.call(expand.dots = FALSE)
   # Record mu link

R/clean_taxa_names.R

@@ -23,3 +23,30 @@ clean_taxa_names <- function(x, name = "OTU") {
     warn_phyloseq()
   }
 }
+
+#' Rename taxa
+#'
+#' Renames taxa to have short human-readable names
+#'
+#' @param x Object of class \code{SummarizedExperiment}
+#' @param name Character, defaults to \code{"OTU"}. Optional. String to use in every taxa name.
+#'
+#' @details The original taxa names are saved as the \code{original_names} attribute. See the example for an example of how to access the original names.
+#'
+#' @return Object of class \code{SummarizedExperiment}, with taxa renamed (defaults to OTU1, OTU2, ...), with the original taxa names saved as an attribute.
+#'
+#' @export
+clean_taxa_names_se <- function(x, name = "OTU") {
+  if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
+    if (inherits(x, "SummarizedExperiment")) {
+      attr(x, "original_names") <- row.names(x)
+      row.names(x) <- paste0(name, seq_len(nrow(x)))
+      return(x)
+    } else {
+      stop("clean_taxa_names_se is intended for SummarizedExperiment objects!")
+    }
+  } else {
+    warn_sumexp()
+  }
+
+}

R/contrastsTest.R

@@ -95,6 +95,14 @@ contrastsTest <- function(formula, phi.formula,
     } else {
       warn_phyloseq()
     }
+  } else if (inherits(data, "SummarizedExperiment")) {
+    if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
+      # Set up response
+      taxanames <- row.names(data)
+      sample_data <- SummarizedExperiment::colData(data)
+    } else {
+      warn_sumexp()
+    }
   } else if (is.matrix(data) || is.data.frame(data)) {
 
     # # use phyloseq
@@ -118,7 +126,7 @@ contrastsTest <- function(formula, phi.formula,
     M <- rowSums(data)
 
   } else {
-    stop("Input must be either data frame, matrix, or phyloseq object!")
+    stop("Input must be either data frame, matrix, phyloseq object, or SummarizedExperiment object!")
   }
 
   # Set up output
@@ -147,6 +155,8 @@ contrastsTest <- function(formula, phi.formula,
     # Subset data to only select that taxa
     if ("phyloseq" %in% class(data)) {
       data_i <- convert_phylo(data, select = taxanames[i])
+    } else if (inherits(data, "SummarizedExperiment")) {
+      data_i <- convert_sumexp(data, select = taxanames[i])
     } else {
       response_i <- data.frame(W = data[, taxanames[i]], M = M)
       data_i <- cbind(response_i, sample_data)

R/convert_sumexp.R

@@ -0,0 +1,22 @@
+#' Function to subset and convert SummarizedExperiment data
+#'
+#' @param data a \code{SummarizedExperiment} object
+#' @param select Name of OTU or taxa to select, must match taxa name in \code{data}
+#'
+#' @return A \code{data.frame} object, with elements \code{W} as the observed counts, \code{M} as the sequencing depth, and the sample data with their original names.
+#'
+#' @export
+convert_sumexp <- function(data, select) {
+  if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
+    subsamp <- data[select, ]
+    W_tmp <- matrix(t(SummarizedExperiment::assay(subsamp)), ncol = 1)
+
+
+    out <- data.frame(W = W_tmp,
+                      M = colSums(SummarizedExperiment::assay(data)),
+                      SummarizedExperiment::colData(subsamp))
+    return(out)
+  } else {
+    warn_sumexp()
+  }
+}

R/differentialTest.R

@@ -4,7 +4,7 @@
 #' @param phi.formula an object of class \code{formula} without the response: a symbolic description of the model to be fitted to the dispersion
 #' @param formula_null Formula for mean under null, without response
 #' @param phi.formula_null Formula for overdispersion under null, without response
-#' @param data a data frame containing the OTU table, or \code{phyloseq} object containing the variables in the models
+#' @param data a data frame containing the OTU table, \code{phyloseq}, or \code{SummarizedExperiment} object containing the variables in the models
 #' @param link link function for abundance covariates, defaults to \code{"logit"}
 #' @param phi.link link function for dispersion covariates, defaults to \code{"logit"}
 #' @param test Character. Hypothesis testing procedure to use. One of \code{"Wald"}, \code{"LRT"} (likelihood ratio test), or \code{"Rao"}.
@@ -95,6 +95,15 @@ differentialTest <- function(formula, phi.formula,
     } else {
       warn_phyloseq()
     }
+  } else if (inherits(data, "SummarizedExperiment")) {
+    if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
+      # Set up response
+      taxanames <- row.names(data)
+      sample_data <- SummarizedExperiment::colData(data)
+    } else {
+      warn_sumexp()
+    }
+
   } else if (is.matrix(data) || is.data.frame(data)) {
 
     # # use phyloseq
@@ -118,7 +127,7 @@ differentialTest <- function(formula, phi.formula,
     M <- rowSums(data)
 
   } else {
-    stop("Input must be either data frame, matrix, or phyloseq object!")
+    stop("Input must be either data frame, matrix, phyloseq object or SummarizedExperiment!")
   }
 
   # Set up output
@@ -158,6 +167,8 @@ differentialTest <- function(formula, phi.formula,
     # Subset data to only select that taxa
     if ("phyloseq" %in% class(data)) {
       data_i <- convert_phylo(data, select = taxanames[i])
+    } else if (inherits(data, "SummarizedExperiment")) {
+      data_i <- convert_sumexp(data, select = taxanames[i])
     } else {
       response_i <- data.frame(W = data[, taxanames[i]], M = M)
       data_i <- cbind(response_i, sample_data)
@@ -275,6 +286,8 @@ We *strongly recommend* running `bbdml` on a single taxon (especially before pos
     i <- (try_only[!(try_only %in% ind_disc)])[1]
     if ("phyloseq" %in% class(data)) {
       data_i <- convert_phylo(data, select = taxanames[i])
+    } else if (inherits(data, "SummarizedExperiment")) {
+      data_i <- convert_sumexp(data, select = taxanames[i])
     } else {
       response_i <- data.frame(W = data[, taxanames[i]], M = M)
       data_i <- cbind(response_i, sample_data)

R/genInits.R

@@ -33,7 +33,7 @@ genInits <- function(W, M,
                      nstart = 1, use = TRUE) {
 
 
-  init.glm <- eval(parse(text = paste("quasibinomial(link =", link,")")))
+  init.glm <- eval(parse(text = paste("quasibinomial(link =", link, ")")))
   tmp <- stats::glm.fit(x = X, y = cbind(W, M - W), family = init.glm)
   b_init <- stats::coef(tmp)
   # Just use 0.5 for phi_init

R/getRestrictionTerms.R

@@ -18,9 +18,9 @@ getRestrictionTerms <- function(mod, mod_null = NULL, restrictions = NULL, restr
   assigns.mu <- attr(mod$X.mu, "assign")
   assigns.phi <- attr(mod$X.phi, "assign")
 
-  sortInteraction <- function(x){
+  sortInteraction <- function(x) {
     # Sorts interaction effects so they always match
-    sapply(lapply(strsplit(x,":"), sort), paste, collapse = ":")
+    sapply(lapply(strsplit(x, ":"), sort), paste, collapse = ":")
   }
 
   if (!is.null(mod_null)) {

R/otu_to_taxonomy.R

@@ -24,3 +24,31 @@ otu_to_taxonomy <- function(OTU, data, level = NULL) {
     warn_phyloseq()
   }
 }
+
+
+#' Transform OTUs to their taxonomic label
+#'
+#' @param OTU String vector. Names of OTU labels in \code{data}
+#' @param data \code{phyloseq} object with a taxonomy table
+#' @param level (Optional). Character vector. Desired taxonomic levels for output.
+#'
+#' @importFrom magrittr %>%
+#'
+#' @return String vector. Names of taxonomic labels matching labels of \code{OTU}.
+#'
+#' @export
+otu_to_taxonomy_se <- function(OTU, data, level = NULL) {
+  if (!inherits(data, "SummarizedExperiment")) {
+    stop("This function currently only works for SummarizedExperiment objects.")
+  }
+
+  if (!requireNamespace("SummarizedExperiment", quietly = TRUE)) {
+    warn_sumexp()
+  }
+
+  if (is.null(level)) {
+    apply(SummarizedExperiment::rowData(data)[OTU,], 1, function(x) {paste(stats::na.omit(x), collapse = "_")})
+  } else {
+    sapply(SummarizedExperiment::rowData(data)[OTU, level], function(x) {paste(stats::na.omit(x), collapse = "_")})
+  }
+}

R/plot_differentialTest.R

@@ -35,6 +35,17 @@ plot.differentialTest <- function(x, level = NULL, data_only = FALSE, ...) {
       }
     }
   }
+
+  if (inherits(x$data, "SummarizedExperiment")) {
+    if (!(nrow(SummarizedExperiment::rowData(x$data)) == 0) && "kingdom" %in% tolower(colnames(SummarizedExperiment::rowData(x$data)))) {
+      signif_taxa <- otu_to_taxonomy_se(signif_taxa, x$data, level = level)
+      if (length(unique(signif_taxa)) != length(unique(x$significant_taxa))) {
+        # Make sure if repeated taxa add unique otu identifiers
+        signif_taxa <- paste0(signif_taxa, " (", x$significant_taxa, ")")
+      }
+    }
+  }
+
   if (length(x$significant_models) != 0) {
     var_per_mod <- length(x$restrictions_DA) + length(x$restrictions_DV)
     total_var_count <- length(signif_taxa) * var_per_mod

R/warn_phyloseq.R

@@ -4,3 +4,10 @@
 warn_phyloseq <- function() {
   stop("You are trying to use a `phyloseq` data object or `phyloseq` helper function without having the `phyloseq` package installed. Please either install the package or use a standard data frame.")
 }
+
+#' Function to throw error if the `SummarizedExperiment` package is called but it is not installed
+#'
+#' @export
+warn_sumexp <- function() {
+  stop("You are trying to use a `SummarizedExperiment` data object or `SummarizedExperiment` helper function without having the `SummarizedExperiment` package installed. Please either install the package or use a standard data frame.")
+}