@@ -8,7 +8,8 @@ Authors@R: c(person("Justin", "Guinney", role=c("aut", "cre"), email="justin.gui

 Depends: R (>= 3.5.0)

 Imports: methods, stats, utils, graphics, BiocGenerics, S4Vectors,

          Biobase, SummarizedExperiment, GSEABase, Matrix,

-         parallel, BiocParallel, fastmatch, SingleCellExperiment

+         parallel, BiocParallel, fastmatch, SingleCellExperiment,

+         BiocSingular

 Suggests: RUnit, limma, RColorBrewer, genefilter, edgeR, GSVAdata,

           shiny, shinythemes, ggplot2, data.table, plotly

 Description: Gene Set Variation Analysis (GSVA) is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the samples of a expression data set. GSVA performs a change in coordinate systems, transforming the data from a gene by sample matrix to a gene-set by sample matrix, thereby allowing the evaluation of pathway enrichment for each sample. This new matrix of GSVA enrichment scores facilitates applying standard analytical methods like functional enrichment, survival analysis, clustering, CNV-pathway analysis or cross-tissue pathway analysis, in a pathway-centric manner.

@@ -29,6 +29,8 @@ importMethodsFrom(BiocParallel, bpiterate,

                                 bplapply,

                                 "bpprogressbar<-")

+importMethodsFrom(BiocSingular, runExactSVD)

+

 importFrom(graphics, plot)

 importFrom(stats, ecdf,

                   na.omit)

@@ -41,7 +41,12 @@ setMethod("gsva", signature(expr="SingleCellExperiment", gset.idx.list="list"),

   ## filter genes according to verious criteria,

   ## e.g., constant expression

-  expr <- .filterFeatures(expr, method)

+  if(is(expr, "dgCMatrix")){

+    expr <- .filterFeaturesSparse(expr, method)

+  } else {

+    expr <- .filterFeatures(expr, method)

+  }

+  

   if (nrow(expr) < 2)

     stop("Less than two genes in the input ExpressionSet object\n")

@@ -918,29 +923,48 @@ zscore <- function(X, geneSets, parallel.sz, verbose=TRUE,

 }

 rightsingularsvdvectorgset <- function(gSetIdx, Z) {

+  if(is(Z, "dgCMatrix")){

+    s <- BiocSingular::runExactSVD(Z[gSetIdx, ])

+  } else {

     s <- svd(Z[gSetIdx, ])

+  }

   s$v[, 1]

 }

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

                   BPPARAM=SerialParam(progressbar=verbose)) {

-

-  p <- nrow(X)

-  n <- ncol(X)

-

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

-

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

-                 BPPARAM=BPPARAM)

-

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

-

-  if (length(geneSets) == 1)

-    es <- matrix(es, nrow=1)

-

-  rownames(es) <- names(geneSets)

-  colnames(es) <- colnames(X)

-

+  if(is(X, "dgCMatrix")){

+    message("Please bear in mind that this method first scale the values of the gene expression data.

+            In order to take advantage of the sparse Matrix type, the scaling will only be applied

+            to the non-zero values of the data.")

+    

+    Z <- Matrix::t(X)

+    Z <- .dgCapply(Z, scale, 2)

+    Z <- Matrix::t(Z)

+    

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

+                   BPPARAM=BPPARAM)

+    

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

+    

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

+    

+  } else {

+    

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

+    

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

+                   BPPARAM=BPPARAM)

+    

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

+    

+    if (length(geneSets) == 1)

+      es <- matrix(es, nrow=1)

+    

+    rownames(es) <- names(geneSets)

+    colnames(es) <- colnames(X)

+  }

+  

   es

 }

@@ -35,6 +35,12 @@

   result

 }

+.dgCapply<-function(m,f, MARGIN){

+  x <- lapply(.sparseToList(m, MARGIN), f)

+  m@x <- unlist(x, use.names=FALSE)

+  m

+}

+

 ## filter out genes which non-zero values have

 ## constant expression values

 .filterFeaturesSparse <- function(expr, method) {