@@ -23,15 +23,15 @@

   ## because we keep the argument 'parallel.sz' for backwards compatibility

   ## we need to harmonize it with the contents of BPPARAM

-  if (parallel.sz > 1L && class(BPPARAM) == "SerialParam") {

+  if (parallel.sz > 1L && inherits(BPPARAM, "SerialParam")) {

     BPPARAM=MulticoreParam(progressbar=verbose, workers=parallel.sz, tasks=100)

-  } else if (parallel.sz == 1L && class(BPPARAM) != "SerialParam") {

+  } else if (parallel.sz == 1L && !inherits(BPPARAM, "SerialParam")) {

     parallel.sz <- bpnworkers(BPPARAM)

-  } else if (parallel.sz > 1L && class(BPPARAM) != "SerialParam") {

+  } else if (parallel.sz > 1L && !inherits(BPPARAM, "SerialParam")) {

     bpworkers(BPPARAM) <- parallel.sz

   }

-  if (class(BPPARAM) != "SerialParam" && verbose)

+  if (!inherits(BPPARAM, "SerialParam") && verbose)

     cat(sprintf("Setting parallel calculations through a %s back-end\nwith workers=%d and tasks=100.\n",

                 class(BPPARAM), parallel.sz))

@@ -139,7 +139,7 @@ rankHDF5 <- function(X){

     write_block(sink, grid, block)

   }

-  sink <- viewportReduce(colRanks_byBlock, grid, sink)

+  sink <- gridReduce(colRanks_byBlock, grid, sink)

   close(sink)

   res <- as(sink, "DelayedArray")

   res

@@ -62,7 +62,8 @@

     return(plageDelayed(expr, gset.idx.list, parallel.sz, verbose, BPPARAM=BPPARAM))

   }

-  

+

+  stop("Not yet implemented for 'method=\"gsva\"'")

 }

 h5BackendRealization <- function(gSetIdx, FUN, Z, bpp) {

@@ -149,7 +149,7 @@ rankHDF5 <- function(X){

 .fastRndWalk2 <- function(gSetIdx, geneRanking, ra_block) {

   n <- length(geneRanking)

   k <- length(gSetIdx)

-  idxs <- sort.int(fastmatch::fmatch(gSetIdx, geneRanking))

+  idxs <- sort.int(match(gSetIdx, geneRanking))

   stepCDFinGeneSet2 <-

     sum(ra_block[geneRanking[idxs]] * (n - idxs + 1)) /

     sum((ra_block[geneRanking[idxs]]))    

@@ -116,15 +116,13 @@ zscoreDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   Z <- t(DelayedArray::scale(t(X)))

-  es <- bplapply(geneSets, h5BackendRealization, combinezDelayed,

-                 Z, BPPARAM = BPPARAM)

-  

-  es <- do.call(rbind, es)

-  rownames(es) <- names(geneSets)

-  colnames(es) <- colnames(X)

-  

+  es <- bplapply(geneSets, function(gSetIdx, Z){

+    x <- combinezDelayed(gSetIdx, Z)

+    x <- matrix(x, 1, length(x))

+    x <- writeHDF5Array(x)

+  }, Z, BPPARAM = BPPARAM)

+  es <- do.call(DelayedArray::rbind, es)

   es <- as(es, "HDF5Array")

-  

   es

 }

@@ -183,7 +181,7 @@ ssgseaDelayed <- function(X, geneSets, alpha=0.25, parallel.sz,

     res

   }, BPPARAM=BPPARAM)

-  es <- do.call(cbind, es)

+  es <- do.call(DelayedArray::cbind, es)

   es

   if (normalization) {

@@ -85,7 +85,7 @@ h5BackendRealization <- function(gSetIdx, FUN, Z, bpp) {

 }

-rightsingularsvdvectorgset <- function(gSetIdx, Z, bpp) {

+rightsvdDelayed <- function(gSetIdx, Z, bpp) {

   s <-runRandomSVD(Z[gSetIdx,], k=1, nu=0, BPPARAM = bpp)

   s$v[, 1]

 }

@@ -95,7 +95,7 @@ plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   Z <- t(DelayedArray::scale(t(X)))

-  es <- bplapply(geneSets, h5BackendRealization, rightsingularsvdvectorgset,

+  es <- bplapply(geneSets, h5BackendRealization, rightsvdDelayed,

                  Z, bpp=BPPARAM, BPPARAM=BPPARAM)

   es <- do.call(rbind, es)

@@ -171,8 +171,6 @@ ssgseaDelayed <- function(X, geneSets, alpha=0.25, parallel.sz,

   Ra <- abs(R)^alpha

-  Ra <- as(Ra, "HDF5Array")

-  

   es <- bplapply(as.list(1:n), function(j) {

     geneRanking <- order(R[, j], decreasing=TRUE)

     colRa <- Ra[,j]

@@ -65,7 +65,7 @@

 }

-h5BackendRealization <- function(gSetIdx, FUN, Z) {

+h5BackendRealization <- function(gSetIdx, FUN, Z, bpp) {

   FUN <- match.fun(FUN)

@@ -74,7 +74,7 @@ h5BackendRealization <- function(gSetIdx, FUN, Z) {

   # step 2: create grid over sink

   sink_grid <- rowAutoGrid(sink, nrow = 1)

   # step 3: create block using FUN and write it on sink

-  block <- FUN(gSetIdx, Z)

+  block <- FUN(gSetIdx, Z, bpp)

   block <- matrix(block, 1, length(block))

   sink <- DelayedArray::write_block(sink, sink_grid[[1L]], block)

   # step 4: close the sink as an hdf5Array

@@ -85,8 +85,8 @@ h5BackendRealization <- function(gSetIdx, FUN, Z) {

 }

-rightsingularsvdvectorgset <- function(gSetIdx, Z) {

-  s <- svd(Z[gSetIdx, ])

+rightsingularsvdvectorgset <- function(gSetIdx, Z, bpp) {

+  s <-runRandomSVD(Z[gSetIdx,], k=1, nu=0, BPPARAM = bpp)

   s$v[, 1]

 }

@@ -96,7 +96,7 @@ plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   Z <- t(DelayedArray::scale(t(X)))

   es <- bplapply(geneSets, h5BackendRealization, rightsingularsvdvectorgset,

-                 Z, BPPARAM=BPPARAM)

+                 Z, bpp=BPPARAM, BPPARAM=BPPARAM)

   es <- do.call(rbind, es)

   rownames(es) <- names(geneSets)

@@ -70,9 +70,9 @@ h5BackendRealization <- function(gSetIdx, FUN, Z) {

   FUN <- match.fun(FUN)

   # step 1: create realization sink

-  sink <- HDF5Array::HDF5RealizationSink(dim = c(1L, ncol(Z)))

+  sink <- HDF5RealizationSink(dim = c(1L, ncol(Z)))

   # step 2: create grid over sink

-  sink_grid <- DelayedArray::rowAutoGrid(sink, nrow = 1)

+  sink_grid <- rowAutoGrid(sink, nrow = 1)

   # step 3: create block using FUN and write it on sink

   block <- FUN(gSetIdx, Z)

   block <- matrix(block, 1, length(block))

@@ -131,26 +131,22 @@ zscoreDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

 #### rank function for hdf5 files using sink and grid methods

 rankHDF5 <- function(X){

   sink <- HDF5RealizationSink(dim(X))

-  sink_grid <- colAutoGrid(sink, ncol = 1)

-  X_grid <- colAutoGrid(X, ncol=1)

+  grid <- defaultAutoGrid(sink, block.shape="first-dim-grows-first")

-  

-  FUN <- function(sink_grid, sink){

-    bid <- currentBlockId()

-    X_viewport <- X_grid[[bid]]

-    block <- read_block(X, X_viewport)

-    block <- as.integer(rank(block))

-    block <- matrix(block, length(block), 1)

-    write_block(sink, sink_grid, block)

+  colRanks_byBlock <- function(grid, sink){

+    block <- read_block(X, grid)

+    block <- t(colRanks(block, ties.method = "average"))

+    mode(block) <- "integer"

+    write_block(sink, grid, block)

   }

-  sink <- viewportReduce(FUN, sink_grid, sink)

+  sink <- viewportReduce(colRanks_byBlock, grid, sink)

   close(sink)

   res <- as(sink, "DelayedArray")

   res

 }

-## slightly modified .fastRndWalk() for porpuse of only 

+## slightly modified .fastRndWalk() for porpoise of only 

 ## receiving a vector and not a matrix column for sums

 .fastRndWalk2 <- function(gSetIdx, geneRanking, ra_block) {

   n <- length(geneRanking)

@@ -164,58 +160,45 @@ rankHDF5 <- function(X){

   walkStat

 }

-delayedGeneRanking <- function(r_block, ra_block, geneSets, BPPARAM){

-  geneRanking <- order(r_block, decreasing=TRUE)

-  res <- bplapply(geneSets, .fastRndWalk2, geneRanking, ra_block, BPPARAM = BPPARAM)

-  return(unlist(res))

-}

 ssgseaDelayed <- function(X, geneSets, alpha=0.25, parallel.sz,

                    normalization=TRUE, verbose=TRUE,

                    BPPARAM=SerialParam(progressbar=verbose)) {

-

+  

   n <- ncol(X)

   R <- rankHDF5(X)

   Ra <- abs(R)^alpha

-  # step1: creating a grid for both DelayedArray objects

-  # that will be iterated

-  R_grid <- colAutoGrid(R, ncol=1)

-  Ra_grid <- colAutoGrid(Ra, ncol=1)

-  

-  # step2: create a sink and a grid for it

-  sink <- HDF5RealizationSink(dim = c(1L, ncol(R)))

-  sink_grid <- colAutoGrid(sink, ncol=1)

-  

-  # step3: function that will read blocks of R and

-  # Ra grids, apply transformation and write into

-  # sink object

-  FUN <- function(sink_grid, sink, geneSets){

-    bid <- currentBlockId()

-    r_block <- read_block(R, R_grid[[bid]])

-    ra_block <- read_block(Ra, Ra_grid[[bid]])

-    res_block <- delayedGeneRanking(r_block, ra_block, geneSets, BPPARAM)

-    res_block <- matrix(res_block, 1, length(res_block))

-    write_block(sink, sink_grid, res_block)

-  }

-  

-  sink <- viewportReduce(FUN, sink_grid, sink, geneSets)

+  Ra <- as(Ra, "HDF5Array")

-  # step4: close sink as a DelayedArray object

-  close(sink)

-  es <- as(sink, "DelayedArray")

+  es <- bplapply(as.list(1:n), function(j) {

+    geneRanking <- order(R[, j], decreasing=TRUE)

+    colRa <- Ra[,j]

+    sink <- HDF5RealizationSink(c(length(names(geneSets)), 1L))

+    sink_grid <- colAutoGrid(sink, ncol=1)

+    es_sample <- lapply(geneSets, .fastRndWalk2, geneRanking, colRa)

+    sink <- DelayedArray::write_block(sink, sink_grid[[1L]], do.call("rbind", es_sample))

+    DelayedArray::close(sink)

+    res <- as(sink, "DelayedArray")

+    res

+  }, BPPARAM=BPPARAM)

+  

+  es <- do.call(cbind, es)

+  es

   if (normalization) {

     ## normalize enrichment scores by using the entire data set, as indicated

     ## by Barbie et al., 2009, online methods, pg. 2

     sink <- HDF5RealizationSink(dim(es))

-    sink_grid <- colAutoGrid(sink, ncol=1)

-    es_grid <- colAutoGrid(es, ncol=1)

+    sink_grid <- defaultAutoGrid(sink, block.shape="first-dim-grows-first")

+    es_grid <- defaultAutoGrid(es, block.shape="first-dim-grows-first")

+    fin <- range(es)[2]

+    ini <- range(es)[1]

     for(bid in seq_along(sink_grid)){

       block <- read_block(es, es_grid[[bid]])

-      block <- block / (range(es)[2] - range(es)[1])

+      block <- apply(block, 2, function(x) x / ( fin - ini))

       write_block(sink, sink_grid[[bid]], block)

     }

     close(sink)

@@ -34,6 +34,14 @@

   if (class(BPPARAM) != "SerialParam" && verbose)

     cat(sprintf("Setting parallel calculations through a %s back-end\nwith workers=%d and tasks=100.\n",

                 class(BPPARAM), parallel.sz))

+  

+  if (method == "ssgsea") {

+    if(verbose)

+      cat("Estimating ssGSEA scores for", length(gset.idx.list),"gene sets.\n")

+    

+    return(ssgseaDelayed(expr, gset.idx.list, alpha=tau, parallel.sz=parallel.sz,

+                  normalization=ssgsea.norm, verbose=verbose, BPPARAM=BPPARAM))

+  }

   if (method == "zscore") {

     if (rnaseq)

@@ -118,4 +126,106 @@ zscoreDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   es <- as(es, "HDF5Array")

   es

-}

\ No newline at end of file

+}

+

+#### rank function for hdf5 files using sink and grid methods

+rankHDF5 <- function(X){

+  sink <- HDF5RealizationSink(dim(X))

+  sink_grid <- colAutoGrid(sink, ncol = 1)

+  X_grid <- colAutoGrid(X, ncol=1)

+  

+  

+  FUN <- function(sink_grid, sink){

+    bid <- currentBlockId()

+    X_viewport <- X_grid[[bid]]

+    block <- read_block(X, X_viewport)

+    block <- as.integer(rank(block))

+    block <- matrix(block, length(block), 1)

+    write_block(sink, sink_grid, block)

+  }

+  

+  sink <- viewportReduce(FUN, sink_grid, sink)

+  close(sink)

+  res <- as(sink, "DelayedArray")

+  res

+}

+

+## slightly modified .fastRndWalk() for porpuse of only 

+## receiving a vector and not a matrix column for sums

+.fastRndWalk2 <- function(gSetIdx, geneRanking, ra_block) {

+  n <- length(geneRanking)

+  k <- length(gSetIdx)

+  idxs <- sort.int(fastmatch::fmatch(gSetIdx, geneRanking))

+  stepCDFinGeneSet2 <-

+    sum(ra_block[geneRanking[idxs]] * (n - idxs + 1)) /

+    sum((ra_block[geneRanking[idxs]]))    

+  stepCDFoutGeneSet2 <- (n * (n + 1) / 2 - sum(n - idxs + 1)) / (n - k)

+  walkStat <- stepCDFinGeneSet2 - stepCDFoutGeneSet2

+  walkStat

+}

+

+delayedGeneRanking <- function(r_block, ra_block, geneSets, BPPARAM){

+  geneRanking <- order(r_block, decreasing=TRUE)

+  res <- bplapply(geneSets, .fastRndWalk2, geneRanking, ra_block, BPPARAM = BPPARAM)

+  return(unlist(res))

+}

+

+ssgseaDelayed <- function(X, geneSets, alpha=0.25, parallel.sz,

+                   normalization=TRUE, verbose=TRUE,

+                   BPPARAM=SerialParam(progressbar=verbose)) {

+

+  n <- ncol(X)

+  

+  R <- rankHDF5(X)

+  

+  Ra <- abs(R)^alpha

+  

+  # step1: creating a grid for both DelayedArray objects

+  # that will be iterated

+  R_grid <- colAutoGrid(R, ncol=1)

+  Ra_grid <- colAutoGrid(Ra, ncol=1)

+  

+  # step2: create a sink and a grid for it

+  sink <- HDF5RealizationSink(dim = c(1L, ncol(R)))

+  sink_grid <- colAutoGrid(sink, ncol=1)

+  

+  # step3: function that will read blocks of R and

+  # Ra grids, apply transformation and write into

+  # sink object

+  FUN <- function(sink_grid, sink, geneSets){

+    bid <- currentBlockId()

+    r_block <- read_block(R, R_grid[[bid]])

+    ra_block <- read_block(Ra, Ra_grid[[bid]])

+    res_block <- delayedGeneRanking(r_block, ra_block, geneSets, BPPARAM)

+    res_block <- matrix(res_block, 1, length(res_block))

+    write_block(sink, sink_grid, res_block)

+  }

+  

+  sink <- viewportReduce(FUN, sink_grid, sink, geneSets)

+  

+  # step4: close sink as a DelayedArray object

+  close(sink)

+  es <- as(sink, "DelayedArray")

+  

+  if (normalization) {

+    ## normalize enrichment scores by using the entire data set, as indicated

+    ## by Barbie et al., 2009, online methods, pg. 2

+    sink <- HDF5RealizationSink(dim(es))

+    sink_grid <- colAutoGrid(sink, ncol=1)

+    es_grid <- colAutoGrid(es, ncol=1)

+    for(bid in seq_along(sink_grid)){

+      block <- read_block(es, es_grid[[bid]])

+      block <- block / (range(es)[2] - range(es)[1])

+      write_block(sink, sink_grid[[bid]], block)

+    }

+    close(sink)

+    es <- as(sink, "DelayedArray")

+  }

+  

+  rownames(es) <- names(geneSets)

+  colnames(es) <- colnames(X)

+  

+  es <- as(es, "HDF5Array")

+  es

+}

+

@@ -57,27 +57,19 @@

 }

-rightsingularsvdvectorgset <- function(gSetIdx, Z) {

-  s <- svd(Z[gSetIdx, ])

-  s$v[, 1]

-}

-

-svdDelayed <- function(gSetIdx, Z) {

+h5BackendRealization <- function(gSetIdx, FUN, Z) {

-  # step 1: creating a HDF5 sink object

-  sink <- HDF5Array::HDF5RealizationSink(dim = c(1L, ncol(Z)))

+  FUN <- match.fun(FUN)

-  #step 2: creating the sink's grid

+  # step 1: create realization sink

+  sink <- HDF5Array::HDF5RealizationSink(dim = c(1L, ncol(Z)))

+  # step 2: create grid over sink

   sink_grid <- DelayedArray::rowAutoGrid(sink, nrow = 1)

-  

-  # step 3: creating the block and writing it in the sink

-  # object using the sink's grid

-  block <- rightsingularsvdvectorgset(gSetIdx, Z)

+  # step 3: create block using FUN and write it on sink

+  block <- FUN(gSetIdx, Z)

   block <- matrix(block, 1, length(block))

   sink <- DelayedArray::write_block(sink, sink_grid[[1L]], block)

-  

-  # step 4: closing the sink and realizing it

-  # in a hdf5 array object

+  # step 4: close the sink as an hdf5Array

   DelayedArray::close(sink)

   res <- as(sink, "DelayedArray")

@@ -85,13 +77,18 @@ svdDelayed <- function(gSetIdx, Z) {

 }

+rightsingularsvdvectorgset <- function(gSetIdx, Z) {

+  s <- svd(Z[gSetIdx, ])

+  s$v[, 1]

+}

+

 plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

                   BPPARAM=SerialParam(progressbar=verbose)) {

   Z <- t(DelayedArray::scale(t(X)))

-  es <- bplapply(geneSets, svdDelayed, Z,

-                 BPPARAM=BPPARAM)

+  es <- bplapply(geneSets, h5BackendRealization, rightsingularsvdvectorgset,

+                 Z, BPPARAM=BPPARAM)

   es <- do.call(rbind, es)

   rownames(es) <- names(geneSets)

@@ -102,26 +99,23 @@ plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   es

 }

-combinez2 <- function(gSetIdx, Z){

+combinezDelayed <- function(gSetIdx, Z){

   DelayedMatrixStats::colSums2(Z[gSetIdx,]) / sqrt(length(gSetIdx))

 }

 zscoreDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

-                   BPPARAM=SerialParam(progressbar=verbose)) {

+                          BPPARAM=SerialParam(progressbar=verbose)){

   Z <- t(DelayedArray::scale(t(X)))

-  sink <- HDF5Array::HDF5RealizationSink(dim=c(length(names(geneSets)), ncol(X)),

-                                         dimnames = list(names(geneSets), colnames(X)))

+  es <- bplapply(geneSets, h5BackendRealization, combinezDelayed,

+                 Z, BPPARAM = BPPARAM)

-  sink_grid <- DelayedArray::rowAutoGrid(sink, nrow = 1)

+  es <- do.call(rbind, es)

+  rownames(es) <- names(geneSets)

+  colnames(es) <- colnames(X)

+  

+  es <- as(es, "HDF5Array")

-  for(bid in seq_along(sink_grid)){

-    block <- combinez2(geneSets[[bid]], Z)

-    block <- matrix(block, 1, length(block))

-    sink <- DelayedArray::write_block(sink, sink_grid[[bid]], block)

-  }

-  DelayedArray::close(sink)

-  es <- as(sink, "DelayedArray")

   es

-}

+}

\ No newline at end of file

@@ -94,6 +94,8 @@ plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

                  BPPARAM=BPPARAM)

   es <- do.call(rbind, es)

+  rownames(es) <- names(geneSets)

+  colnames(es) <- colnames(X)

   es <- as(es, "HDF5Array")

@@ -35,6 +35,15 @@

     cat(sprintf("Setting parallel calculations through a %s back-end\nwith workers=%d and tasks=100.\n",

                 class(BPPARAM), parallel.sz))

+  if (method == "zscore") {

+    if (rnaseq)

+      stop("rnaseq=TRUE does not work with method='zscore'.")

+    

+    if(verbose)

+      cat("Estimating combined z-scores for", length(gset.idx.list), "gene sets.\n")

+    

+    return(zscoreDelayed(expr, gset.idx.list, parallel.sz, verbose, BPPARAM=BPPARAM))

+  }

   if (method == "plage") {

     if (rnaseq)

@@ -91,3 +100,26 @@ plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

   es

 }

+combinez2 <- function(gSetIdx, Z){

+  DelayedMatrixStats::colSums2(Z[gSetIdx,]) / sqrt(length(gSetIdx))

+}

+

+zscoreDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

+                   BPPARAM=SerialParam(progressbar=verbose)) {

+  

+  Z <- t(DelayedArray::scale(t(X)))

+  

+  sink <- HDF5Array::HDF5RealizationSink(dim=c(length(names(geneSets)), ncol(X)),

+                                         dimnames = list(names(geneSets), colnames(X)))

+  

+  sink_grid <- DelayedArray::rowAutoGrid(sink, nrow = 1)

+  

+  for(bid in seq_along(sink_grid)){

+    block <- combinez2(geneSets[[bid]], Z)

+    block <- matrix(block, 1, length(block))

+    sink <- DelayedArray::write_block(sink, sink_grid[[bid]], block)

+  }

+  DelayedArray::close(sink)

+  es <- as(sink, "DelayedArray")

+  es

+}

@@ -46,7 +46,6 @@

     return(plageDelayed(expr, gset.idx.list, parallel.sz, verbose, BPPARAM=BPPARAM))

   }

-

 }

 rightsingularsvdvectorgset <- function(gSetIdx, Z) {

@@ -68,7 +68,7 @@ svdDelayed <- function(gSetIdx, Z) {

   block <- matrix(block, 1, length(block))

   sink <- DelayedArray::write_block(sink, sink_grid[[1L]], block)

-  # step 4: closing the sink and realizating it

+  # step 4: closing the sink and realizing it

   # in a hdf5 array object

   DelayedArray::close(sink)

   res <- as(sink, "DelayedArray")

@@ -79,7 +79,7 @@ svdDelayed <- function(gSetIdx, Z) {

 plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

                   BPPARAM=SerialParam(progressbar=verbose)) {

-  

+

   Z <- t(DelayedArray::scale(t(X)))

   es <- bplapply(geneSets, svdDelayed, Z,

new file mode 100644

@@ -0,0 +1,94 @@

+.gsvaDelayedArray <- function(expr, gset.idx.list,

+                  method=c("gsva", "ssgsea", "zscore", "plage"),

+                  kcdf=c("Gaussian", "Poisson", "none"),

+                  rnaseq=FALSE,

+                  abs.ranking=FALSE,

+                  parallel.sz=1L,

+                  mx.diff=TRUE,

+                  tau=1,

+                  kernel=TRUE,

+                  ssgsea.norm=TRUE,

+                  verbose=TRUE,

+                  BPPARAM=SerialParam(progressbar=verbose)) {

+  

+  if (length(gset.idx.list) == 0)

+    stop("The gene set list is empty! Filter may be too stringent.")

+  

+  if (any(lengths(gset.idx.list) == 1))

+    warning("Some gene sets have size one. Consider setting 'min.sz > 1'.")

+  

+  parallel.sz <- as.integer(parallel.sz)

+  if (parallel.sz < 1L)

+    parallel.sz <- 1L

+  

+  ## because we keep the argument 'parallel.sz' for backwards compatibility

+  ## we need to harmonize it with the contents of BPPARAM

+  if (parallel.sz > 1L && class(BPPARAM) == "SerialParam") {

+    BPPARAM=MulticoreParam(progressbar=verbose, workers=parallel.sz, tasks=100)

+  } else if (parallel.sz == 1L && class(BPPARAM) != "SerialParam") {

+    parallel.sz <- bpnworkers(BPPARAM)

+  } else if (parallel.sz > 1L && class(BPPARAM) != "SerialParam") {

+    bpworkers(BPPARAM) <- parallel.sz

+  }

+  

+  if (class(BPPARAM) != "SerialParam" && verbose)

+    cat(sprintf("Setting parallel calculations through a %s back-end\nwith workers=%d and tasks=100.\n",

+                class(BPPARAM), parallel.sz))

+

+  

+  if (method == "plage") {

+    if (rnaseq)

+      stop("rnaseq=TRUE does not work with method='plage'.")

+    

+    if(verbose)

+      cat("Estimating PLAGE scores for", length(gset.idx.list),"gene sets.\n")

+    

+    return(plageDelayed(expr, gset.idx.list, parallel.sz, verbose, BPPARAM=BPPARAM))

+  }

+  

+

+}

+

+rightsingularsvdvectorgset <- function(gSetIdx, Z) {

+  s <- svd(Z[gSetIdx, ])

+  s$v[, 1]

+}

+

+svdDelayed <- function(gSetIdx, Z) {

+  

+  # step 1: creating a HDF5 sink object

+  sink <- HDF5Array::HDF5RealizationSink(dim = c(1L, ncol(Z)))

+  

+  #step 2: creating the sink's grid

+  sink_grid <- DelayedArray::rowAutoGrid(sink, nrow = 1)

+  

+  # step 3: creating the block and writing it in the sink

+  # object using the sink's grid

+  block <- rightsingularsvdvectorgset(gSetIdx, Z)

+  block <- matrix(block, 1, length(block))

+  sink <- DelayedArray::write_block(sink, sink_grid[[1L]], block)

+  

+  # step 4: closing the sink and realizating it

+  # in a hdf5 array object

+  DelayedArray::close(sink)

+  res <- as(sink, "DelayedArray")

+  

+  res

+  

+}

+

+plageDelayed <- function(X, geneSets, parallel.sz, verbose=TRUE,

+                  BPPARAM=SerialParam(progressbar=verbose)) {

+  

+  Z <- t(DelayedArray::scale(t(X)))

+  

+  es <- bplapply(geneSets, svdDelayed, Z,

+                 BPPARAM=BPPARAM)

+  

+  es <- do.call(rbind, es)

+  

+  es <- as(es, "HDF5Array")

+  

+  es

+}

+