feat: support DESeq2

R/deseq2.R

@@ -21,22 +21,58 @@
 #' moo <- run_deseq2(moo, ~condition)
 #' }
 #' @family moo methods
-run_deseq2 <- S7::new_generic("run_deseq2", "moo", function(moo, design, ...) {
+run_deseq <- S7::new_generic("run_deseq2", "moo", function(moo, design, ...) {
   return(S7::S7_dispatch())
 })
 
-S7::method(run_deseq2, multiOmicDataSet) <- function(moo, design, feature_id_colname = "gene_id", min_count = 10, ...) {
-  if (is.null(moo@counts$filt)) {
-    stop("moo must contain filtered counts for DESeq2. Hint: Did you forget to run filter_counts()?")
+S7::method(run_deseq, multiOmicDataSet) <- function(moo,
+                                                    design,
+                                                    count_type = "filt",
+                                                    sub_count_type = NULL,
+                                                    feature_id_colname = NULL,
+                                                    min_count = 10, ...) {
+  sample_metadata <- moo@sample_meta
+  message(glue::glue("* differential counts with DESeq2"))
+  # select correct counts matrix
+  if (!(count_type %in% names(moo@counts))) {
+    stop(glue::glue("count_type {count_type} not in moo@counts"))
   }
+  counts_dat <- moo@counts[[count_type]]
+  if (!is.null(sub_count_type)) {
+    if (!(inherits(counts_dat, "list"))) {
+      stop(
+        glue::glue(
+          "{count_type} counts is not a named list. To use {count_type} counts, set sub_count_type to NULL"
+        )
+      )
+    } else if (!(sub_count_type %in% names(counts_dat))) {
+      stop(
+        glue::glue(
+          "sub_count_type {sub_count_type} is not in moo@counts[[{count_type}]]"
+        )
+      )
+    }
+    counts_dat <- moo@counts[[count_type]][[sub_count_type]]
+  }
+  if (is.null(sample_id_colname)) {
+    sample_id_colname <- colnames(sample_metadata)[1]
+  }
+  if (is.null(feature_id_colname)) {
+    feature_id_colname <- colnames(counts_dat)[1]
+  }
+
   dds <- DESeq2::DESeqDataSetFromMatrix(
-    countData = moo@counts$filt %>%
+    countData = counts_dat %>%
       dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round)) %>% # DESeq2 requires integer counts
       counts_dat_to_matrix(feature_id_colname = feature_id_colname),
-    colData = moo@sample_meta,
+    colData = sample_metadata,
     design = design
-  )
-  moo@analyses$deseq2_ds <- DESeq2::DESeq(dds, ...)
-  moo@analyses$deseq2_results <- DESeq2::results(moo@analyses$deseq2_ds)
+  ) %>%
+    DESeq2::DESeq(...)
+
+  moo@analyses$deseq$dds <- dds
+  moo@analyses$deseq$colData <- DESeq2::colData(dds)
+  moo@analyses$deseq$res <- DESeq2::results(dds)
+  moo@analyses$deseq$vsd <- DESeq2::vst(dds, blind = FALSE)
   return(moo)
 }

R/differential.R

@@ -5,6 +5,7 @@
 #' @inheritParams normalize_counts
 #' @inheritParams option_params
 #'
+#' @param diff_method method for differential analysis (Default: "limma")
 #' @param sub_count_type if `count_type` is a list, specify the sub count type within the list. (Default: `NULL`)
 #' @param covariates_colnames The column name(s) from the sample metadata containing variable(s) of interest, such as
 #'   phenotype. Most commonly this will be the same column selected for your Groups Column. Some experimental designs
@@ -17,7 +18,7 @@
 #' @param return_mean_and_sd if TRUE, return Mean and Standard Deviation of groups in addition to DEG estimates for
 #'   contrast(s)
 #'
-#' @returns `multiOmicDataSet` with `diff` added to the `analyses` slot (i.e. `moo@analyses$diff`)
+#' @returns `multiOmicDataSet` with `diff` added to the `analyses` slot (i.e. `moo@analyses$limma$diff`)
 #' @export
 #'
 #' @family moo methods
@@ -42,8 +43,9 @@
 #'     contrasts = c("B-A", "C-A", "B-C"),
 #'     voom_normalization_method = "quantile",
 #'   )
-#' head(moo@analyses$diff)
+#' head(moo@analyses$limma$diff)
 diff_counts <- function(moo,
+                        diff_method = "limma",
                         count_type = "filt",
                         sub_count_type = NULL,
                         sample_id_colname = NULL,
@@ -396,6 +398,7 @@ diff_counts <- function(moo,
 
   # Mean-variance Plot.
   mv_plot <- plot_mean_variance(voom_elist = v)
+  plots_subdir <- file.path(diff_method, plots_subdir)
   print_or_save_plot(
     mv_plot,
     filename = file.path(plots_subdir, "mean-variance.png"),
@@ -423,7 +426,7 @@ diff_counts <- function(moo,
 
   names(df_list) <- contrasts
 
-  moo@analyses[["diff"]] <- df_list
+  moo@analyses[[diff_method]][["diff"]] <- df_list
   return(moo)
 }
 
@@ -509,6 +512,7 @@ plot_mean_variance <- function(voom_elist) {
 #'
 #' @inheritParams option_params
 #' @inheritParams filter_counts
+#' @param diff_method where to access differential results from the analysis slot (default: "limma")
 #' @param significance_column Column name for significance, e.g. `"pval"` or `"pvaladj"` (default)
 #' @param significance_cutoff Features will only be kept if their `significance_column` is less then this cutoff
 #'   threshold
@@ -567,8 +571,9 @@ plot_mean_variance <- function(voom_elist) {
 #'     voom_normalization_method = "quantile",
 #'   ) %>%
 #'   filter_diff()
-#' head(moo@analyses$diff_filt)
+#' head(moo@analyses[["limma"]][["diff_filt"]])
 filter_diff <- function(moo,
+                        diff_method = "limma",
                         feature_id_colname = NULL,
                         significance_column = "adjpval",
                         significance_cutoff = 0.05,
@@ -597,7 +602,7 @@ filter_diff <- function(moo,
   Count <- Count_format <- L1 <- Label <- Percent <- Significant <- Var1 <- Var2 <- value <- NULL
 
   # from NIDAP DEG_Gene_List template - filters DEG table
-  diff_dat <- moo@analyses$diff %>%
+  diff_dat <- moo@analyses[[diff_method]]$diff %>%
     join_dfs_wide() %>%
     as.data.frame()
 
@@ -810,6 +815,8 @@ filter_diff <- function(moo,
       ) +
       ggplot2::xlab("") +
       ggplot2::scale_y_continuous(name = "", expand = c(y_axis_expansion, 0))
+
+    plots_subdir <- file.path(diff_method, plots_subdir)
     print_or_save_plot(
       pp,
       filename = file.path(plots_subdir, glue::glue("{plot_type}chart.png")),
@@ -977,6 +984,6 @@ filter_diff <- function(moo,
     }
     ### PH: END Create DEG summary PieChart
   }
-  moo@analyses$diff_filt <- out
+  moo@analyses[[diff_method]][["diff_filt"]] <- out
   return(moo)
 }

R/plot_volcano.R

@@ -2,13 +2,13 @@
 #   S7::S7_dispatch()
 # })
 #
-# S7::method(plot_volcano_summary, multiOmicDataSet) <- function(moo_diff,
+# S7::method(plot_volcano_summary, multiOmicDataSet) <- function(moo_diff, diff_method = "limma",
 #                                                                ...) {
-#   diff_dat <- moo_diff@analyses$diff %>% join_dfs_wide()
+#   diff_dat <- moo_diff@analyses[[diff_method]][["diff"]] %>% join_dfs_wide()
 #   plot_volcano_summary(diff_dat, ...)
 # }
 
 # S7::method(plot_volcano_enhanced, multiOmicDataSet) <- function(moo_diff, ...) {
-#   diff_dat <- moo_diff@analyses$diff %>% join_dfs_wide()
+#   diff_dat <- moo_diff@analyses[[diff_method]][["diff"]] %>% join_dfs_wide()
 #   plot_volcano_enhanced(diff_dat, ...)
 # }

tests/testthat/test-differential.R

@@ -26,13 +26,13 @@ test_that("differential analysis works for NIDAP", {
   #   dplyr::arrange(Gene) %>%
   #   dplyr::select(order(colnames(.))) %>% dplyr::select(-C2) %>%
   #   as.data.frame()
-  # y <- moo@analyses$diff %>%
+  # y <- moo@analyses$limma$diff %>%
   #   dplyr::arrange(Gene) %>%
   #   dplyr::select(order(colnames(.)))
   # equal_dfs(x, y)
 
   expect_equal(
-    deg_moo@analyses$diff %>%
+    deg_moo@analyses$limma$diff %>%
       join_dfs_wide() %>%
       dplyr::arrange(Gene) %>%
       dplyr::select(order(colnames(.))),
@@ -71,7 +71,7 @@ test_that("diff_counts works for RENEE", {
       return_mean_and_sd = TRUE,
       input_in_log_counts = TRUE
     )
-  actual <- moo_renee@analyses$diff[[1]] %>%
+  actual <- moo_renee@analyses$limma$diff[[1]] %>%
     head() %>%
     dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, digits = 0)))
   expected <- tibble::tibble(
@@ -211,5 +211,5 @@ test_that("filter_diff works for NIDAP", {
       plot_type = "bar",
       plot_titles_fontsize = 12
     )
-  expect_equal(moo@analyses$diff_filt, nidap_deg_gene_list)
+  expect_equal(moo@analyses$limma$diff_filt, nidap_deg_gene_list)
 })

tests/testthat/test-plot_volcano.R

@@ -8,7 +8,7 @@
 #     "norm" = list("voom" = as.data.frame(nidap_norm_counts))
 #   )
 # )
-# moo_nidap@analyses$diff <- nidap_deg_analysis_2
+# moo_nidap@analyses$limma$diff <- nidap_deg_analysis_2
 #
 # test_that("volcano plots work on MOO", {
 #   volc_sum <- plot_volcano_summary(moo_nidap)