@@ -1,14 +1,13 @@

 Package: GSVA

-Version: 1.35.1

+Version: 1.35.2

 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"),

              person("Joan", "Fernandez", role="ctb", email="joanfernandez1331@gmail.com"))

 Depends: R (>= 3.5.0)

 Imports: methods, BiocGenerics, Biobase, GSEABase (>= 1.17.4),

-         geneplotter, shiny, shinythemes

-Suggests: limma, RColorBrewer, genefilter,

-          edgeR, snow, parallel, GSVAdata

+         parallel, BiocParallel, geneplotter, shiny, shinythemes

+Suggests: limma, RColorBrewer, genefilter, edgeR, snow, GSVAdata

 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.

 License: GPL (>= 2)

 LazyLoad: yes

@@ -9,11 +9,17 @@ importClassesFrom(GSEABase, GeneSetCollection)

 importMethodsFrom(Biobase, featureNames,

                            phenoData,

-                           experimentData)

+                           experimentData,

+                           esApply,

+                           exprs,

+                           annotation)

 importMethodsFrom(GSEABase, geneIds,

                             incidence)

+importMethodsFrom(BiocParallel, bpiterate,

+                                "bpworkers<-")

+

 importFrom(graphics, plot)

 importFrom(stats, ecdf,

                   na.omit)

@@ -21,11 +27,14 @@ importFrom(utils, setTxtProgressBar,

                   txtProgressBar,

                   read.csv,

                   write.csv)

-importFrom(Biobase, exprs,

-                    annotation)

 importFrom(GSEABase, AnnoOrEntrezIdentifier,

                      mapIdentifiers,

                      getGmt)

+importFrom(parallel, splitIndices)

+importFrom(BiocParallel, SerialParam,

+                         MulticoreParam,

+                         multicoreWorkers,

+                         bpnworkers)

 importFrom(geneplotter, multidensity)

 importFrom(shinythemes, shinytheme)

@@ -12,18 +12,19 @@ setMethod("gsva", signature(expr="ExpressionSet", gset.idx.list="list"),

   abs.ranking=FALSE,

   min.sz=1,

   max.sz=Inf,

-  parallel.sz=0, 

+  parallel.sz=1L, 

   parallel.type="SOCK",

   mx.diff=TRUE,

   tau=switch(method, gsva=1, ssgsea=0.25, NA),

   ssgsea.norm=TRUE,

-  verbose=TRUE)

+  verbose=TRUE,

+  BPPARAM=SerialParam(progressbar=verbose))

 {

   method <- match.arg(method)

   kcdf <- match.arg(kcdf)

   ## filter out genes with constant expression values

-  sdGenes <- Biobase::esApply(expr, 1, sd)

+  sdGenes <- esApply(expr, 1, sd)

   if (any(sdGenes == 0) || any(is.na(sdGenes))) {

     warning(sum(sdGenes == 0 | is.na(sdGenes)),

             " genes with constant expression values throuhgout the samples.")

@@ -62,7 +63,8 @@ setMethod("gsva", signature(expr="ExpressionSet", gset.idx.list="list"),

   }

   eSco <- .gsva(exprs(expr), mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking,

-                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm, verbose)

+                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm,

+                verbose, BPPARAM) 

   eScoEset <- new("ExpressionSet", exprs=eSco, phenoData=phenoData(expr),

                   experimentData=experimentData(expr), annotation="")

@@ -79,18 +81,19 @@ setMethod("gsva", signature(expr="ExpressionSet", gset.idx.list="GeneSetCollecti

   abs.ranking=FALSE,

   min.sz=1,

   max.sz=Inf,

-  parallel.sz=0, 

+  parallel.sz=1L, 

   parallel.type="SOCK",

   mx.diff=TRUE,

   tau=switch(method, gsva=1, ssgsea=0.25, NA),

   ssgsea.norm=TRUE,

-  verbose=TRUE)

+  verbose=TRUE,

+  BPPARAM=SerialParam(progressbar=verbose))

 {

   method <- match.arg(method)

   kcdf <- match.arg(kcdf)

   ## filter out genes with constant expression values

-  sdGenes <- Biobase::esApply(expr, 1, sd)

+  sdGenes <- esApply(expr, 1, sd)

   if (any(sdGenes == 0) || any(is.na(sdGenes))) {

     warning(sum(sdGenes == 0 | is.na(sdGenes)),

             " genes with constant expression values throuhgout the samples.")

@@ -107,8 +110,10 @@ setMethod("gsva", signature(expr="ExpressionSet", gset.idx.list="GeneSetCollecti

     cat("Mapping identifiers between gene sets and feature names\n")

   ## map gene identifiers of the gene sets to the features in the chip

-  mapped.gset.idx.list <- GSEABase::mapIdentifiers(gset.idx.list,

-                                                   GSEABase::AnnoOrEntrezIdentifier(Biobase::annotation(expr)))

+  ## Biobase::annotation() is necessary to disambiguate from the

+  ## 'annotation' argument

+  mapped.gset.idx.list <- mapIdentifiers(gset.idx.list,

+                                         AnnoOrEntrezIdentifier(Biobase::annotation(expr)))

   ## map to the actual features for which expression data is available

   tmp <- lapply(geneIds(mapped.gset.idx.list),

@@ -133,7 +138,8 @@ setMethod("gsva", signature(expr="ExpressionSet", gset.idx.list="GeneSetCollecti

   }

   eSco <- .gsva(exprs(expr), mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking,

-                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm, verbose)

+                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm,

+                verbose, BPPARAM)

   eScoEset <- new("ExpressionSet", exprs=eSco, phenoData=phenoData(expr),

                   experimentData=experimentData(expr), annotation="")

@@ -150,12 +156,13 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="GeneSetCollection"),

   abs.ranking=FALSE,

   min.sz=1,

   max.sz=Inf,

-  parallel.sz=0, 

+  parallel.sz=1L, 

   parallel.type="SOCK",

   mx.diff=TRUE,

   tau=switch(method, gsva=1, ssgsea=0.25, NA),

   ssgsea.norm=TRUE,

-  verbose=TRUE)

+  verbose=TRUE,

+  BPPARAM=SerialParam(progressbar=verbose))

 {

   method <- match.arg(method)

   kcdf <- match.arg(kcdf)

@@ -180,8 +187,8 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="GeneSetCollection"),

     if (verbose)

       cat("Mapping identifiers between gene sets and feature names\n")

-    mapped.gset.idx.list <- GSEABase::mapIdentifiers(gset.idx.list,

-                                                     GSEABase::AnnoOrEntrezIdentifier(annotation))

+    mapped.gset.idx.list <- mapIdentifiers(gset.idx.list,

+                                           AnnoOrEntrezIdentifier(annotation))

   }

   ## map to the actual features for which expression data is available

@@ -211,7 +218,8 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="GeneSetCollection"),

   }

   rval <- .gsva(expr, mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking,

-                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm, verbose)

+                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm,

+                verbose, BPPARAM)

   rval

 })

@@ -223,12 +231,13 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

   abs.ranking=FALSE,

   min.sz=1,

   max.sz=Inf,

-  parallel.sz=0, 

+  parallel.sz=1L, 

   parallel.type="SOCK",

   mx.diff=TRUE,

   tau=switch(method, gsva=1, ssgsea=0.25, NA),

   ssgsea.norm=TRUE,

-  verbose=TRUE)

+  verbose=TRUE,

+  BPPARAM=SerialParam(progressbar=verbose))

 {

   method <- match.arg(method)

   kcdf <- match.arg(kcdf)

@@ -272,7 +281,8 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

   }

   rval <- .gsva(expr, mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking,

-                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm, verbose)

+                parallel.sz, parallel.type, mx.diff, tau, kernel, ssgsea.norm,

+                verbose, BPPARAM)

   rval

 })

@@ -282,25 +292,49 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

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

   rnaseq=FALSE,

   abs.ranking=FALSE,

-  parallel.sz=0, 

+  parallel.sz=1L,

   parallel.type="SOCK",

   mx.diff=TRUE,

   tau=1,

   kernel=TRUE,

   ssgsea.norm=TRUE,

-  verbose=TRUE)

-{

-	if(length(gset.idx.list) == 0){

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

+  verbose=TRUE,

+  BPPARAM=SerialParam(progressbar=verbose)) {

+

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

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

 	}

+

+  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)

+    if (verbose)

+      message(sprintf("Setting parallel calculations through a multicore back-end with workers=%d and tasks=100.",

+              parallel.sz))

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

+    parallel.sz <- bpnworkers(BPPARAM)

+    if (verbose)

+      message(sprintf("Setting parallel calculations through a %s back-end with %d workers.",

+              class(BPPARAM), parallel.sz))

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

+    bpworkers(BPPARAM) <- parallel.sz

+    if (verbose)

+      message(sprintf("Setting parallel calculations through a %s back-end with %d workers.",

+              class(BPPARAM), parallel.sz))

+  }

+

   if (method == "ssgsea") {

 	  if(verbose)

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

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

                   parallel.type=parallel.type, normalization=ssgsea.norm,

-                  verbose=verbose))

+                  verbose=verbose, BPPARAM=BPPARAM))

   }

   if (method == "zscore") {

@@ -310,7 +344,8 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

 	  if(verbose)

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

-    return(zscore(expr, gset.idx.list, parallel.sz, parallel.type, verbose))

+    return(zscore(expr, gset.idx.list, parallel.sz, parallel.type,

+                  verbose, BPPARAM=BPPARAM))

   }

   if (method == "plage") {

@@ -320,7 +355,8 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

 	  if(verbose)

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

-    return(plage(expr, gset.idx.list, parallel.sz, parallel.type, verbose))

+    return(plage(expr, gset.idx.list, parallel.sz, parallel.type,

+                 verbose, BPPARAM=BPPARAM))

   }

 	if(verbose)

@@ -335,11 +371,12 @@ setMethod("gsva", signature(expr="matrix", gset.idx.list="list"),

 	rownames(es.obs) <- names(gset.idx.list)

 	if (verbose)

-    cat("Computing observed enrichment scores\n")

+    cat("Computing GSVA enrichment scores\n")

 	es.obs <- compute.geneset.es(expr, gset.idx.list, 1:n.samples,

-                               rnaseq=rnaseq, abs.ranking=abs.ranking, parallel.sz=parallel.sz,

-                               parallel.type=parallel.type, mx.diff=mx.diff, tau=tau,

-                               kernel=kernel, verbose=verbose)

+                               rnaseq=rnaseq, abs.ranking=abs.ranking,

+                               parallel.sz=parallel.sz, parallel.type=parallel.type,

+                               mx.diff=mx.diff, tau=tau, kernel=kernel,

+                               verbose=verbose, BPPARAM=BPPARAM)

 	colnames(es.obs) <- colnames(expr)

 	rownames(es.obs) <- names(gset.idx.list)

@@ -377,8 +414,9 @@ compute.gene.density <- function(expr, sample.idxs, rnaseq=FALSE, kernel=TRUE){

 }

 compute.geneset.es <- function(expr, gset.idx.list, sample.idxs, rnaseq=FALSE,

-                               abs.ranking, parallel.sz=0, parallel.type="SOCK",

-                               mx.diff=TRUE, tau=1, kernel=TRUE, verbose=TRUE){

+                               abs.ranking, parallel.sz=1L, parallel.type="SOCK",

+                               mx.diff=TRUE, tau=1, kernel=TRUE,

+                               verbose=TRUE, BPPARAM=SerialParam(progressbar=verbose)) {

 	num_genes <- nrow(expr)

 	if (verbose) {

     if (kernel) {

@@ -392,7 +430,7 @@ compute.geneset.es <- function(expr, gset.idx.list, sample.idxs, rnaseq=FALSE,

   if (parallel.sz > 0 && length(sample.idxs > 100) && nrow(expr) > 100) {

     cat(sprintf("Using %d cores in parallel for ECDFs estimation\n", parallel.sz))

     iter <- function(Y, n_chunks=BiocParallel::multicoreWorkers()) {

-      idx <- parallel::splitIndices(nrow(Y), min(nrow(Y), n_chunks))

+      idx <- splitIndices(nrow(Y), min(nrow(Y), n_chunks))

       i <- 0L

       function() {

         if (i == length(idx))

@@ -401,13 +439,12 @@ compute.geneset.es <- function(expr, gset.idx.list, sample.idxs, rnaseq=FALSE,

         Y[idx[[i]], , drop=FALSE]

       }

     }

-    gene.density <- BiocParallel::bpiterate(iter(expr, 100),

-                        GSVA:::compute.gene.density,

-                        sample.idxs=sample.idxs, rnaseq=rnaseq, kernel=kernel,

-                        REDUCE=rbind, reduce.in.order=TRUE,

-                        BPPARAM=BiocParallel::MulticoreParam(workers=parallel.sz,

-                                                             progressbar=TRUE,

-                                                             tasks=100))

+    gene.density <- bpiterate(iter(expr, 100),

+                              compute.gene.density,

+                              sample.idxs=sample.idxs,

+                              rnaseq=rnaseq, kernel=kernel,

+                              REDUCE=rbind, reduce.in.order=TRUE,

+                              BPPARAM=BPPARAM)

   } else 

 	  gene.density <- compute.gene.density(expr, sample.idxs, rnaseq, kernel)

@@ -613,7 +650,8 @@ setCores <- function(nCores, parallel.sz) {

 }

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

-                   parallel.type, normalization=TRUE, verbose) {

+                   parallel.type, normalization=TRUE,

+                   verbose=TRUE, BPPARAM=SerialParam(progressbar=verbose)) {

   p <- nrow(X)

   n <- ncol(X)

@@ -704,7 +742,8 @@ ssgsea <- function(X, geneSets, alpha=0.25, parallel.sz,

 combinez <- function(gSetIdx, j, Z) sum(Z[gSetIdx, j]) / sqrt(length(gSetIdx))

-zscore <- function(X, geneSets, parallel.sz, parallel.type, verbose) {

+zscore <- function(X, geneSets, parallel.sz, parallel.type,

+                   verbose=TRUE, BPPARAM=SerialParam(progressbar=verbose)) {

   p <- nrow(X)

   n <- ncol(X)

@@ -788,7 +827,8 @@ rightsingularsvdvectorgset <- function(gSetIdx, Z) {

   s$v[, 1]

 }

-plage <- function(X, geneSets, parallel.sz, parallel.type, verbose) {

+plage <- function(X, geneSets, parallel.sz, parallel.type,

+                  verbose=TRUE, BPPARAM=SerialParam(progressbar=verbose)) {

   p <- nrow(X)

   n <- ncol(X)

@@ -929,7 +969,7 @@ setMethod("computeGeneSetsOverlap", signature(gSets="list", uniqGenes="character

 setMethod("computeGeneSetsOverlap", signature(gSets="list", uniqGenes="ExpressionSet"),

           function(gSets, uniqGenes, min.sz=1, max.sz=Inf) {

-  uniqGenes <- Biobase::featureNames(uniqGenes)

+  uniqGenes <- featureNames(uniqGenes)

   totalGenes <- length(uniqGenes)

   ## map to the actual features for which expression data is available

@@ -958,9 +998,11 @@ setMethod("computeGeneSetsOverlap", signature(gSets="GeneSetCollection", uniqGen

 setMethod("computeGeneSetsOverlap", signature(gSets="GeneSetCollection", uniqGenes="ExpressionSet"),

           function(gSets, uniqGenes, min.sz=1, max.sz=Inf) {

   ## map gene identifiers of the gene sets to the features in the chip

-  gSets <- GSEABase::mapIdentifiers(gSets, GSEABase::AnnoOrEntrezIdentifier(Biobase::annotation(uniqGenes)))

+  ## Biobase::annotation() is necessary to disambiguate from the

+  ## 'annotation' argument

+  gSets <- mapIdentifiers(gSets, AnnoOrEntrezIdentifier(Biobase::annotation(uniqGenes)))

-  uniqGenes <- Biobase::featureNames(uniqGenes)

+  uniqGenes <- featureNames(uniqGenes)

   gSetsMembershipMatrix <- incidence(gSets)

   gSetsMembershipMatrix <- t(gSetsMembershipMatrix[, colnames(gSetsMembershipMatrix) %in% uniqGenes])

@@ -43,6 +43,6 @@ install_github("rcastelo/GSVA")

 For questions and bug reports regarding the __release__ version of **GSVA**

 please use the [Bioconductor support site](https://support.bioconductor.org "Bioconductor support site").

-For bug reports and issues regarding this __development__ version of **GSVA**

+For feature requests or bug reports and issues regarding this __development__ version of **GSVA**

 please use the GitHub issues link at the top-right of this page

 ([https://github.com/rcastelo/GSVA/issues](https://github.com/rcastelo/GSVA/issues)).

@@ -20,48 +20,52 @@ Estimates GSVA enrichment scores.

     abs.ranking=FALSE,

     min.sz=1,

     max.sz=Inf,

-    parallel.sz=0,

+    parallel.sz=1L,

     parallel.type="SOCK",

     mx.diff=TRUE,

     tau=switch(method, gsva=1, ssgsea=0.25, NA),

     ssgsea.norm=TRUE,

-    verbose=TRUE)

+    verbose=TRUE,

+    BPPARAM=SerialParam(progressbar=verbose))

 \S4method{gsva}{ExpressionSet,GeneSetCollection}(expr, gset.idx.list, annotation,

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

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

     abs.ranking=FALSE,

     min.sz=1,

     max.sz=Inf,

-    parallel.sz=0,

+    parallel.sz=1L,

     parallel.type="SOCK",

     mx.diff=TRUE,

     tau=switch(method, gsva=1, ssgsea=0.25, NA),

     ssgsea.norm=TRUE,

-    verbose=TRUE)

+    verbose=TRUE,

+    BPPARAM=SerialParam(progressbar=verbose))

 \S4method{gsva}{matrix,GeneSetCollection}(expr, gset.idx.list, annotation,

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

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

     abs.ranking=FALSE,

     min.sz=1,

     max.sz=Inf,

-    parallel.sz=0,

+    parallel.sz=1L,

     parallel.type="SOCK",

     mx.diff=TRUE,

     tau=switch(method, gsva=1, ssgsea=0.25, NA),

     ssgsea.norm=TRUE,

-    verbose=TRUE)

+    verbose=TRUE,

+    BPPARAM=SerialParam(progressbar=verbose))

 \S4method{gsva}{matrix,list}(expr, gset.idx.list, annotation,

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

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

     abs.ranking=FALSE,

     min.sz=1,

     max.sz=Inf,

-    parallel.sz=0,

+    parallel.sz=1L,

     parallel.type="SOCK",

     mx.diff=TRUE,

     tau=switch(method, gsva=1, ssgsea=0.25, NA),

     ssgsea.norm=TRUE,

-    verbose=TRUE)

+    verbose=TRUE,

+    BPPARAM=SerialParam(progressbar=verbose))

 }

 \arguments{

   \item{expr}{Gene expression data which can be given either as an \code{ExpressionSet}

@@ -99,14 +103,9 @@ Estimates GSVA enrichment scores.

             enriched on either extreme (high or low) will be regarded as 'highly' activated.}

   \item{min.sz}{Minimum size of the resulting gene sets.}

   \item{max.sz}{Maximum size of the resulting gene sets.}

-  \item{parallel.sz}{Number of processors to use when doing the calculations in parallel.

-                     This requires to previously load either the \code{parallel} or the

-                     \code{snow} library. If \code{parallel} is loaded and this argument

-                     is left with its default value (\code{parallel.sz=0}) then it will use

-                     all available core processors unless we set this argument with a

-                     smaller number. If \code{snow} is loaded then we must set this argument

-                     to a positive integer number that specifies the number of processors to

-                     employ in the parallel calculation.}

+  \item{parallel.sz}{Number of threads of execution to use when doing the calculations in parallel.

+                     The argument \code{BPPARAM} allows one to set the parallel back-end and fine

+                     tune its configuration.}

   \item{parallel.type}{Type of cluster architecture when using \code{snow}.}

   \item{mx.diff}{Offers two approaches to calculate the enrichment statistic (ES)

                  from the KS random walk statistic. \code{mx.diff=FALSE}: ES is calculated as

@@ -122,6 +121,7 @@ Estimates GSVA enrichment scores.

                      the minimum and the maximum, as described in their paper. When \code{ssgsea.norm=FALSE}

                      this last normalization step is skipped.}

   \item{verbose}{Gives information about each calculation step. Default: \code{FALSE}.}

+  \item{BPPARAM}{An object of class \code{\link[BiocParallel]{BiocParallelParam}} specifiying parameters related to the parallel execution of some of the tasks and calculations within this function.}

 }

 \details{