Merge pull request #288 from stemangiola/custom-dispersion-DE

Custom dispersion de
stemangiola · Aug 29, 2023 · fd8b3aa · fd8b3aa
2 parents c83614f + cb96f32
commit fd8b3aa
Show file tree

Hide file tree

Showing 6 changed files with 178 additions and 92 deletions.
diff --git a/R/functions.R b/R/functions.R
@@ -685,13 +685,15 @@ get_differential_transcript_abundance_glmmSeq <- function(.data,
                                                                omit_contrast_in_colnames = FALSE,
                                                                prefix = "",
                                                                .sample_total_read_count = NULL,
+                                                          .dispersion = NULL,
                                                           ...) {
   # Get column names
   .sample = enquo(.sample)
   .transcript = enquo(.transcript)
   .abundance = enquo(.abundance)
   .sample_total_read_count = enquo(.sample_total_read_count)
-
+  .dispersion = enquo(.dispersion)
+
   # Check if omit_contrast_in_colnames is correctly setup
   if(omit_contrast_in_colnames & length(.contrasts) > 1){
     warning("tidybulk says: you can omit contrasts in column names only when maximum one contrast is present")
@@ -755,13 +757,27 @@ get_differential_transcript_abundance_glmmSeq <- function(.data,
 
   # Reorder counts
   counts = counts[,rownames(metadata),drop=FALSE]
-
-  glmmSeq_object =
+
+  if(quo_is_symbolic(.dispersion))
+    dispersion = .data |> pivot_transcript(!!.transcript) |> select(!!.transcript, !!.dispersion) |> deframe()
+  else
+    dispersion = setNames(edgeR::estimateDisp(counts)$tagwise.dispersion, rownames(counts))
+
+  # # Check dispersion
+  # if(!names(dispersion) |> sort() |> identical(
+  #   rownames(counts) |>
+  #   sort()
+  # )) stop("tidybulk says: The features in the dispersion vector do not overlap with the feature in the assay")
+
+  # Make sure the order matches the counts
+  dispersion = dispersion[rownames(counts)]
+
+  glmmSeq_object = 
     glmmSeq( .formula,
           countdata = counts ,
           metadata =   metadata,
-          dispersion = setNames(edgeR::estimateDisp(counts)$tagwise.dispersion, rownames(counts)),
-          progress = TRUE,
+          dispersion = dispersion,
+          progress = TRUE, 
           method = method |> str_remove("(?i)^glmmSeq_"),
           ...
   )