Split GSVA ranks calculations, enabling distribution of computations

.covrignore

@@ -1,3 +1,5 @@
 inst/*
+R/GSVA-pkg-deprecated.R
+R/GSVA-pkg-defunct.R
 R/methods-ShinyApp.R
 R/zzz.R

DESCRIPTION

@@ -30,4 +30,5 @@ BugReports: https://github.com/rcastelo/GSVA/issues
 Encoding: UTF-8
 biocViews: FunctionalGenomics, Microarray, RNASeq, Pathways, GeneSetEnrichment
 Roxygen: list(markdown = TRUE)
+Config/roxygen2/version: 8.0.0
 RoxygenNote: 7.3.3

NAMESPACE

@@ -10,9 +10,12 @@ export(geneSetSizes)
 export(geneSets)
 export(gsva)
 export(gsvaAnnotation)
+export(gsvaColRanks)
+export(gsvaColScores)
 export(gsvaEnrichment)
 export(gsvaParam)
 export(gsvaRanks)
+export(gsvaRowNorm)
 export(gsvaScores)
 export(guessGeneIdType)
 export(igsva)
@@ -32,13 +35,17 @@ exportClasses(plageParam)
 exportClasses(ssgseaParam)
 exportClasses(zscoreParam)
 exportMethods(computeGeneSetsOverlap)
+exportMethods(details)
 exportMethods(filterGeneSets)
 exportMethods(geneSetSizes)
 exportMethods(geneSets)
 exportMethods(gsva)
 exportMethods(gsvaAnnotation)
+exportMethods(gsvaColRanks)
+exportMethods(gsvaColScores)
 exportMethods(gsvaEnrichment)
 exportMethods(gsvaRanks)
+exportMethods(gsvaRowNorm)
 exportMethods(gsvaScores)
 exportMethods(spatCor)
 import(methods)
@@ -60,8 +67,6 @@ importFrom(Biobase,selectSome)
 importFrom(BiocGenerics,"annotation<-")
 importFrom(BiocGenerics,"type<-")
 importFrom(BiocGenerics,annotation)
-importFrom(BiocGenerics,cbind)
-importFrom(BiocGenerics,rbind)
 importFrom(BiocGenerics,type)
 importFrom(BiocParallel,"bpprogressbar<-")
 importFrom(BiocParallel,"bpstopOnError<-")
@@ -78,25 +83,21 @@ importFrom(DelayedArray,chunkGrid)
 importFrom(DelayedArray,chunkdim)
 importFrom(DelayedArray,close)
 importFrom(DelayedArray,colAutoGrid)
-importFrom(DelayedArray,colSums)
-importFrom(DelayedArray,currentBlockId)
 importFrom(DelayedArray,defaultAutoGrid)
 importFrom(DelayedArray,getAutoBlockLength)
 importFrom(DelayedArray,gridReduce)
-importFrom(DelayedArray,read_block)
 importFrom(DelayedArray,rowAutoGrid)
 importFrom(DelayedArray,rowRanges)
 importFrom(DelayedArray,seed)
 importFrom(DelayedArray,setAutoBPPARAM)
-importFrom(DelayedArray,t)
-importFrom(DelayedArray,write_block)
 importFrom(DelayedMatrixStats,rowSds)
 importFrom(GSEABase,AnnoOrEntrezIdentifier)
 importFrom(GSEABase,ENSEMBLIdentifier)
 importFrom(GSEABase,EntrezIdentifier)
 importFrom(GSEABase,NullCollection)
 importFrom(GSEABase,NullIdentifier)
 importFrom(GSEABase,SymbolIdentifier)
+importFrom(GSEABase,details)
 importFrom(GSEABase,geneIdType)
 importFrom(GSEABase,geneIds)
 importFrom(GSEABase,mapIdentifiers)
@@ -112,7 +113,6 @@ importFrom(IRanges,match)
 importFrom(IRanges,ranges)
 importFrom(IRanges,start)
 importFrom(IRanges,width)
-importFrom(Matrix,colSums)
 importFrom(Matrix,nnzero)
 importFrom(MatrixGenerics,colRanks)
 importFrom(MatrixGenerics,colSums)
@@ -123,11 +123,14 @@ importFrom(S4Arrays,DummyArrayGrid)
 importFrom(S4Arrays,is_sparse)
 importFrom(S4Arrays,read_block)
 importFrom(S4Arrays,refdim)
+importFrom(S4Arrays,write_block)
 importFrom(S4Vectors,"metadata<-")
 importFrom(S4Vectors,DataFrame)
 importFrom(S4Vectors,SimpleList)
 importFrom(S4Vectors,metadata)
 importFrom(SingleCellExperiment,SingleCellExperiment)
+importFrom(SingleCellExperiment,altExps)
+importFrom(SingleCellExperiment,reducedDims)
 importFrom(SparseArray,nzcount)
 importFrom(SparseArray,rowMeans)
 importFrom(SparseArray,rowSds)
@@ -151,6 +154,7 @@ importFrom(graphics,plot)
 importFrom(graphics,segments)
 importFrom(memuse,Sys.meminfo)
 importFrom(memuse,howbig)
+importFrom(methods,callNextMethod)
 importFrom(methods,getClass)
 importFrom(methods,is)
 importFrom(methods,isClass)

R/AllClasses.R

@@ -6,14 +6,26 @@
 #' Virtual superclass of expression data classes supported by `GSVA`.
 #'
 #' `GSVA` supports expression data matrices in a growing number of containers
-#' and representations.  This class union allows to store any of these in a slot
-#' of another class as well as defining common methods for all of them.
+#' and representations. This class union allows to store any of these in a slot
+#' of another class as well as defining common methods for all of them. The
+#' current list of supported classes is:
+#'
+#' * `matrix`
+#' * `dgCMatrix` from package `Matrix`
+#' * `SVT_SparseMatrix` from package `SparseArray`
+#' * `DelayedMatrix` from package `DelayedArray`
+#' * `HDF5Matrix` from package `HDF5Array`
+#' * `ExpressionSet` from package `Biobase`
+#' * `SummarizedExperiment` from package `SummarizedExperiment`
+#' * `SingleCellExperiment` from package `SingleCellExperiment`
+#' * `SpatialExperiment` from package `SpatialExperiment`
 #'
 #' @seealso
 #' [`matrix`],
 #' [`dgCMatrix`][Matrix::dgCMatrix-class],
 #' [`SVT_SparseMatrix`][SparseArray::SVT_SparseMatrix-class],
 #' [`DelayedMatrix`][DelayedArray::DelayedMatrix-class],
+#' [`HDF5Array`][HDF5Array::HDF5Array-class],
 #' \code{\link[Biobase]{ExpressionSet}},
 ### we are using the plain Rd above because
 ###  #' [`ExpressionSet`][Biobase::ExpressionSet-class],
@@ -22,8 +34,7 @@
 ###  ‘[Biobase:class.ExpressionSet]{ExpressionSet}’
 #' [`SummarizedExperiment`][SummarizedExperiment::SummarizedExperiment-class],
 #' [`SingleCellExperiment`][SingleCellExperiment::SingleCellExperiment-class],
-#' [`SpatialExperiment`][SpatialExperiment::SpatialExperiment-class],
-#' [`HDF5Array`][HDF5Array::HDF5Array-class]
+#' [`SpatialExperiment`][SpatialExperiment::SpatialExperiment-class]
 #'
 #' @importClassesFrom Matrix dgCMatrix
 #' @importClassesFrom Biobase ExpressionSet

R/AllGenerics.R

@@ -3,6 +3,18 @@
 setGeneric("gsva",
            function(param, ...) standardGeneric("gsva"))
 
+#' @export 
+setGeneric("gsvaRowNorm",
+           function(param, ...) standardGeneric("gsvaRowNorm"))
+
+#' @export 
+setGeneric("gsvaColRanks",
+           function(rowNormExprData, ...) standardGeneric("gsvaColRanks"))
+
+#' @export 
+setGeneric("gsvaColScores",
+           function(rankExprData, ...) standardGeneric("gsvaColScores"))
+
 #' @export 
 setGeneric("gsvaRanks",
            function(param, ...) standardGeneric("gsvaRanks"))
@@ -13,7 +25,7 @@ setGeneric("gsvaScores",
 
 #' @export 
 setGeneric("gsvaEnrichment",
-           function(param, ranks, ...) standardGeneric("gsvaEnrichment"))
+           function(rankExprData, ...) standardGeneric("gsvaEnrichment"))
 
 #' @export 
 setGeneric("filterGeneSets",
@@ -50,7 +62,7 @@ setGeneric("unwrapData",
            function(container, ...) standardGeneric("unwrapData"))
 
 setGeneric("wrapData",
-           function(container, dataMatrix, geneSets) standardGeneric("wrapData"))
+           function(container, ...) standardGeneric("wrapData"))
 
 setGeneric("mapGeneSetsToAnno",
            function(geneSets, anno, ...) standardGeneric("mapGeneSetsToAnno"))

R/GSVA-package.R

@@ -17,7 +17,6 @@
 #' @importFrom SingleCellExperiment SingleCellExperiment
 #' @importFrom SpatialExperiment SpatialExperiment
 #' @importFrom DelayedArray rowAutoGrid colAutoGrid defaultAutoGrid
-#' currentBlockId read_block gridReduce write_block close t colSums
 #' @importFrom HDF5Array HDF5RealizationSink writeHDF5Array
 #' @importFrom cli cli_abort cli_alert_info cli_alert_warning cli_alert_success
 NULL

R/GSVA-pkg-deprecated.R

@@ -6,3 +6,202 @@
 #' @name GSVA-pkg-deprecated
 #' @keywords internal
 NULL
+
+#' @description The `gsvaRanks()` method is deprecated. Please use `gsvaRowNorm()`
+#' and `gsvaColRanks()` instead.
+#'
+#' @aliases gsvaRanks,gsvaParam-method
+#' @name gsvaRanks
+#' @rdname GSVA-pkg-deprecated
+#'
+#' @importFrom cli cli_alert_info cli_alert_success
+#' @exportMethod gsvaRanks
+setMethod("gsvaRanks", signature(param="gsvaParam"),
+          function(param,
+                   verbose=TRUE,
+                   BPPARAM=SerialParam(progressbar=verbose),
+                   maxmem="auto") {
+
+              .Deprecated(new="gsvaRowNorm() and gsvaColRanks()",
+                          package="GSVA",
+                          msg=paste("The 'gsvaRanks()' method is deprecated.",
+                                    "Please use 'gsvaRowNorm()' and",
+                                    "'gsvaColRanks()'."))
+
+              if (verbose && gsva_global$show_start_and_end_messages) {
+                  pkgversion <- packageDescription("GSVA")[["Version"]]
+                  cli_alert_info("GSVA version {pkgversion}")
+              }
+
+              .check_bpparam(BPPARAM)
+
+              exprData <- get_exprData(param)
+              dataMatrix <- unwrapData(exprData, get_assay(param))
+              maxmem <- .check_maxmem(param, maxmem=maxmem, verbose=verbose)
+              ondisk <- .check_ondisk(param, maxmem=maxmem, first=NA, last=NA,
+                                      whdim=2, verbose=verbose)
+
+              dataMatrix <- .check_sparse_load_input_expr(dataMatrix, "GSVA",
+                                                          first=NA, last=NA,
+                                                          whdim=2, ondisk,
+                                                          verbose)
+
+              filtDataMatrix <- dataMatrix
+              BPPARAM <- .check_open_parallelism(filtDataMatrix, BPPARAM,
+                                                 minparrows=100, minparcols=100,
+                                                 verbose)
+
+              if (.get_filterRows(param))
+                  filtDataMatrix <- .filterGenes(dataMatrix, anyNA(param),
+                                                 removeConstant=TRUE,
+                                                 removeNzConstant=TRUE,
+                                                 verbose, BPPARAM=BPPARAM,
+                                                 maxmem=maxmem)
+              else if (verbose) {
+                  msg <- "Skipping filtering of constant rows (filterRows=FALSE)"
+                  cli_alert_warning(msg)
+              }
+
+              if (verbose)
+                  cli_alert_info(sprintf("Calculating GSVA ranks"))
+
+              kcdfminssize <- .get_kcdfNoneMinSampleSize(param)
+              gsvarownr <- .compute_row_norm(expr=filtDataMatrix,
+                                             kcdf=.get_kcdf(param),
+                                             kcdf.min.ssize=kcdfminssize,
+                                             sparse=.get_sparse(param),
+                                             any_na=anyNA(param),
+                                             na_use=.get_NAuse(param),
+                                             verbose=verbose,
+                                             BPPARAM=BPPARAM,
+                                             maxmem=maxmem)
+
+              gsvarnks <- .compute_gsva_ranks(Z=gsvarownr,
+                                              verbose=verbose,
+                                              BPPARAM=BPPARAM,
+                                              maxmem=maxmem)
+
+              rownames(gsvarnks) <- rownames(filtDataMatrix)
+              colnames(gsvarnks) <- colnames(filtDataMatrix)
+
+              rnkscontainer <- wrapData(get_exprData(param), gsvarnks, param,
+                                        "gsvaranks", first=NA, last=NA, whdim=2,
+                                        dropAssays=FALSE)
+              rval <- new("gsvaRanksParam",
+                          exprData=rnkscontainer, geneSets=get_geneSets(param),
+                          assay="gsvaranks", annotation=get_annotation(param),
+                          minSize=get_minSize(param), maxSize=get_maxSize(param),
+                          kcdf=.get_kcdf(param),
+                          kcdfNoneMinSampleSize=.get_kcdfNoneMinSampleSize(param),
+                          tau=.get_tau(param), maxDiff=.get_maxDiff(param),
+                          absRanking=.get_absRanking(param),
+                          sparse=.get_sparse(param), checkNA=.get_checkNA(param),
+                          didCheckNA=.get_didCheckNA(param), anyNA=anyNA(param),
+                          use=.get_NAuse(param), filterRows=.get_filterRows(param),
+                          nzcount=nzcount(param), ondisk=.get_ondisk(param))
+
+              if (verbose && gsva_global$show_start_and_end_messages)
+                  cli_alert_success("Calculations finished")
+
+              return(rval)
+          })
+
+
+#' @description The `gsvaScores()` method is deprecated. Please use
+#' `gsvaColScores()` instead.
+#'
+#' @param param A parameter object of the [`gsvaRanksParam-class`] class.
+#'
+#' @aliases gsvaScores,gsvaRanksParam-method
+#' @name gsvaScores
+#' @rdname GSVA-pkg-deprecated
+#'
+#' @importFrom S4Arrays is_sparse
+#' @importFrom cli cli_alert_info cli_alert_success
+#' @exportMethod gsvaScores
+setMethod("gsvaScores", signature(param="gsvaRanksParam"),
+          function(param, verbose=TRUE,
+                   BPPARAM=SerialParam(progressbar=verbose),
+                   maxmem="auto") {
+
+              .Deprecated(new="gsvaColScores()",
+                          package="GSVA",
+                          msg=paste("The 'gsvaScores()' method is deprecated.",
+                                    "Please use 'gsvaColScores()'."))
+
+              if (verbose && gsva_global$show_start_and_end_messages) {
+                  pkgversion <- packageDescription("GSVA")[["Version"]]
+                  cli_alert_info("GSVA version {pkgversion}")
+              }
+
+              .check_bpparam(BPPARAM)
+
+              ## assuming rows in the rank data have been already filtered
+              exprData <- get_exprData(param)
+              filtDataMatrix <- unwrapData(exprData, get_assay(param))
+
+              filtMappedGeneSets <- .filterAndMapGeneSets(param=param,
+                                           filteredDataMatrix=filtDataMatrix,
+                                           verbose=verbose)
+
+              sparse <- .get_sparse(param)
+              if (sparse && !is_sparse(filtDataMatrix))
+                  sparse <- FALSE
+
+              if (verbose) {
+                if (sparse)
+                    cli_alert_info("GSVA sparse algorithm")
+                  else
+                    cli_alert_info("GSVA dense (classical) algorithm")
+              }
+
+              maxmem <- .check_maxmem(param, maxmem=maxmem, verbose=verbose)
+              ondisk <- .check_ondisk(param, maxmem=maxmem, first=NA, last=NA,
+                                      whdim=2, verbose=verbose)
+
+              filtDataMatrix <- .check_sparse_load_input_expr(filtDataMatrix,
+                                                              "GSVA", first=NA,
+                                                              last=NA, whdim=2,
+                                                              ondisk, verbose)
+
+              BPPARAM <- .check_open_parallelism(filtDataMatrix, BPPARAM,
+                                                 minparrows=100, minparcols=100,
+                                                 verbose)
+
+              ondisk <- .check_es_memory_requirements(filtDataMatrix,
+                                                      filtMappedGeneSets,
+                                                      ondisk, maxmem)
+              if (verbose) {
+                  n <- length(filtMappedGeneSets)
+                  cli_alert_info("Calculating GSVA scores for {n} gene sets")
+              }
+
+              gsva_es <- .processMatrixCols(filtDataMatrix,
+                                            FUN=.compute_gsva_scores,
+                                            geneSetsIdx=filtMappedGeneSets,
+                                            tau=.get_tau(param),
+                                            maxDiff=.get_maxDiff(param),
+                                            absRanking=.get_absRanking(param),
+                                            sparse=sparse, any_na=anyNA(param),
+                                            na_use=.get_NAuse(param),
+                                            minSize=get_minSize(param),
+                                            ondisk=ondisk, verbose=verbose,
+                                            minparrows=100, minparcols=100,
+                                            BPPARAM=BPPARAM,
+                                            maxmem=ceiling(maxmem/100)) ## use
+                                            ## of memory increases here about
+                                            ## 10-fold over block size memory
+
+              rownames(gsva_es) <- names(filtMappedGeneSets)
+              colnames(gsva_es) <- colnames(filtDataMatrix)
+
+              gs <- .geneSetsIndices2Names(indices=filtMappedGeneSets,
+                                           names=rownames(filtDataMatrix))
+              rval <- wrapData(get_exprData(param), gsva_es, param, "es",
+                               first=NA, last=NA, whdim=2, dropAssays=FALSE, gs)
+
+              if (verbose && gsva_global$show_start_and_end_messages)
+                  cli_alert_success("Calculations finished")
+
+              return(rval)
+          })