@@ -1,5 +1,5 @@

 Package: GSVA

-Version: 1.39.10

+Version: 1.39.11

 Title: Gene Set Variation Analysis for microarray and RNA-seq data

 Authors@R: c(person("Justin", "Guinney", role=c("aut", "cre"), email="justin.guinney@sagebase.org"),

              person("Robert", "Castelo", role="aut", email="robert.castelo@upf.edu"),

@@ -52,6 +52,10 @@ importFrom(BiocParallel, SerialParam,

 importFrom(SingleCellExperiment, SingleCellExperiment)

 importFrom(sparseMatrixStats, colRanks)

 importFrom(DelayedArray, rowAutoGrid,

+			  colAutoGrid,

+			  currentBlockId,

+			  read_block, 

+			  viewportReduce,

 			  write_block,

 			  close)

 importFrom(HDF5Array, HDF5RealizationSink)

@@ -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

+}

+