Browse code

adding support to DelayedArray backend

pablo-rodr-bio2 authored on 19/11/2020 10:20:24
Showing 2 changed files

1 1
new file mode 100644
... ...
@@ -0,0 +1,512 @@
1
+x1 <- t(scale(t(m)))
2
+x1
3
+m[gene.sets[[1]],]
4
+x1[gene.sets[[1]], ]
5
+x1 <- t(scale(t(m)))
6
+x1 <- svd(x1)
7
+x1 <- x1$v[,1]
8
+x1
9
+x2 <- runExactSVD(m, scale = TRUE, center = TRUE)
10
+x2 <- x2$v[,1]
11
+x2
12
+x2 <- runExactSVD(t(m), scale = TRUE, center = TRUE)
13
+x2 <- x2$v[,1]
14
+x2
15
+x1
16
+t(m)
17
+scale(t(m))
18
+x2 <- runExactSVD(m, scale = TRUE, center = TRUE)
19
+x1 <- t(scale(t(m)))
20
+x1 <- svd(x1)
21
+x1 <- x1$v[,1]
22
+x1
23
+x2 <- runExactSVD(m, scale = colMeans(m), center = TRUE)
24
+x2 <- x2$v[,1]
25
+x2
26
+colMeans(m)
27
+x2 <- runExactSVD(m, scale = colSds(m), center = TRUE)
28
+x2 <- x2$v[,1]
29
+x2
30
+x1
31
+x2
32
+x2 <- runExactSVD(m, scale = colSds(m), center = colMeans(m))
33
+x2 <- x2$v[,1]
34
+x1
35
+x2
36
+colSds(m)
37
+colMeans(m)
38
+x2 <- runExactSVD(t(m), scale = colSds(m), center = colMeans(m))
39
+x2 <- x2$v[,1]
40
+x2
41
+x1
42
+x2 <- runExactSVD(m, scale = colSds(t(m)), center = colMeans(t(m)))
43
+x2 <- x2$v[,1]
44
+x2
45
+x1
46
+x1 <- t(scale(t(m)))
47
+x1
48
+colSds(m)
49
+rowSds(m)
50
+x1 <- t(scale(t(m)))
51
+x1 <- svd(x1)
52
+x1 <- x1$v[,1]
53
+x1
54
+x2 <- runExactSVD(m, center= colMeans(m), scale= rowSds(m))
55
+x2 <- x2$v[,1]
56
+x2
57
+x2 <- runExactSVD(m, center= colMeans(m), scale= rowSds(m))
58
+x2$v[,1]
59
+x3 <- runSVD(m, k=Inf)
60
+x4 <- runExactSVD(m)
61
+x3
62
+x4 <- runExactSVD(m)
63
+x4
64
+x1 <- t(scale(t(m)))
65
+x1 <- svd(x1[gene.sets[[1]],])
66
+x1 <- x1$v[,1]
67
+x1
68
+x2 <- runExactSVD(m[gene.sets[[1]],], center= colMeans(m), scale= rowSds(m))
69
+x2 <- x2$v[,1]
70
+x2
71
+x3 <- runSVD(m[gene.sets[[1]],], k=Inf)
72
+x3$v[,1]
73
+x2
74
+x2 <- runExactSVD(m[gene.sets[[1]],], center= colMeans(m), scale= rowSds(m))
75
+x2
76
+x1
77
+x1 <- t(scale(t(m)))
78
+x1
79
+x1 <- svd(x1[gene.sets[[1]],])
80
+x1
81
+x2 <- runExactSVD(m[gene.sets[[1]],], center= colMeans(m), scale= rowSds(m))
82
+x2
83
+rowSds(m)
84
+x1 <- t(scale(t(m)))
85
+x1 <- svd(x1[gene.sets[[1]],])
86
+x1 <- x1$v[,1]
87
+x1
88
+x2 <- runExactSVD(m[gene.sets[[1]],], center= colMeans(m), scale= rowSds(t(m)))
89
+x2 <- x2$v[,1]
90
+x2
91
+colMeans(m)
92
+apply(m, 2, sd)
93
+x1 <- t(scale(t(m)))
94
+x1 <- svd(x1[gene.sets[[1]],])
95
+x1 <- x1$v[,1]
96
+x1
97
+x2 <- runExactSVD(m[gene.sets[[1]],], center= colMeans(m), scale= TRUE)
98
+x2 <- x2$v[,1]
99
+x2
100
+x1 <- scale(m)
101
+x1 <- t(scale(t(m)))
102
+x1 <- svd(x1[1,])
103
+x1
104
+x1 <- svd(x1[1,,drop=FALSE])
105
+x1 <- svd(x1[1, , drop=FALSE])
106
+x1[1, , drop=FALSE]
107
+x1[1,  drop=FALSE]
108
+x1 <- t(scale(t(m)))
109
+x1 <- t(scale(t(m)))
110
+x1 <- svd(x1[1,  ,drop=FALSE])
111
+x1
112
+x2 <- runExactSVD(m[1, ])
113
+x2 <- runExactSVD(m[1, , drop=FALSE])
114
+x2
115
+x2 <- runExactSVD(m[1, , drop=FALSE], center=TRUE, scale=TRUE)
116
+x2
117
+x1 <- t(scale(t(m)))
118
+x1 <- svd(x1[1:2,])
119
+x1
120
+x2 <- runExactSVD(m[1:2,], center=TRUE, scale=TRUE)
121
+x2
122
+x2 <- runExactSVD(m[1:2,], center=colMeans(t(m)), scale=TRUE)
123
+x2
124
+x <- gene.sets[[1]]
125
+x
126
+x1 <- svd(x1[x,])
127
+x <- gene.sets[[1]]
128
+x1 <- t(scale(t(m)))
129
+x1 <- svd(x1[x,])
130
+x1 <- x1$v[,1]
131
+x1
132
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= colSds(t(m)))
133
+x2 <- x2$v[,1]
134
+x2
135
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= colSds(m))
136
+x2 <- x2$v[,1]
137
+x2
138
+x1 <- t(scale(t(m)))
139
+x1 <- svd(x1[x,])
140
+x1 <- x1$v[,1]
141
+x1
142
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(t(m)))
143
+x2 <- x2$v[,1]
144
+x2
145
+x2 <- runExactSVD(m[x,], center= colMeans(m), scale= rowSds(t(m)))
146
+x2 <- x2$v[,1]
147
+x2
148
+x2 <- runExactSVD(m[x,], center= colMeans(m), scale= colSds(t(m)))
149
+x2 <- x2$v[,1]
150
+x2
151
+x1
152
+x2 <- runExactSVD(m[x,], center= colMeans(m), scale= colSds(m))
153
+x2 <- x2$v[,1]
154
+x2
155
+x1
156
+scale(t(m))
157
+colMeans(t(m))
158
+colSds(m)
159
+colSds(t(m))
160
+rowSds(m)
161
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(m))
162
+x2 <- x2$v[,1]
163
+x2
164
+x1
165
+x1 <- t(scale(t(m)))
166
+x1 <- svd(x1[x,])
167
+x1
168
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(m))
169
+x2
170
+min(dim(m[x,]))
171
+x2 <- runSvd(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m))
172
+x2 <- runSVD(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m))
173
+x2
174
+x1$v[,1]
175
+x2$v[,1]
176
+x2 <- runSVD(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m), BSPARAM = ExactParam())
177
+x2$v[,1]
178
+x2 <- runSVD(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m), BSPARAM = FastAutoParam())
179
+x2$v[,1]
180
+x2 <- runSVD(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m), BSPARAM = bsparam())
181
+x2$v[,1]
182
+x2 <- runSVD(m[x,], k=2, center= colMeans(t(m)), scale= rowSds(m), BSPARAM = ExactParam(deferred=TRUE))
183
+x2$v[,1]
184
+x1 <- t(scale(t(m)))
185
+x1
186
+colMeans(t(m))
187
+owSds(m)
188
+rowSds(m)
189
+m
190
+m <- matrix(sample.int(10, 25, T), 5, 5)
191
+colnames(m) <- paste0("cell_", 1:5)
192
+rownames(m) <- paste0("gene_", 1:5)
193
+set.seed(123)
194
+m <- matrix(sample.int(10, 25, T), 5, 5)
195
+colnames(m) <- paste0("cell_", 1:5)
196
+rownames(m) <- paste0("gene_", 1:5)
197
+gene.sets <- list("my_list1"= c(1,2))
198
+m <- matrix(sample.int(10, 25, T), 5, 5)
199
+colnames(m) <- paste0("cell_", 1:5)
200
+rownames(m) <- paste0("gene_", 1:5)
201
+x <- c(1,2)
202
+x1 <- t(scale(t(m)))
203
+x1 <- svd(x1[x,])
204
+x1$v[,1]
205
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(m))
206
+x2$v[,1]
207
+rowSds(m)
208
+scale(t(m))
209
+t(scale(t(m)))
210
+x1 <- t(scale(t(m)))
211
+x1@scaled
212
+str(x)
213
+str(x1)
214
+x1$scaled
215
+x1@scaled
216
+attr(x1)
217
+attr(x1, "scaled:scale")
218
+x2 <- runExactSVD(t(m[x,]), center= colMeans(t(m)), scale= rowSds(m))
219
+x2$v[,1]
220
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(m))
221
+x2$v[,1]
222
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(t(m)))
223
+x1$v[,1]
224
+x1 <- t(scale(t(m)))
225
+x1 <- svd(x1[x,])
226
+x1$v[,1]
227
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(t(m)))
228
+x2$v[,1]
229
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= rowSds(m))
230
+x2$v[,1]
231
+C <- colMeans(t(m))
232
+i <- 1
233
+sqrt(sum((m[,i] - C[i])^2)/(ncol(m)-1))
234
+x2 <- runExactSVD(m[x,], center= colMeans(t(m)), scale= TRUE)
235
+x2$v[,1]
236
+x1$v[,1]
237
+colMeans(t(m))
238
+rowMeans(m)
239
+rowSds(m)
240
+x1 <- t(scale(t(m)))
241
+x1 <- svd(x1[x,])
242
+x1$v[,1]
243
+x2 <- runExactSVD(m[x,], center= rowMeans(m), scale= rowSds(m))
244
+x2$v[,1]
245
+library(BiocSingular)
246
+set.seed(123)
247
+m <- matrix(sample.int(10, 25, T), 5, 5)
248
+x <- c(1,2)
249
+x1 <- t(scale(t(m)))
250
+x1 <- svd(x1[x,])
251
+x1$v[,1]
252
+x2 <- runExactSVD(t(scale(t(m))))
253
+x2$v[,1]
254
+x1 <- t(scale(t(m)))
255
+res1 <- svd(x1[x,])
256
+res1$v[,1]
257
+x2 <- runExactSVD(x1[x,])
258
+x2$v[,1]
259
+x3 <- runExactSVD(m[x,], center=colMeans(t(m)))
260
+x3
261
+library(matrixStats)
262
+x1 <- t(scale(t(m)))
263
+x2 <- t( (t(m) - colMeans(t(m))) /  colSds(t(m)) )
264
+x1
265
+x2
266
+x2 <- t( (t(m) - colMeans(t(m))) /  rowSds(t(m)) )
267
+x2
268
+x1 <- t(scale(t(m)))
269
+x2 <- t( (t(m) - colMeans(t(m))) /  apply(t(m), 2, sd) )
270
+x1
271
+x2
272
+apply(t(m), 2, sd)
273
+m
274
+x1
275
+t(m)
276
+scale(m)
277
+m
278
+x1 <- t(scale(t(m)))
279
+x1
280
+x1 <- t(scale(t(m)))
281
+x2 <- t( (t(m) - colMeans(t(m))) /  apply(t(m), 2, sd) )
282
+x1
283
+x2
284
+colMeans(t(m))
285
+apply(t(m), 2, sd)
286
+m
287
+x1 <- scale(t(m))
288
+x2 <- (t(m) - colMeans(t(m))) /  apply(t(m), 2, sd)
289
+x1
290
+x2
291
+x2 <- (t(m) - colMeans(t(m))) /  (apply(t(m), 2, sd))
292
+x1
293
+x2
294
+x1 <- scale(m)
295
+x2 <- (m - colMeans(m)) /  (apply(m, 2, sd))
296
+x1
297
+x2
298
+library(bench)
299
+m <- matrix(sample.int(10, 25, T), 10, 10)
300
+genes <- c(1,2)
301
+x <- t(scale(t(m)))
302
+m <- matrix(sample.int(10, 25, T), 10, 10)
303
+genes <- c(1,2)
304
+bench::mark(
305
+svd(x[genes,]),
306
+runExactSVD(x[genes,])
307
+)
308
+library(Matrix)
309
+m<-rsparsematrix(10,10,.5)
310
+m
311
+# m <- matrix(sample.int(10, 25, T), 10, 10)
312
+m<-rsparsematrix(1000,10000,.5)
313
+genes <- c(1:400)
314
+x <- t(scale(t(m)))
315
+X
316
+x
317
+x1 <- t(scale(t(m)))
318
+x2 <- as(x1, "dgCMatrix")
319
+bench::mark(
320
+svd(x1[genes,]),
321
+runExactSVD(x2[genes,])
322
+)
323
+bench::mark(
324
+svd(x1[genes,]),
325
+runExactSVD(m[genes,])
326
+)
327
+bench::mark(
328
+svd(x1[genes,]),
329
+runExactSVD(m[genes,]),
330
+check=FALSE
331
+)
332
+bench::mark(
333
+svd(x1[genes,]),
334
+runExactSVD(x1[genes,]),
335
+check=FALSE
336
+)
337
+x1 <- NULL
338
+x2 <- NULL
339
+m <- NULL
340
+x <- NULL
341
+library(BiocSingular)
342
+set.seed(123)
343
+m <- matrix(sample.int(10, 25, T), 10, 10)
344
+genes <- 1:2
345
+genes
346
+genes <- c(1:2)
347
+genes
348
+genes <- 1:2
349
+m <- matrix(sample.int(10, 25, T), 10, 10)
350
+genes <- 1:2
351
+x1 <- t(scale(t(m)))
352
+res1 <- svd(x1[genes,])
353
+res$v[,1]
354
+res1$v[,1]
355
+m
356
+t(m)
357
+apply(t(m),2,sd)
358
+res2 <- runExactSVD(m, center=colMeans(t(m)), scale=apply(t(m),2,sd))
359
+res2$v[,1]
360
+res1$v[,1]
361
+res2 <- runExactSVD(m, center=colMeans(t(m)), scale=apply(m,2,sd))
362
+res2$v[,1]
363
+svd( t( (t(m) - C) / colSds(t(m)) )[genes,] )
364
+t(m) - colMeans(t(m))
365
+t(m)
366
+colMeans(t(m))
367
+colSds(t(m))
368
+( t(m) - colMeans(t(m)) ) / colSds(t(m))
369
+sweep(t(m), 2, colMeanst(m))
370
+sweep(t(m), 2, colMeans(t(m)))
371
+svd( t( sweep(t(m), 2, colMeans(t(m))) / colSds(t(m)) ) )
372
+svd( t( sweep(t(m), 2, colMeans(t(m))) / colSds(t(m)) )[genes,] )
373
+svd( t( sweep(t(m), 2, colMeans(t(m))) / colSds(t(m)) )[genes,] )$v[,1]
374
+x1 <- t(scale(t(m)))
375
+res1 <- svd(x1[genes,])
376
+res1$v[,1]
377
+svd( t( sweep(t(m), 2, colMeans(t(m))) / rowSds(t(m)) )[genes,] )$v[,1]
378
+res1$v[,1]
379
+X <- m
380
+Z <- Matrix::t(X)
381
+Z <- .dgCapply(Z, scale, 2)
382
+Z <- Matrix::t(Z)
383
+Z
384
+library(Matrix)
385
+library(SingleCellExperiment)
386
+library(BiocParallel)
387
+library(sparseMatrixStats)
388
+m<-rsparsematrix(10,10,.5)
389
+colnames(m) <- paste0("cell_", 1:10)
390
+rownames(m) <- paste0("gene_", 1:10)
391
+# gset.idx.list <- list("my_genes" = sample(rownames(m), 2))
392
+geneSets <- list("my_genes1" = c(1,3))
393
+library(GSVA)
394
+x <- gsva(m, geneSets, method="plage")
395
+colnames(m) <- paste0("cell_", 1:10)
396
+rownames(m) <- paste0("gene_", 1:10)
397
+x <- gsva(m, geneSets, method="plage")
398
+# gset.idx.list <- list("my_genes" = sample(rownames(m), 2))
399
+geneSets <- list("my_genes1" = c("gene_1","gene_3"))
400
+x <- gsva(m, geneSets, method="plage")
401
+x
402
+x <- gsva(m, geneSets, method="zscore")
403
+x
404
+x <- gsva(m, geneSets, method="ssgsea")
405
+x
406
+X <- m
407
+n <- ncol(X)
408
+if(is(X, "dgCMatrix")){
409
+R <- t(sparseMatrixStats::colRanks(X, ties.method = "average"))
410
+mode(R) <- "integer"
411
+} else {
412
+R <- apply(X, 2, function(x, p) as.integer(rank(x)), p)
413
+}
414
+R
415
+R <- t(sparseMatrixStats::colRanks(X, ties.method = "average"))
416
+mode(R) <- "integer"
417
+R
418
+R <- apply(X, 2, function(x, p) as.integer(rank(x)), p)
419
+R
420
+R <- t(sparseMatrixStats::colRanks(X, ties.method = "average"))
421
+source('~/curro/gsva-devel/devel3_ssgsea.R', echo=TRUE)
422
+mode(R) <- "integer"
423
+R
424
+Ra <- abs(R)^alpha
425
+Ra
426
+es <- bplapply(as.list(1:n), function(j) {
427
+geneRanking <- order(R[, j], decreasing=TRUE)
428
+es_sample <- lapply(geneSets, .fastRndWalk, geneRanking, j, Ra)
429
+unlist(es_sample)
430
+}, BPPARAM=BPPARAM)
431
+.fastRndWalk <- function(gSetIdx, geneRanking, j, Ra) {
432
+n <- length(geneRanking)
433
+k <- length(gSetIdx)
434
+idxs <- sort.int(fastmatch::fmatch(gSetIdx, geneRanking))
435
+stepCDFinGeneSet2 <-
436
+sum(Ra[geneRanking[idxs], j] * (n - idxs + 1)) /
437
+sum((Ra[geneRanking[idxs], j]))
438
+stepCDFoutGeneSet2 <- (n * (n + 1) / 2 - sum(n - idxs + 1)) / (n - k)
439
+walkStat <- stepCDFinGeneSet2 - stepCDFoutGeneSet2
440
+walkStat
441
+}
442
+n <- ncol(X)
443
+es <- bplapply(as.list(1:n), function(j) {
444
+geneRanking <- order(R[, j], decreasing=TRUE)
445
+es_sample <- lapply(geneSets, .fastRndWalk, geneRanking, j, Ra)
446
+unlist(es_sample)
447
+}, BPPARAM=BPPARAM)
448
+es
449
+geneSets <- 1,3
450
+geneSets <- c(1,3)
451
+es <- bplapply(as.list(1:n), function(j) {
452
+geneRanking <- order(R[, j], decreasing=TRUE)
453
+es_sample <- lapply(geneSets, .fastRndWalk, geneRanking, j, Ra)
454
+unlist(es_sample)
455
+}, BPPARAM=BPPARAM)
456
+es
457
+es <- do.call("cbind", es)
458
+es
459
+geneSets <- 1
460
+es <- bplapply(as.list(1:n), function(j) {
461
+geneRanking <- order(R[, j], decreasing=TRUE)
462
+es_sample <- lapply(geneSets, .fastRndWalk, geneRanking, j, Ra)
463
+unlist(es_sample)
464
+}, BPPARAM=BPPARAM)
465
+es
466
+es <- do.call("cbind", es)
467
+es
468
+if (normalization) {
469
+## normalize enrichment scores by using the entire data set, as indicated
470
+## by Barbie et al., 2009, online methods, pg. 2
471
+es <- apply(es, 2, function(x, es) x / (range(es)[2] - range(es)[1]), es)
472
+}
473
+normalization=TRUE
474
+if (normalization) {
475
+## normalize enrichment scores by using the entire data set, as indicated
476
+## by Barbie et al., 2009, online methods, pg. 2
477
+es <- apply(es, 2, function(x, es) x / (range(es)[2] - range(es)[1]), es)
478
+}
479
+es
480
+rownames(es) <- names(geneSets)
481
+colnames(es) <- colnames(X)
482
+geneSets <- c(1,3)
483
+es <- bplapply(as.list(1:n), function(j) {
484
+geneRanking <- order(R[, j], decreasing=TRUE)
485
+es_sample <- lapply(geneSets, .fastRndWalk, geneRanking, j, Ra)
486
+unlist(es_sample)
487
+}, BPPARAM=BPPARAM)
488
+es <- do.call("cbind", es)
489
+if (normalization) {
490
+## normalize enrichment scores by using the entire data set, as indicated
491
+## by Barbie et al., 2009, online methods, pg. 2
492
+es <- apply(es, 2, function(x, es) x / (range(es)[2] - range(es)[1]), es)
493
+}
494
+es
495
+if(is(X, "dgCMatrix")){
496
+es <- as(es, "dgCMatrix")
497
+}
498
+es
499
+rownames(es) <- names(geneSets)
500
+colnames(es) <- colnames(X)
501
+if(is(X, "dgCMatrix")){
502
+es <- as(es, "dgCMatrix")
503
+}
504
+es
505
+setwd("~/curro/gsva-devel/GSVA")
506
+devtools::check(vignettes = FALSE)
507
+devtools::check(vignettes = FALSE)
508
+devtools::build(vignettes = FALSE)
509
+devtools::install(build_vignettes = FALSE)
510
+x <- gsva(sce, gset.idx.list, method="zscore")
511
+x <- gsva(sce, gset.idx.list, method="ssgsea")
512
+X
... ...
@@ -1,7 +1,7 @@
1 1
 .filterFeatures <- function(expr, method) {
2 2
 
3 3
   ## filter out genes with constant expression values
4
-  sdGenes <- apply(expr, 1, sd)
4
+  sdGenes <- DelayedMatrixStats::rowSds(expr)
5 5
   if (any(sdGenes == 0) || any(is.na(sdGenes))) {
6 6
     warning(sum(sdGenes == 0 | is.na(sdGenes)),
7 7
             " genes with constant expression values throuhgout the samples.")