Greatly speed up calculation of permutation ES

Switched to a method that eliminates the need to compute matrix dot
products when deriving permutation ES. This new version can generate
results faster than the previous version when R was linked with the
OpenBLAS library, and the new approach also requires less memory.
Unfortunately, this change breaks support for directional databases,
such as PTMsigDB.

Author: TylerSagendorf

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: fast.ssgsea
 Type: Package
 Title: Fast Single-Sample Gene Set Enrichment Analysis (ssGSEA)
-Version: 0.1.0.9017
-Date: 2025-09-14
+Version: 0.1.0.9018
+Date: 2025-10-30
 Authors@R: 
       person(given = "Tyler", family = "Sagendorf", 
              email = "tyler.sagendorf@pnnl.gov", 

R/RcppExports.R

@@ -1,6 +1,36 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' @title Multiplication of an unseen binary matrix and real-valued matrix
+#'
+#' @description Compute the product of an unseen binary and real-valued matrix
+#'   when the binary matrix is in a specific format. See details.
+#'
+#' @param set_sizes integer vector of unique set sizes, sorted in ascending
+#'   order.
+#' @param Y real-valued matrix. The number of rows is equal to
+#'   \code{max(set_sizes)}.
+#'
+#' @details The set_sizes vector acts as a stand-in for the unseen binary
+#'   matrix, which has the following characteristics:
+#'
+#'   1. The first `set_sizes[i]` elements of the i-th row are all 1 and the
+#'   remaining entries of that row are 0.
+#'
+#'   2. If the i-th row has n ones, then the i + 1 row will have at least (n +
+#'   1) ones.
+#'
+#'   3. All elements of the last row are 1.
+#'
+#' @author Tyler Sagendorf
+#'
+#' @references Sanderson, C., & Curtin, R. (2016). Armadillo: A template-based
+#'   C++ library for linear algebra. The Journal of Open Source Software, 1(2),
+#'   26. \url{https://doi.org/10.21105/joss.00026}
+#'
+#' @noRd
+NULL
+
 #' @title Find Index of First Positive Value in Sorted Vector
 #'
 #' @description Given a sorted vector of real-valued numbers, find the index of
@@ -170,8 +200,8 @@ NULL
 #'
 #' @noRd
 #'
-.Rcpp_calcESPermCore <- function(alpha, Y_perm, R_perm, sumRanks_i, A_perm, theta_m_i, theta_w_i) {
-    .Call(`_fast_ssgsea_Rcpp_calcESPermCore`, alpha, Y_perm, R_perm, sumRanks_i, A_perm, theta_m_i, theta_w_i)
+.Rcpp_calcESPermCore <- function(alpha, Y_perm, R_perm, sumRanks_i, theta_m_i, theta_w_i) {
+    .Call(`_fast_ssgsea_Rcpp_calcESPermCore`, alpha, Y_perm, R_perm, sumRanks_i, theta_m_i, theta_w_i)
 }
 
 #' @title Extract Information from a Permutation Enrichment Score Matrix

R/globals.R

@@ -1,11 +1,15 @@
 utils::globalVariables(
   c(
+    "A",
+    "A_d",
     "adj_p_value",
     "A_perm",
     "A_perm_d",
     "direction_down",
     "ES",
+    "ES_d",
     "ES_idx",
+    "ES_u",
     "i",
     "j",
     "M",

R/utils.R

@@ -213,7 +213,16 @@
   dims <- lengths(dim_names)
 
   # Keep genes expected to be "up"
-  dt_u <- dt[direction_down == FALSE]
+  if (any(dt[["direction_down"]])) {
+    stop(
+      "fast_ssgsea does not currently support directional databases. ",
+      "This will be fixed in a future update."
+    )
+
+    dt_u <- dt[direction_down == FALSE]
+  } else {
+    dt_u <- dt
+  }
 
   # Incidence matrix where a 1 indicates that the element is in the set. If x
   # is a directional database, then A will only contain elements that are
@@ -533,13 +542,16 @@
                         y_i,
                         r_i,
                         sumRanks_i,
-                        A_perm,
                         theta_m_i,
                         theta_w_i,
-                        A_perm_d = NULL,
                         theta_m_d_i = NULL,
                         theta_w_d_i = NULL) {
-  max_set_size <- ncol(A_perm)
+  # max_set_size <- ncol(A_perm)
+  if (is.null(theta_m_d_i)) {
+    max_set_size <- max(theta_m_i)
+  } else {
+    max_set_size <- max(theta_m_i + theta_m_d_i)
+  }
 
   n_elements <- length(element_indices) # number of nonmissing values
 
@@ -571,18 +583,16 @@
     Y_perm,
     R_perm,
     sumRanks_i,
-    A_perm,
     theta_m_i,
     theta_w_i
   )
 
-  if (!is.null(A_perm_d)) { # directional sets
+  if (!is.null(theta_m_d_i)) { # directional sets
     ES_perm_d <- .Rcpp_calcESPermCore(
       alpha,
       Y_perm,
       R_perm,
       sumRanks_i,
-      A_perm_d,
       theta_m_d_i,
       theta_w_d_i
     )
@@ -787,22 +797,22 @@
 
     A_perm_d <- NULL
 
-    A_perm <- .Rcpp_calcAPerm(
-      end = theta_m_i,
-      MAX_SET_SIZE = max_set_size,
-      check = FALSE
-    )
+    # A_perm <- .Rcpp_calcAPerm(
+    #   end = theta_m_i,
+    #   MAX_SET_SIZE = max_set_size,
+    #   check = FALSE
+    # )
 
     rep_idx <- match(x = m_i, table = theta_m_i)
   }
 
   out <- list(
     "rep_idx" = rep_idx,
-    "A_perm" = A_perm,
+    # "A_perm" = A_perm,
     "theta_m_i" = theta_m_i,
     "theta_w_i" = theta_w_i,
     # These may be NULL
-    "A_perm_d" = A_perm_d,
+    # "A_perm_d" = A_perm_d,
     "theta_m_d_i" = theta_m_d_i,
     "theta_w_d_i" = theta_w_d_i
   )
@@ -895,6 +905,7 @@
   tab_i <- data.table(
     set = sets,
     set_size = m_i,
+    m_d_i = m_d_i,
     ES = ES_i,
     # Initialize vectors of 0's. These 3 vectors will be updated using the
     # results from each batch of permutations.
@@ -905,6 +916,16 @@
     stringsAsFactors = FALSE
   )
 
+  setorderv(
+    x = tab_i,
+    cols = c("set_size", "ES", "row_order"),
+    order = c(1L, 1L, 1L)
+  )
+
+  # Reordered
+  m_i <- tab_i[["set_size"]]
+  m_d_i <- tab_i[["m_d_i"]]
+
   if (nperm != 0L) {
     # Incidence matrices and other information for permutations
     A_list_perm <- .permIncidenceMatrix(
@@ -913,14 +934,12 @@
       m_d_i = m_d_i
     )
 
-    # Extract list components: rep_idx, A_perm, theta_m_i, theta_w_i,
-    # A_perm_d, theta_m_d_i, theta_w_d_i
+    # Extract list components: rep_idx, theta_m_i, theta_w_i, theta_m_d_i,
+    # theta_w_d_i
     list2env(x = A_list_perm, envir = environment())
 
     tab_i[, rep_idx := rep_idx]
 
-    setorderv(x = tab_i, cols = c("rep_idx", "ES"), order = c(1L, 1L))
-
     ES_ls <- split(
       x = tab_i[["ES"]],
       f = tab_i[["rep_idx"]]
@@ -941,10 +960,8 @@
         y_i = y_i,
         r_i = r_i,
         sumRanks_i = sumRanks_i,
-        A_perm = A_perm,
         theta_m_i = theta_m_i,
         theta_w_i = theta_w_i,
-        A_perm_d = A_perm_d,
         theta_m_d_i = theta_m_d_i,
         theta_w_d_i = theta_w_d_i
       )
@@ -954,7 +971,7 @@
         ES_perm = ES_perm
       )
 
-      perm_dt <- rbindlist(perm_dt)
+      perm_dt <- rbindlist(perm_dt, use.names = FALSE)
 
       # Update summary vectors
       tab_i[, `:=`(

README.Rmd

@@ -75,8 +75,8 @@ X <- matrix(
 
 ## Simulate list of gene sets
 n_sets <- 20000L # number of gene sets
-min_size <- 10L # size of smallest gene set
-max_size <- 500L # size of largest gene set
+min_size <- 5L # size of smallest gene set
+max_size <- 1000L # size of largest gene set
 
 size_range <- max_size - min_size + 1L
 n_reps <- ceiling(n_sets / size_range)
@@ -91,7 +91,7 @@ names(gene_sets) <- paste0("set", seq_along(gene_sets))
 
 ### Runtime and Results
 
-This shows the runtime of `fast_ssgsea` with the reference Basic Linear Algebra Subprograms (BLAS) library [@blas] (single-threaded) running on an AMD Ryzen 5 7600X CPU with a clock speed of 4.7 GHz.
+This shows the runtime of `fast_ssgsea` running on an AMD Ryzen 5 7600X CPU with a clock speed of 4.7 GHz.
 
 ```{r time-results}
 library(fast.ssgsea)
@@ -123,78 +123,17 @@ print(sessionInfo(), locale = FALSE, tzone = FALSE)
 
 ## Performance
 
-The `fast.ssgsea` R package utilizes linear algebra and ideas from Fast Gene Set Enrichment Analysis [@korotkevich-fast-2021] to greatly reduce the runtime of ssGSEA and PTM-SEA while also properly controlling the type I error rate.
+The `fast.ssgsea` R package utilizes linear algebra and ideas from Fast Gene Set Enrichment Analysis [@korotkevich-fast-2021] to greatly reduce the runtime of gene permutation GSEA and PTM-SEA.
 
 Tests were performed on a desktop computer with an AMD Ryzen 5 7600X CPU (6 cores, 12 threads) at 4.7 GHz. Different combinations of the number of samples, gene sets, maximum gene set size, number of permutations, and value of the $\alpha$ parameter (the weighting exponent) were tested in a random order (3 replicates each) to minimize the influence of previous runs.
 
 ```{r, echo=FALSE}
-fig1_cap <- "Runtime of fast_ssgsea with A) 1,000 or B) 10,000 permutations. R was linked to the reference BLAS library, so only a single thread was used."
-
-fig2_cap <- "Runtime of fast_ssgsea with A) 1,000 or B) 10,000 permutations. R was linked to OpenBLAS, and all 12 threads were used."
+fig_cap <- "Runtime of fast_ssgsea with A) 1,000 or B) 10,000 permutations."
 ```
 
 
-```{r, echo=FALSE, fig.cap=fig1_cap}
+```{r, echo=FALSE, fig.cap=fig_cap}
 knitr::include_graphics("./man/figures/README-figure-1.png")
 ```
 
-## Switching the BLAS Library
-
-Linking R to an external BLAS, such as the optimized, open-source OpenBLAS library [@openblas-1; @openblas-2], can greatly reduce the runtime (though it may affect precision):
-
-```{r, echo=FALSE, fig.cap=fig2_cap}
-knitr::include_graphics("./man/figures/README-figure-2.png")
-```
-
-### Linux
-
-Linux users can follow [these instructions](https://docs.posit.co/resources/install-r-source.html#optional-configure-r-to-use-a-different-blas-library) from Posit to easily switch the BLAS library.
-
-### macOS
-
-For Mac computers with an Apple Silicon chip (M1 and beyond), follow [these instructions](https://cran.r-project.org/bin/macosx/RMacOSX-FAQ.html#Which-BLAS-is-used-and-how-can-it-be-changed_003f) to switch to Apple's optimized BLAS library.
-
-### Windows
-
-To install OpenBLAS for R on Windows, users may follow [this tutorial](https://github.com/david-cortes/R-openblas-in-windows).
-
-
-## Further Improvements
-
-If using an external BLAS library, the runtime can be reduced by ~1/3 by switching from double-precision to single-precision floating point arithmetic for the permutation tests. If the BLAS supports multi-threading, the difference in runtime will be negligible, so the switch is largely unnecessary. While using floats will slightly affect the precision of the normalized enrichment scores (NES), the differences are small compared to differences observed from changing the value of the `seed` parameter. Unfortunately, the R BLAS does not support floats (see [this issue](https://github.com/RcppCore/RcppArmadillo/issues/197) in RcppArmadillo). Since Windows has no default BLAS/LAPACK library, it was not possible to implement this in `fast.ssgsea` without complicating the installation process or making it impossible for Windows users.
-
-For users looking to implement this change, please follow these instructions:
-
-1. Link R to an external BLAS library, such as OpenBLAS. Verify success by examining the result of `sessionInfo()["BLAS"]`.
-2. Clone pnnl/fast.ssgsea (e.g., with `git clone https://github.com/pnnl/fast.ssgsea` in a terminal) and open the fast.ssgsea.Rproj file.
-3. In src/Rcpp_functions.cpp, replace the `Rcpp_calcESPermCore` C++ function with the following:
-
-```
-arma::fmat Rcpp_calcESPermCore(const float alpha,
-                               const arma::fmat& Y_perm,
-                               const arma::fmat& R_perm,
-                               const float sumRanks_i,
-                               const arma::fmat& A_perm,
-                               const arma::fvec& theta_m_i,
-                               const arma::fvec& theta_w_i)
-{
-   arma::fmat AR_perm = A_perm * R_perm;
-
-   arma::fmat ES_perm(A_perm.n_rows, Y_perm.n_cols, arma::fill::zeros);
-
-   if (alpha == 0.0f) {
-      ES_perm = arma::diagmat(1.0f / theta_m_i) * AR_perm;
-   } else {
-      ES_perm = (A_perm * (Y_perm % R_perm)) / (A_perm * Y_perm);
-   }
-
-   ES_perm += arma::diagmat(1.0f / theta_w_i) * (AR_perm - sumRanks_i);
-
-   return ES_perm;
-}
-```
-
-4. Build and install `fast.ssgea`.
-
-
 ## References