Browse code

Bugfixes and tests for the new quickCorrect function.

LTLA authored on 06/03/2021 07:55:59
Showing 4 changed files

... ...
@@ -72,7 +72,7 @@ intersectRows <- function(..., subset.row=NULL, keep.all=FALSE) {
72 72
             stop("only character 'subset.row' is allowed when '...' have different rownames")
73 73
         }
74 74
         if (!keep.all) {
75
-            all.batches[[i]] <- lapply(all.batches, function(x) x[subset.row,,drop=FALSE])
75
+            all.batches <- lapply(all.batches, function(x) x[subset.row,,drop=FALSE])
76 76
         }
77 77
     }
78 78
 
... ...
@@ -66,11 +66,12 @@ quickCorrect <- function(...,
66 66
     assay.type="counts", 
67 67
     PARAM=FastMnnParam(), 
68 68
     multi.norm.args=list(),
69
-    precomputed.var=NULL,
69
+    precomputed=NULL,
70 70
     model.var.args=list(),
71 71
     hvg.args=list(n=5000))
72 72
 {
73 73
     all.batches <- intersectRows(...)
74
+    single.batch <- length(all.batches)==1L
74 75
 
75 76
     # Perform multibatchNorm.
76 77
     all.batches <- do.call(multiBatchNorm, c(
... ...
@@ -80,9 +81,9 @@ quickCorrect <- function(...,
80 81
     ))
81 82
 
82 83
     # Perform HVG calling within each batch.
83
-    if (is.null(precomputed.var)) {
84
-        if (length(all.batches) == 1) {
85
-            dec <- do.call(scran::modelGeneVar, c(list(all.batches[[1]], block=batch), model.var.args))
84
+    if (is.null(precomputed)) {
85
+        if (single.batch) { 
86
+            dec <- do.call(scran::modelGeneVar, c(list(all.batches, block=batch), model.var.args))
86 87
         } else {
87 88
             sub.dec <- vector("list", length(all.batches))
88 89
             for (i in seq_along(all.batches)) {
... ...
@@ -91,14 +92,17 @@ quickCorrect <- function(...,
91 92
             dec <- do.call(scran::combineVar, sub.dec)
92 93
         }
93 94
     } else {
94
-        if (length(precomputed.var)==length(all.batches)) {
95
-            stop("'length(precomputed.var)' is not the same as the number of objects in '...'")
96
-        }
97
-        if (length(all.batches) == 1) {
98
-            dec <- precomputed.var[[1]]
95
+        if (single.batch) {
96
+            if (length(precomputed)!=1L){ 
97
+                stop("'length(precomputed)' is not the same as the number of objects in '...'")
98
+            }
99
+            dec <- precomputed[[1]]
99 100
         } else {
101
+            if (length(precomputed)!=length(all.batches)) {
102
+                stop("'length(precomputed)' is not the same as the number of objects in '...'")
103
+            }
100 104
             universe <- rownames(all.batches[[1]])
101
-            sub.dec <- lapply(precomputed.var, function(x) x[universe,,drop=FALSE])
105
+            sub.dec <- lapply(precomputed, function(x) x[universe,,drop=FALSE])
102 106
             dec <- do.call(scran::combineVar, sub.dec)
103 107
         }
104 108
     }
... ...
@@ -12,7 +12,7 @@ quickCorrect(
12 12
   assay.type = "counts",
13 13
   PARAM = FastMnnParam(),
14 14
   multi.norm.args = list(),
15
-  precomputed.var = NULL,
15
+  precomputed = NULL,
16 16
   model.var.args = list(),
17 17
   hvg.args = list(n = 5000)
18 18
 )
... ...
@@ -45,14 +45,14 @@ and \linkS4class{RegressParam} will dispatch to \code{\link{regressBatches}}.}
45 45
 
46 46
 \item{multi.norm.args}{Named list of further arguments to pass to \code{\link{multiBatchNorm}}.}
47 47
 
48
+\item{precomputed}{List of \link{DataFrame}s containing precomputed variance modelling results.
49
+This should be a list of the same length as the number of entries in \code{...},
50
+and each should have the same row names as the corresponding entry of \code{...}.}
51
+
48 52
 \item{model.var.args}{Named list of further arguments to pass to \code{\link[scran]{modelGeneVar}}.}
49 53
 
50 54
 \item{hvg.args}{Named list of further arguments to pass to \code{\link[scran]{getTopHVGs}}.
51 55
 By default, we take the top 5000 genes with the highest variances.}
52
-
53
-\item{precomputed}{List of \link{DataFrame}s containing precomputed variance modelling results.
54
-This should be a list of the same length as the number of entries in \code{...},
55
-and each should have the same row names as the corresponding entry of \code{...}.}
56 56
 }
57 57
 \value{
58 58
 A list containing:
59 59
new file mode 100644
... ...
@@ -0,0 +1,52 @@
1
+# This tests the various quickCorrect functions.
2
+# library(testthat); library(batchelor); source("test-quick-correct.R")
3
+
4
+set.seed(100)
5
+d1 <- matrix(rnbinom(50000, mu=10, size=1), ncol=100)
6
+sce1 <- SingleCellExperiment(list(counts=d1))
7
+sizeFactors(sce1) <- runif(ncol(d1))
8
+rownames(sce1) <- paste0("GENE", 1:500)
9
+
10
+d2 <- matrix(rnbinom(20000, mu=50, size=1), ncol=40)
11
+sce2 <- SingleCellExperiment(list(counts=d2))
12
+sizeFactors(sce2) <- runif(ncol(d2))
13
+rownames(sce2) <- paste0("GENE", 201:700)
14
+
15
+set.seed(1000)
16
+output <- quickCorrect(sce1, sce2)
17
+universe <- intersect(rownames(sce1), rownames(sce2))
18
+
19
+test_that("quickCorrect works as expected", {
20
+    expect_identical(rownames(output$corrected), universe)
21
+
22
+    # Same results.
23
+    normed <- multiBatchNorm(sce1[universe,], sce2[universe,])
24
+    dec1 <- scran::modelGeneVar(normed[[1]])
25
+    dec2 <- scran::modelGeneVar(normed[[2]])
26
+
27
+    set.seed(1000)
28
+    pre <- quickCorrect(sce1, sce2, precomputed=list(dec1, dec2))
29
+    expect_equal(reducedDim(output$corrected), reducedDim(pre$corrected))
30
+
31
+    expect_error(quickCorrect(sce1, sce2, precomputed=list(dec2)), "not the same")
32
+})
33
+
34
+test_that("quickCorrect works with a single object", {
35
+    com <- cbind(sce1[universe,], sce2[universe,])
36
+    b <- rep(1:2, c(ncol(sce1), ncol(sce2)))
37
+
38
+    set.seed(1000)
39
+    single <- quickCorrect(com, batch=b)
40
+    expect_equal(reducedDim(output$corrected), reducedDim(single$corrected))
41
+
42
+    # Same results.
43
+    normed <- multiBatchNorm(com, batch=b)
44
+    dec <- scran::modelGeneVar(normed, block=b)
45
+
46
+    set.seed(1000)
47
+    pre <- quickCorrect(com, batch=b, precomputed=list(dec))
48
+    expect_equal(reducedDim(output$corrected), reducedDim(pre$corrected))
49
+
50
+    expect_error(quickCorrect(com, batch=b, precomputed=list(dec, dec)), "not the same")
51
+})
52
+