Few changes to ensure backward compability on the GSVA classic algorithm

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: GSVA
-Version: 1.53.25
+Version: 1.53.26
 Title: Gene Set Variation Analysis for Microarray and RNA-Seq Data
 Authors@R: c(person("Robert", "Castelo", role=c("aut", "cre"), email="[email protected]"),
              person("Justin", "Guinney", role="aut", email="[email protected]"),

R/gsva.R

@@ -322,7 +322,8 @@ compute.geneset.es <- function(expr, gset.idx.list, sample.idxs, kcdf,
 .gsva_rnd_walk <- function(gsetIdx, geneRanking, rankStat) {
   stopifnot(is.integer(gsetIdx)) ## QC
   stopifnot(is.integer(geneRanking)) ## QC
-  stopifnot(is.integer(rankStat)) ## QC
+  ## stopifnot(is.integer(rankStat)) ## QC
+  stopifnot(is.numeric(rankStat)) ## QC
   .Call("gsva_rnd_walk_R", gsetIdx, geneRanking, rankStat)
 }
 
@@ -332,7 +333,8 @@ compute.geneset.es <- function(expr, gset.idx.list, sample.idxs, kcdf,
   stopifnot(length(geneSetsRankIdx) > 0) ## QC
   stopifnot(is.integer(geneSetsRankIdx[[1]])) ## QC
   stopifnot(is.integer(geneRanking)) ## QC
-  stopifnot(is.integer(rankStat)) ## QC
+  ## stopifnot(is.integer(rankStat)) ## QC
+  stopifnot(is.numeric(rankStat)) ## QC
   stopifnot(is.logical(maxDiff)) ## QC
   stopifnot(is.logical(absRnk)) ## QC
   stopifnot(is.numeric(tau)) ## QC

R/gsvaParam.R

@@ -55,7 +55,7 @@
 #' this parameter decides at what minimum sample size `kcdf="none"`, i.e., the
 #' estimation of the empirical cumulative distribution function (ECDF) of
 #' expression levels across samples is performed directly without using a
-#' kernel; see the `kcdf` slot.
+#' kernel. By default, this value is set to 200; see the `kcdf` slot.
 #'
 #' @param tau Numeric vector of length 1.  The exponent defining the weight of
 #' the tail in the random walk performed by the `GSVA` (Hänzelmann et al.,
@@ -115,7 +115,7 @@ gsvaParam <- function(exprData, geneSets,
                       assay=NA_character_, annotation=NULL,
                       minSize=1, maxSize=Inf,
                       kcdf=c("auto", "Gaussian", "Poisson", "none"),
-                      kcdfNoneMinSampleSize=50, tau=1, maxDiff=TRUE,
+                      kcdfNoneMinSampleSize=200, tau=1, maxDiff=TRUE,
                       absRanking=FALSE, sparse=TRUE) {
     kcdf <- match.arg(kcdf)
     kcdfNoneMinSampleSize <- as.integer(kcdfNoneMinSampleSize)

src/ks_test.c

@@ -104,12 +104,12 @@ ks_matrix_R(SEXP XR, SEXP sidxsR, SEXP n_genesR, SEXP geneset_idxsR,
 }
 
 void
-gsva_rnd_walk(int* gsetidx, int k, int* generanking, int* rankstat, int n,
+gsva_rnd_walk(int* gsetidx, int k, int* generanking, double* rankstat, int n,
               double* walkstat, double* walkstatpos, double* walkstatneg) {
-  int*    stepcdfingeneset;
+  double* stepcdfingeneset;
   int*    stepcdfoutgeneset;
 
-  stepcdfingeneset = R_Calloc(n, int);  /* assuming zeroes are set */
+  stepcdfingeneset = R_Calloc(n, double);  /* assuming zeroes are set */
   stepcdfoutgeneset = R_Calloc(n, int);
   for (int i=0; i < n; i++)
     stepcdfoutgeneset[i] = 1;
@@ -153,7 +153,7 @@ gsva_rnd_walk_R(SEXP gsetidxR, SEXP generankingR, SEXP rankstatR) {
   int     k = length(gsetidxR);
   int*    gsetidx;
   int*    generanking;
-  int*    rankstat;
+  double* rankstat;
   SEXP    walkstatR;
   double* walkstat;
   double  walkstatpos, walkstatneg;
@@ -165,7 +165,7 @@ gsva_rnd_walk_R(SEXP gsetidxR, SEXP generankingR, SEXP rankstatR) {
 
   gsetidx = INTEGER(gsetidxR);
   generanking = INTEGER(generankingR);
-  rankstat = INTEGER(rankstatR);
+  rankstat = REAL(rankstatR);
   walkstat = REAL(walkstatR);
 
   gsva_rnd_walk(gsetidx, k, generanking, rankstat, n,
@@ -177,16 +177,16 @@ gsva_rnd_walk_R(SEXP gsetidxR, SEXP generankingR, SEXP rankstatR) {
 }
 
 void
-gsva_rnd_walk_nonunittau(int* gsetidx, int k, int* generanking, int* rankstat,
+gsva_rnd_walk_nonunittau(int* gsetidx, int k, int* generanking, double* rankstat,
                          int n, double tau,
                          double* walkstat, double* walkstatpos, double* walkstatneg) {
   double* stepcdfingeneset;
-  double* stepcdfoutgeneset;
+  int* stepcdfoutgeneset;
 
   stepcdfingeneset = R_Calloc(n, double);  /* assuming zeroes are set */
-  stepcdfoutgeneset = R_Calloc(n, double);
+  stepcdfoutgeneset = R_Calloc(n, int);
   for (int i=0; i < n; i++)
-    stepcdfoutgeneset[i] = 1.0;
+    stepcdfoutgeneset[i] = 1;
 
   for (int i=0; i < k; i++) {
     /* convert 1-based gene indices to 0-based ! */
@@ -230,7 +230,7 @@ gsva_score_genesets_R(SEXP genesetsrankidxR, SEXP generankingR, SEXP rankstatR,
   Rboolean absrnk=asLogical(absrnkR);
   double   tau=REAL(tauR)[0];
   int*     generanking;
-  int*     rankstat;
+  double*  rankstat;
   SEXP     esR;
   double*  es;
 
@@ -240,7 +240,7 @@ gsva_score_genesets_R(SEXP genesetsrankidxR, SEXP generankingR, SEXP rankstatR,
   PROTECT(esR = allocVector(REALSXP, m));
 
   generanking = INTEGER(generankingR);
-  rankstat = INTEGER(rankstatR);
+  rankstat = REAL(rankstatR);
   es = REAL(esR);
 
   for (int i=0; i < m; i++) {

src/utils.c

@@ -95,15 +95,15 @@ SEXP
 order_rankstat_R(SEXP xR);
 
 void
-order_rankstat(double* x, int n, int* ord, int* rst) {
+order_rankstat(double* x, int n, int* ord, double* rst) {
   for (int i=0; i < n; i++)
     ord[i] = i + 1; /* 1-based to comply with R upstream */
 
   global_dbl_p = x;
   qsort(ord, n, sizeof(int), indirect_dbl_cmp_dec);
 
   for (int i=0; i < n; i++)
-    rst[ord[i]-1] = abs(n - i - ((int) (n / 2)));
+    rst[ord[i]-1] = fabs(((double) n) - ((double) i) - (((double ) n) / 2.0));
 }
 
 SEXP
@@ -112,16 +112,16 @@ order_rankstat_R(SEXP xR) {
   double* x;
   SEXP ordR, rstR, ansR;
   int* ord;
-  int* rst;
+  double* rst;
 
   PROTECT(xR);
   x = REAL(xR);
 
   PROTECT(ordR = allocVector(INTSXP, n));
-  PROTECT(rstR = allocVector(INTSXP, n));
+  PROTECT(rstR = allocVector(REALSXP, n));
 
   ord = INTEGER(ordR);
-  rst = INTEGER(rstR);
+  rst = REAL(rstR);
 
   order_rankstat(x, n, ord, rst);
 
@@ -146,7 +146,7 @@ order_rankstat_sparse_to_dense_R(SEXP XCspR, SEXP jR) { /* column in jR is 1-bas
   int     nr;
   SEXP    ordR, rstR, ansR;
   int*    ord;
-  int*    rst;
+  double* rst;
 
   PROTECT(XCspR);
   XCsp_dim = INTEGER(GET_SLOT(XCspR, Matrix_DimSym));
@@ -161,10 +161,10 @@ order_rankstat_sparse_to_dense_R(SEXP XCspR, SEXP jR) { /* column in jR is 1-bas
     x[XCsp_i[i]] = XCsp_x[i];
 
   PROTECT(ordR = allocVector(INTSXP, nr));
-  PROTECT(rstR = allocVector(INTSXP, nr));
+  PROTECT(rstR = allocVector(REALSXP, nr));
 
   ord = INTEGER(ordR);
-  rst = INTEGER(rstR);
+  rst = REAL(rstR);
 
   order_rankstat(x, nr, ord, rst);
 
@@ -195,7 +195,7 @@ order_rankstat_sparse_to_sparse_R(SEXP XCspR, SEXP jR) { /* column in jR is 1-ba
   int     nr;
   SEXP    ordR, rstR, ansR;
   int*    ord;
-  int*    rst;
+  double* rst;
   int*    ord2;
   int*    allord;
   int     allord_i;
@@ -220,11 +220,11 @@ order_rankstat_sparse_to_sparse_R(SEXP XCspR, SEXP jR) { /* column in jR is 1-ba
     allord[i] = i + 1;
 
   PROTECT(ordR = allocVector(INTSXP, nr));
-  PROTECT(rstR = allocVector(INTSXP, nr));
+  PROTECT(rstR = allocVector(REALSXP, nr));
 
   ord2 = R_Calloc(nnz_j, int);
   ord = INTEGER(ordR);
-  rst = INTEGER(rstR);
+  rst = REAL(rstR);
 
   for (int i=0; i < nnz_j; i++)
     ord2[i] = i+1; /* indirect_dbl_cmp_dec() assumes 1-based! */
@@ -245,10 +245,10 @@ order_rankstat_sparse_to_sparse_R(SEXP XCspR, SEXP jR) { /* column in jR is 1-ba
     }
 
   for (int i=0; i < nr; i++)
-    rst[i] = (int) (nnz_j / 2) + 1; /* zero entries get the same symmtric rank */
+    rst[i] = (nnz_j / 2.0) + 1.0; /* zero entries get the same symmetric rank */
 
   for (int i=0; i < nnz_j; i++)
-    rst[ord[i]-1] = abs(nnz_j - i - ((int) (nnz_j / 2)));
+    rst[ord[i]-1] = fabs(((double) nnz_j) - ((double) i) - (((double) nnz_j) / 2.0));
 
   R_Free(ord2);
   R_Free(allord);