From 4855baabdf1d580353c477fedd851d7eef45f208 Mon Sep 17 00:00:00 2001
From: Stevie Ped <stephen.pederson.au@gmail.com>
Date: Tue, 23 Jul 2024 23:39:05 +0930
Subject: [PATCH] Test after more bug fixing in WaldTest

---
 R/fitAssayDiff.R                   | 25 ++++++++++++++-----------
 man/fitAssayDiff-methods.Rd        |  5 +++--
 tests/testthat/test_fitAssayDiff.R |  2 +-
 3 files changed, 18 insertions(+), 14 deletions(-)

diff --git a/R/fitAssayDiff.R b/R/fitAssayDiff.R
index cb5cfa2..49f05c3 100644
--- a/R/fitAssayDiff.R
+++ b/R/fitAssayDiff.R
@@ -19,7 +19,7 @@
 #' No normalisation is applied when using the limma-trend model, as this allows
 #' for previous normalisation strategies to be performed on the data.
 #'
-#' When applying the \link[DESeq2]{nbinomWaldTest}, applying RLE normalisation
+#' If testing with \link[DESeq2]{nbinomWaldTest}, applying RLE normalisation
 #' without groups, and using colSums for library sizes (instead of total
 #' alignments), the standard normalisation factors from
 #' \link[DESeq2]{estimateSizeFactors} will be used.
@@ -27,7 +27,8 @@
 #' \link[edgeR]{normLibSizes} will be used.
 #' The fitType is set to 'local' when estimating dispersions, and this can be
 #' easily modified by passing fitType via the dot arguments.
-#' Results are additionally returned after applying \link[DESeq2]{lfcShrink}.
+#' Results are additionally returned after applying \link[DESeq2]{lfcShrink},
+#' including the svalue returned by this approach.
 #'
 #' Normalising to ChIP Input samples is not yet implemented.
 #' Similarly, the use of offsets when applying the Wald test is not yet
@@ -193,19 +194,21 @@ setMethod(
             dds <- .se2Wald(
                 x, assay, design, lib.size, norm, groups, offset, weighted, ...
             )
-            res <- DESeq2::results(dds)
-            p_mu0 <- res$pvalue
-            pval <- p_mu0
-            ## Apply lfcShrink as the default
-            res <- DESeq2::lfcShrink(
+            res <- DESeq2::results(dds, name = coef, lfcThreshold = 0)
+            if (abs(lfc) > 0) {
+                p_mu0 <- res$pvalue
+                res <- DESeq2::results(dds, name = coef, lfcThreshold = abs(lfc))
+            }
+            ## Always return shrunken estimates. Maybe make optional later
+            shrink <- DESeq2::lfcShrink(
                 dds, coef = coef, res = res, type = type,
                 lfcThreshold = abs(lfc), svalue = TRUE
             )
-            if (abs(lfc) > 0) pval <- res$svalue
-            ## Reformat for standard edgeR column layout
+            ## Collate into consistent output
             res <- data.frame(
-                logFC = res$log2FoldChange, logCPM = log2(res$baseMean),
-                PValue = pval
+                logFC = shrink$log2FoldChange, logCPM = log2(shrink$baseMean),
+                svalue = shrink$svalue, PValue = res$pvalue,
+                row.names = rownames(shrink)
             )
 
         }
diff --git a/man/fitAssayDiff-methods.Rd b/man/fitAssayDiff-methods.Rd
index c418057..6536a2d 100644
--- a/man/fitAssayDiff-methods.Rd
+++ b/man/fitAssayDiff-methods.Rd
@@ -112,7 +112,7 @@ called by column name.
 No normalisation is applied when using the limma-trend model, as this allows
 for previous normalisation strategies to be performed on the data.
 
-When applying the \link[DESeq2]{nbinomWaldTest}, applying RLE normalisation
+If testing with \link[DESeq2]{nbinomWaldTest}, applying RLE normalisation
 without groups, and using colSums for library sizes (instead of total
 alignments), the standard normalisation factors from
 \link[DESeq2]{estimateSizeFactors} will be used.
@@ -120,7 +120,8 @@ In all other scenarios, normalisation factors as returned by
 \link[edgeR]{normLibSizes} will be used.
 The fitType is set to 'local' when estimating dispersions, and this can be
 easily modified by passing fitType via the dot arguments.
-Results are additionally returned after applying \link[DESeq2]{lfcShrink}.
+Results are additionally returned after applying \link[DESeq2]{lfcShrink},
+including the svalue returned by this approach.
 
 Normalising to ChIP Input samples is not yet implemented.
 Similarly, the use of offsets when applying the Wald test is not yet
diff --git a/tests/testthat/test_fitAssayDiff.R b/tests/testthat/test_fitAssayDiff.R
index 9280ac7..d7cf005 100644
--- a/tests/testthat/test_fitAssayDiff.R
+++ b/tests/testthat/test_fitAssayDiff.R
@@ -48,7 +48,7 @@ test_that("Results appear structurally correct", {
 
     new_se <- fitAssayDiff(se, method = "wald", design = X)
     row_data <- rowData(new_se)
-    cols <- c("range", "logFC", "logCPM", "PValue", "FDR")
+    cols <- c("range", "logFC", "logCPM", "svalue", "PValue", "FDR")
     expect_equal(colnames(row_data), cols)
 
 })