Merge pull request #190 from luhann/tse

Add full SummarizedExperiment support

Author: Sarah Teichman

.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