edit dispersion for SE methods

stemangiola · Dec 7, 2023 · 0693ee1 · 0693ee1
1 parent 33c5dcd
commit 0693ee1
Show file tree

Hide file tree

Showing 5 changed files with 68 additions and 3 deletions.
diff --git a/R/functions_SE.R b/R/functions_SE.R
@@ -1133,10 +1133,17 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data,
     .data %>%
     assay(my_assay)
 
+  # Create design matrix for dispersion, removing random effects
+  design =
+    model.matrix(
+      object = .formula |> eliminate_random_effects(),
+      data = metadata
+    )
+
   if(quo_is_symbolic(.dispersion))
     dispersion = rowData(.data)[,quo_name(.dispersion),drop=FALSE] |> as_tibble(rownames = feature__$name) |> deframe()
   else
-    dispersion = setNames(edgeR::estimateDisp(counts)$tagwise.dispersion, rownames(counts))
+    dispersion = counts |> edgeR::estimateDisp(design = design) %$% tagwise.dispersion |> setNames(rownames(counts))
 
   # # Check dispersion
   # if(!names(dispersion) |> sort() |> identical(

diff --git a/R/methods.R b/R/methods.R
@@ -2600,6 +2600,7 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
 #' All methods use raw counts, irrespective of if scale_abundance or adjust_abundance have been calculated, therefore it is essential to add covariates such as batch effects (if applicable) in the formula.
 #'
 #' Underlying method for edgeR framework:
+#' 
 #' 	.data |>
 #'
 #' 	# Filter
@@ -2623,7 +2624,10 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
 #'			edgeR::glmQLFit(design) |> // or glmFit according to choice
 #'			edgeR::glmQLFTest(coef = 2, contrast = my_contrasts) // or glmLRT according to choice
 #'
+#'
+#'
 #'	Underlying method for DESeq2 framework:
+#'	
 #'	keep_abundant(
 #'			factor_of_interest = !!as.symbol(parse_formula(.formula)[[1]]),
 #'			minimum_counts = minimum_counts,
@@ -2636,6 +2640,31 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
 #'	DESeq2::results()
 #'
 #'
+#'
+#' Underlying method for glmmSeq framework:
+#'
+#' counts =
+#' .data %>%
+#'   assay(my_assay)
+#' 
+#' # Create design matrix for dispersion, removing random effects
+#' design =
+#'   model.matrix(
+#'     object = .formula |> eliminate_random_effects(),
+#'     data = metadata
+#'   )
+#' 
+#' dispersion = counts |> edgeR::estimateDisp(design = design) %$% tagwise.dispersion |> setNames(rownames(counts))
+#' 
+#'   glmmSeq( .formula,
+#'            countdata = counts ,
+#'            metadata =   metadata |> as.data.frame(),
+#'            dispersion = dispersion,
+#'            progress = TRUE,
+#'            method = method |> str_remove("(?i)^glmmSeq_" ),
+#'   )
+#'   
+#' 
 #' @return A consistent object (to the input) with additional columns for the statistics from the test (e.g.,  log fold change, p-value and false discovery rate).
 #'
 #'

diff --git a/man/test_differential_abundance-methods.Rd b/man/test_differential_abundance-methods.Rd
diff --git a/tests/testthat/test-bulk_methods.R b/tests/testthat/test-bulk_methods.R
@@ -855,7 +855,7 @@ test_that("differential trancript abundance - random effects",{
     pull(P_condition_adjusted) |>
     head(4) |>
     expect_equal(
-      c(0.02381167, 0.01097209, 0.01056741, 0.02381167),
+      c(0.03643414, 0.02938584, 0.02938584, 0.03643414),
       tolerance=1e-3
     )
 

diff --git a/tests/testthat/test-bulk_methods_SummarizedExperiment.R b/tests/testthat/test-bulk_methods_SummarizedExperiment.R
@@ -493,7 +493,7 @@ test_that("differential trancript abundance - random effects SE",{
  rowData(res)[,"P_condition_adjusted"] |>
     head(4) |>
     expect_equal(
-      c(0.1065866, 0.1109067, 0.1116562 , NA),
+      c(0.03394914, 0.03394914, 0.03394914,  NA),
       tolerance=1e-2
     )