Merge pull request #306 from stemangiola/improve-glmmseq

Improve glmmseq

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: tidybulk
 Title: Brings transcriptomics to the tidyverse 
-Version: 1.15.3
+Version: 1.15.4
 Authors@R: c(person("Stefano", "Mangiola", email = "[email protected]",
                   role = c("aut", "cre")),
             person("Maria", "Doyle", email = "[email protected]",

R/functions.R

@@ -779,11 +779,16 @@ get_differential_transcript_abundance_glmmSeq <- function(.data,
   # Make sure the order matches the counts
   dispersion = dispersion[rownames(counts)]
 
+  # Scaling
+  sizeFactors <- counts |> edgeR::calcNormFactors(method = scaling_method)
+
+
   glmmSeq_object =
     glmmSeq( .formula,
           countdata = counts ,
           metadata =   metadata,
           dispersion = dispersion,
+          sizeFactors = sizeFactors,
           progress = TRUE,
           method = method |> str_remove("(?i)^glmmSeq_"),
           ...

R/functions_SE.R

@@ -1100,13 +1100,13 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data,
     omit_contrast_in_colnames = FALSE
   }
 
-  # Check if package is installed, otherwise install
-  if (find.package("edgeR", quiet = TRUE) %>% length %>% equals(0)) {
-    message("tidybulk says: Installing edgeR needed for differential transcript abundance analyses")
-    if (!requireNamespace("BiocManager", quietly = TRUE))
-      install.packages("BiocManager", repos = "https://cloud.r-project.org")
-    BiocManager::install("edgeR", ask = FALSE)
-  }
+  # # Check if package is installed, otherwise install
+  # if (find.package("edgeR", quiet = TRUE) %>% length %>% equals(0)) {
+  #   message("tidybulk says: Installing edgeR needed for differential transcript abundance analyses")
+  #   if (!requireNamespace("BiocManager", quietly = TRUE))
+  #     install.packages("BiocManager", repos = "https://cloud.r-project.org")
+  #   BiocManager::install("edgeR", ask = FALSE)
+  # }
 
   # Check if package is installed, otherwise install
   if (find.package("glmmSeq", quiet = TRUE) %>% length %>% equals(0)) {
@@ -1136,7 +1136,7 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data,
   # Create design matrix for dispersion, removing random effects
   design =
     model.matrix(
-      object = .formula |> eliminate_random_effects(),
+      object = .formula |> lme4::nobars(),
       data = metadata
     )
 
@@ -1154,11 +1154,16 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data,
   # Make sure the order matches the counts
   dispersion = dispersion[rownames(counts)]
 
+  # Scaling
+  sizeFactors <- counts |> edgeR::calcNormFactors(method = scaling_method)
+
+
   glmmSeq_object =
     glmmSeq( .formula,
                       countdata = counts ,
                       metadata =   metadata |> as.data.frame(),
                       dispersion = dispersion,
+             sizeFactors = sizeFactors,
                       progress = TRUE,
                       method = method |> str_remove("(?i)^glmmSeq_" ),
                       ...

R/utilities.R

@@ -1524,32 +1524,3 @@ feature__ =  get_special_column_name_symbol(".feature")
 sample__ = get_special_column_name_symbol(".sample")
 
 
-#' Produced using ChatGPT - Eliminate Random Effects from a Formula
-#'
-#' This function takes a mixed-effects model formula and returns a modified
-#' formula with all random effects removed, leaving only the fixed effects.
-#'
-#' @param formula An object of class \code{formula}, representing a mixed-effects model formula.
-#' @return A formula object with random effects parts removed.
-#' @examples
-#' eliminate_random_effects(~ age_days * sex + (1 | file_id) + ethnicity_simplified + assay_simplified + .aggregated_cells + (1 + age_days * sex | tissue))
-#' @noRd
-#' @importFrom stats as.formula
-eliminate_random_effects <- function(formula) {
-  # Convert the formula to a string
-  formula_str <- deparse(formula)
-
-  # Split the string by '+' while keeping the brackets content together
-  split_str <- unlist(strsplit(formula_str, "(?<=\\))\\s*\\+\\s*|\\+\\s*(?=\\()", perl = TRUE))
-
-  # Filter out the random effects parts
-  fixed_effects <- grep("\\|", split_str, value = TRUE, invert = TRUE)
-
-  # Combine the fixed effects parts back into a single string
-  modified_str <- paste(fixed_effects, collapse = " + ")
-
-  # Convert the modified string back to a formula
-  modified_formula <- as.formula(modified_str)
-
-  return(modified_formula)
-}