diff --git a/R/deseq2.R b/R/deseq2.R
index 68661b6e..53912054 100644
--- a/R/deseq2.R
+++ b/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)
 }
diff --git a/R/differential.R b/R/differential.R
index 4e912415..22b95f08 100644
--- a/R/differential.R
+++ b/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)
 }
diff --git a/R/plot_volcano.R b/R/plot_volcano.R
index 116dae36..e09d20a8 100644
--- a/R/plot_volcano.R
+++ b/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, ...)
 # }
diff --git a/man/batch_correct_counts.Rd b/man/batch_correct_counts.Rd
index b1fa1db0..07728149 100644
--- a/man/batch_correct_counts.Rd
+++ b/man/batch_correct_counts.Rd
@@ -110,7 +110,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/clean_raw_counts.Rd b/man/clean_raw_counts.Rd
index e8705312..b585a349 100644
--- a/man/clean_raw_counts.Rd
+++ b/man/clean_raw_counts.Rd
@@ -99,7 +99,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/diff_counts.Rd b/man/diff_counts.Rd
index 3d2df5b6..29b4d825 100644
--- a/man/diff_counts.Rd
+++ b/man/diff_counts.Rd
@@ -6,6 +6,7 @@
 \usage{
 diff_counts(
   moo,
+  diff_method = "limma",
   count_type = "filt",
   sub_count_type = NULL,
   sample_id_colname = NULL,
@@ -25,6 +26,8 @@ diff_counts(
 \arguments{
 \item{moo}{multiOmicDataSet object (see \code{create_multiOmicDataSet_from_dataframes()})}
 
+\item{diff_method}{method for differential analysis (Default: "limma")}
+
 \item{count_type}{the type of counts to use -- must be a name in the counts slot (\code{moo@counts})}
 
 \item{sub_count_type}{if \code{count_type} is a list, specify the sub count type within the list. (Default: \code{NULL})}
@@ -64,7 +67,7 @@ contrast(s)}
 \item{plots_subdir}{subdirectory in where plots will be saved if \code{save_plots} is \code{TRUE}}
 }
 \value{
-\code{multiOmicDataSet} with \code{diff} added to the \code{analyses} slot (i.e. \code{moo@analyses$diff})
+\code{multiOmicDataSet} with \code{diff} added to the \code{analyses} slot (i.e. \code{moo@analyses$limma$diff})
 }
 \description{
 Differential expression analysis
@@ -89,7 +92,7 @@ moo <- multiOmicDataSet(
     contrasts = c("B-A", "C-A", "B-C"),
     voom_normalization_method = "quantile",
   )
-head(moo@analyses$diff)
+head(moo@analyses$limma$diff)
 }
 \seealso{
 Other moo methods: 
@@ -103,7 +106,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/filter_counts.Rd b/man/filter_counts.Rd
index 704ca61c..f9d94f27 100644
--- a/man/filter_counts.Rd
+++ b/man/filter_counts.Rd
@@ -184,7 +184,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/filter_diff.Rd b/man/filter_diff.Rd
index 5897762a..29440d07 100644
--- a/man/filter_diff.Rd
+++ b/man/filter_diff.Rd
@@ -6,6 +6,7 @@
 \usage{
 filter_diff(
   moo,
+  diff_method = "limma",
   feature_id_colname = NULL,
   significance_column = "adjpval",
   significance_cutoff = 0.05,
@@ -36,6 +37,8 @@ filter_diff(
 \arguments{
 \item{moo}{multiOmicDataSet object (see \code{create_multiOmicDataSet_from_dataframes()})}
 
+\item{diff_method}{where to access differential results from the analysis slot (default: "limma")}
+
 \item{feature_id_colname}{The column from the counts data containing the Feature IDs (Usually Gene or Protein ID).
 This is usually the first column of your input Counts Matrix. Only columns of Text type from your input Counts
 Matrix will be available to select for this parameter. (Default: \code{NULL} - first column in the counts matrix will be
@@ -124,7 +127,7 @@ moo <- multiOmicDataSet(
     voom_normalization_method = "quantile",
   ) \%>\%
   filter_diff()
-head(moo@analyses$diff_filt)
+head(moo@analyses[["limma"]][["diff_filt"]])
 }
 \seealso{
 Other moo methods: 
@@ -138,7 +141,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/normalize_counts.Rd b/man/normalize_counts.Rd
index f6667bf6..83baabe0 100644
--- a/man/normalize_counts.Rd
+++ b/man/normalize_counts.Rd
@@ -155,7 +155,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/plot_corr_heatmap.Rd b/man/plot_corr_heatmap.Rd
index d24840f3..bfe2c488 100644
--- a/man/plot_corr_heatmap.Rd
+++ b/man/plot_corr_heatmap.Rd
@@ -93,7 +93,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{heatmaps}
diff --git a/man/plot_expr_heatmap.Rd b/man/plot_expr_heatmap.Rd
index 5b7c3dea..b31bf897 100644
--- a/man/plot_expr_heatmap.Rd
+++ b/man/plot_expr_heatmap.Rd
@@ -154,7 +154,7 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{heatmaps}
diff --git a/man/plot_histogram.Rd b/man/plot_histogram.Rd
index 4bb95955..42386002 100644
--- a/man/plot_histogram.Rd
+++ b/man/plot_histogram.Rd
@@ -74,7 +74,7 @@ Other moo methods:
 \code{\link{plot_expr_heatmap}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/plot_pca.Rd b/man/plot_pca.Rd
index 823afae2..645ff887 100644
--- a/man/plot_pca.Rd
+++ b/man/plot_pca.Rd
@@ -86,7 +86,7 @@ Other moo methods:
 \code{\link{plot_expr_heatmap}()},
 \code{\link{plot_histogram}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{PCA functions}
diff --git a/man/plot_read_depth.Rd b/man/plot_read_depth.Rd
index 1a3b16a9..e7f06efd 100644
--- a/man/plot_read_depth.Rd
+++ b/man/plot_read_depth.Rd
@@ -63,7 +63,7 @@ Other moo methods:
 \code{\link{plot_expr_heatmap}()},
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
-\code{\link{run_deseq2}()},
+\code{\link{run_deseq}()},
 \code{\link{set_color_pal}()}
 }
 \concept{moo methods}
diff --git a/man/run_deseq2.Rd b/man/run_deseq.Rd
similarity index 94%
rename from man/run_deseq2.Rd
rename to man/run_deseq.Rd
index 630aaabe..df5c9037 100644
--- a/man/run_deseq2.Rd
+++ b/man/run_deseq.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/deseq2.R
-\name{run_deseq2}
-\alias{run_deseq2}
+\name{run_deseq}
+\alias{run_deseq}
 \title{Run DESeq2 on a multiOmicDataSet}
 \usage{
-run_deseq2(moo, design, ...)
+run_deseq(moo, design, ...)
 }
 \arguments{
 \item{moo}{multiOmicDataSet object}
diff --git a/man/set_color_pal.Rd b/man/set_color_pal.Rd
index 70c8ec57..8d750066 100644
--- a/man/set_color_pal.Rd
+++ b/man/set_color_pal.Rd
@@ -45,6 +45,6 @@ Other moo methods:
 \code{\link{plot_histogram}()},
 \code{\link{plot_pca}()},
 \code{\link{plot_read_depth}()},
-\code{\link{run_deseq2}()}
+\code{\link{run_deseq}()}
 }
 \concept{moo methods}
diff --git a/tests/testthat/test-differential.R b/tests/testthat/test-differential.R
index 6ef5c520..b1715b8e 100644
--- a/tests/testthat/test-differential.R
+++ b/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)
 })
diff --git a/tests/testthat/test-plot_volcano.R b/tests/testthat/test-plot_volcano.R
index c44bf935..58867faf 100644
--- a/tests/testthat/test-plot_volcano.R
+++ b/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)
diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd
index c1df7d62..524b9ffc 100644
--- a/vignettes/intro.Rmd
+++ b/vignettes/intro.Rmd
@@ -46,11 +46,11 @@ moo_nidap <- create_multiOmicDataSet_from_dataframes(
   ) %>%
   filter_diff()
 
-moo_nidap@analyses$diff %>%
+moo_nidap@analyses$limma$diff %>%
   join_dfs_wide() %>%
   head()
 
-moo_nidap@analyses$diff_filt %>% head()
+moo_nidap@analyses$limma$diff_filt %>% head()
 ```
 
 ## The multiOmicDataSet object structure
diff --git a/vignettes/visualization.Rmd b/vignettes/visualization.Rmd
index 92e484c1..ab66e0fa 100644
--- a/vignettes/visualization.Rmd
+++ b/vignettes/visualization.Rmd
@@ -119,7 +119,7 @@ print(heatmap_plot)
 #### Summary
 
 ```{r volcano_summary}
-dat_volcano_summary <- moo@analyses$diff %>%
+dat_volcano_summary <- moo@analyses$limma$diff %>%
   join_dfs_wide() %>%
   plot_volcano_summary()
 
@@ -129,7 +129,7 @@ head(dat_volcano_summary)
 #### Enhanced
 
 ```{r volcano_enhanced}
-dat_volcano_enhanced <- moo@analyses$diff %>%
+dat_volcano_enhanced <- moo@analyses$limma$diff %>%
   join_dfs_wide() %>%
   plot_volcano_enhanced()
 ```