Merge pull request #39 from fmicompbio/legacy-qlf

Add legacyQLF argument to runEISA (pass to edgeR::glmQLFit)

.github/workflows/R-CMD-check.yaml

@@ -29,7 +29,7 @@ jobs:
       GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       CURL_CONFIG: ${{ matrix.config.curlConfigPath }}curl-config
-      cache-version: v1
+      cache-version: v2
 
     steps:
       - name: checkout branch
@@ -79,13 +79,14 @@ jobs:
           brew install harfbuzz
           brew install fribidi
           Rscript -e 'BiocManager::install(c("GenomeInfoDbData"), type = "source")'
-          Rscript -e 'BiocManager::install(c("GenomicFeatures"), type = "binary")'
+          Rscript -e 'BiocManager::install(c("GenomicFeatures"), type = "source")'
 
 
       - name: Install R package dependencies
         run: |
           local_deps <- remotes::local_package_deps(dependencies = TRUE)
           deps <- remotes::dev_package_deps(dependencies = TRUE, repos = BiocManager::repositories())
+          print(deps)
           BiocManager::install(local_deps[local_deps %in% deps$package[deps$diff != 0]], Ncpu = 2L)
           remotes::install_cran('rcmdcheck', Ncpu = 2L)
         shell: Rscript {0}

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: eisaR
 Title: Exon-Intron Split Analysis (EISA) in R
-Version: 1.15.0
+Version: 1.15.1
 Authors@R: c(person("Michael", "Stadler", email = "[email protected]", role = c("aut", "cre")),
              person("Dimos", "Gaidatzis", email = "[email protected]", role = "aut"),
              person("Lukas", "Burger", email = "[email protected]", role = "aut"),
@@ -22,7 +22,7 @@ biocViews:
     RNASeq
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Suggests: 
     knitr,
     rmarkdown,
@@ -46,7 +46,7 @@ Imports:
     S4Vectors,
     IRanges,
     limma,
-    edgeR,
+    edgeR (>= 4.0),
     methods,
     SummarizedExperiment,
     BiocGenerics,

NEWS.md

@@ -1,3 +1,7 @@
+# eisaR 1.15.1
+
+* Add legacyQLF argument to runEISA (will be passed to edgeR::glmQLFit)
+
 # eisaR 1.5.2
 
 * Add options to collapse introns by gene and restrict introns to feature ranges in getFeatureRanges.

R/runEISA.R

@@ -57,6 +57,10 @@
 #'       \item{\code{"LRT"}: }{Likelihood ratio test using \code{\link[edgeR]{glmFit}}
 #'       and \code{\link[edgeR:glmFit]{glmLRT}}}.
 #'   }
+#' @param legacyQLF Whether to use the 'legacy' version of 
+#'   \code{\link[edgeR:glmQLFTest]{glmQLFit}}. See \code{\link[edgeR:glmQLFTest]{glmQLFit}}
+#'   for more details. If \code{FALSE}, the new method introduced in 
+#'   \code{edgeR} 4.0.0 is used.
 #' @param effects How the effects (contrasts or log2 fold-changes) are calculated.
 #'   One of:\describe{
 #'       \item{\code{"predFC"}: }{(default) Fold-changes are calculated using
@@ -143,6 +147,7 @@ runEISA <- function(cntEx, cntIn, cond, method = NULL,
                     modelSamples = TRUE,
                     geneSelection = c("filterByExpr", "none", "Gaidatzis2015"),
                     statFramework = c("QLF", "LRT"),
+                    legacyQLF = FALSE,
                     effects = c("predFC", "Gaidatzis2015"),
                     pscnt = 2, 
                     sizeFactor = c("exon", "intron", "individual"), 
@@ -346,7 +351,7 @@ runEISA <- function(cntEx, cntIn, cond, method = NULL,
             contr <- as.numeric(colnames(dsgn) == colnames(dsgn)[ncol(dsgn)])
         }
         if (statFramework == "QLF") {
-            fit <- edgeR::glmQLFit(y, dsgn)
+            fit <- edgeR::glmQLFit(y, dsgn, legacy = legacyQLF)
             tst.ExIn <- edgeR::glmQLFTest(fit, contrast = contr)
         } else if (statFramework == "LRT") {
             fit <- edgeR::glmFit(y, dsgn)
@@ -373,7 +378,7 @@ runEISA <- function(cntEx, cntIn, cond, method = NULL,
             # dsgn2 <- model.matrix(~ region * cond2)
             # y2 <- edgeR::estimateDisp(y, dsgn2)
             # if (statFramework == "QLF") {
-            #     fit2 <- edgeR::glmQLFit(y2, dsgn2)
+            #     fit2 <- edgeR::glmQLFit(y2, dsgn2, legacy = legacyQLF)
             # } else if (statFramework == "LRT") {
             #     fit2 <- edgeR::glmFit(y2, dsgn2)
             # }

tests/testthat/test_runEISA.R

@@ -18,11 +18,13 @@ test_that("runEISA() runs", {
     res0 <- runEISA(cntEx[1:1000,], cntIn[1:1000,], cond,
                     modelSamples = TRUE, geneSelection = "none", statFramework = "LRT")
     res0se <- runEISA(cntEx = cntSE2[1:1000,], cntIn = NULL, cond,
-                      modelSamples = TRUE, geneSelection = "none")
+                      modelSamples = TRUE, geneSelection = "none", legacyQLF = TRUE)
     res1 <- runEISA(cntEx, cntIn, cond, method = "Gaidatzis2015")
     res1se <- runEISA(cntEx = cntSE1, cntIn = NULL, cond, method = "Gaidatzis2015")
-    res2 <- runEISA(cntEx, cntIn, cond, method = NULL, modelSamples = FALSE, sizeFactor = "individual")
-    res3 <- runEISA(cntEx, cntIn, cond, recalcLibSizeAfterFilt = TRUE, sizeFactor = "intron")
+    res2 <- runEISA(cntEx, cntIn, cond, method = NULL, modelSamples = FALSE, 
+                    sizeFactor = "individual", legacyQLF = TRUE)
+    res3 <- runEISA(cntEx, cntIn, cond, recalcLibSizeAfterFilt = TRUE, 
+                    sizeFactor = "intron", legacyQLF = TRUE)
     expect_is(res0, "list")
     expect_is(res0se, "list")
     expect_equal(res0$contrasts, res0se$contrasts)
@@ -50,7 +52,8 @@ test_that("runEISA() runs", {
     expect_equal(nrow(res1$tab.ExIn), 0)
     expect_error(plotEISA(res1))
     expect_error(suppressWarnings(runEISA(cntEx[, c(1, 3)], cntIn[, c(1, 3)],
-                                          cond[c(1, 3)], method = NULL)))
+                                          cond[c(1, 3)], method = NULL, 
+                                          legacyQLF = TRUE)))
 })
 
 context("runEISA gives expected results")
@@ -102,9 +105,11 @@ test_that("runEISA() gives expected results", {
 
     # run EISA (use gene filtering that is independent of model)
     res1 <- runEISA(cntEx, cntIn, cond, geneSelection = "Gaidatzis2015",
-                    pscnt = 8, sizeFactor = "individual", modelSamples = FALSE)
+                    pscnt = 8, sizeFactor = "individual", modelSamples = FALSE,
+                    legacyQLF = TRUE)
     res2 <- runEISA(cntEx, cntIn, cond, geneSelection = "Gaidatzis2015",
-                    pscnt = 8, sizeFactor = "individual", modelSamples = TRUE)
+                    pscnt = 8, sizeFactor = "individual", modelSamples = TRUE, 
+                    legacyQLF = TRUE)
 
     # account for filtered genes
     ngenes <- nrow(res1$tab.ExIn)