Add max_size parameter to fast_ssgsea

Gene sets with more than max_size genes will not be analyzed. Decreasing
max_size will decrease the overall runtime proportionally, though this
is less noticeable when the number of permutations is small.

Author: TylerSagendorf

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: fast.ssgsea
 Type: Package
 Title: Fast Single-Sample Gene Set Enrichment Analysis (ssGSEA)
-Version: 0.1.0.9020
+Version: 0.1.0.9021
 Date: 2025-11-11
 Authors@R: 
       person(given = "Tyler", family = "Sagendorf", 

NEWS.md

@@ -4,6 +4,7 @@
 - Greatly reduced runtime of permutation tests.
 - Added `read_gmt` function, which reads a list of gene sets from a Gene Matrix Transposed (GMT) file.
 - Added `alternative` parameter to `fast_ssgsea` to perform one-sided hypothesis tests.
+- Added `max_size` parameter to `fast_ssgsea` to limit the maximum size of sets that will be tested.
 - Updated runtime data and figures in simulation/.
 
 

R/fast_ssgsea.R

@@ -26,6 +26,8 @@
 #'   \emph{all} columns of \code{X}. The default value of 2 is the minimum
 #'   possible set size required for testing, though higher values tend to
 #'   produce more robust results.
+#' @param max_size integer or \code{Inf}; the maximum set size. Recommended to
+#'   set this to 500.
 #' @param sort logical; should the results for each column of \code{X} be sorted
 #'   by p-value? Default is \code{TRUE}.
 #' @param seed integer or \code{NULL}; passed to \code{\link[base]{set.seed}}.
@@ -117,6 +119,7 @@ fast_ssgsea <- function(X,
                         batch_size = 1000L,
                         adjust_globally = FALSE,
                         min_size = 2L,
+                        max_size = Inf,
                         sort = TRUE,
                         seed = NULL,
                         alternative = c("two.sided", "less", "greater")) {
@@ -130,6 +133,7 @@ fast_ssgsea <- function(X,
     batch_size = batch_size,
     adjust_globally = adjust_globally,
     min_size = min_size,
+    max_size = max_size,
     sort = sort,
     seed = seed,
     n_genes = ncol(X)
@@ -171,10 +175,11 @@ fast_ssgsea <- function(X,
     Z_prime = Z[, rownames(A), drop = FALSE], # Z'
     A = A,
     A_d = A_d,
-    min_size = min_size
+    min_size = min_size,
+    max_size = max_size
   )
 
-  # Extract list components: M, W, M_d, W_d, A (optional), A_d (optional)
+  # Extract list components: M, W, M_d, W_d, A, A_d, max_size
   list2env(x = M_list, envir = environment())
 
   # Enrichment score matrices with samples as rows and gene sets as columns

R/utils.R

@@ -52,6 +52,7 @@
                             batch_size = 1000L,
                             adjust_globally = FALSE,
                             min_size = 2L,
+                            max_size = Inf,
                             sort = TRUE,
                             seed = NULL,
                             n_genes) {
@@ -103,6 +104,16 @@
     stop("`min_size` must be >= 2 and < nrow(X).")
   }
 
+  if (
+    !(is.infinite(max_size) || is.vector(max_size, mode = "numeric")) ||
+      length(max_size) != 1L ||
+      is.na(max_size) ||
+      max_size < min_size ||
+      (!is.infinite(max_size) && max_size %% 1 != 0)
+  ) {
+    stop("`max_size` must be Inf or an integer >= min_size.")
+  }
+
   if (
     !is.vector(adjust_globally, mode = "logical") ||
       length(adjust_globally) != 1L ||
@@ -319,6 +330,8 @@
 #'   Indicates which genes are expected to be down-regulated in each sample.
 #' @param min_size integer; minimum gene set size required for testing. Default
 #'   is 2.
+#' @param max_size integer or \code{Inf}; size of the largest gene set that will
+#'   be tested.
 #'
 #' @returns A named list with the following components:
 #'
@@ -350,39 +363,44 @@
                          Z_prime,
                          A,
                          A_d = NULL,
-                         min_size = 2L) {
+                         min_size = 2L,
+                         max_size = Inf) {
   # Matrix with samples as rows and genes as columns. Elements are the number
   # of genes in each set with nonmissing values in the sample.
   M <- .Cpp_matmult_sparse(Z_prime, A) # Z'A
   M[M < min_size] <- 0L
 
-  # Identify sets with too few or too many genes in at least one sample
-  extreme_sets_M <- apply(M, 2L, function(m_j) {
-    any(m_j == 0L | m_j == n)
-  })
+  # It is extremely unlikely that a gene set would consist of all genes with
+  # nonmissing values, but the set size limit will be 1 less than the smallest
+  # number of nonmissing values.
+  max_size <- max(min_size, min(max_size, min(n) - 1L))
 
   if (!is.null(A_d)) {
     M_d <- .Cpp_matmult_sparse(Z_prime, A_d)
     M_d[M_d < min_size] <- 0L
 
-    extreme_sets_M_d <- apply(M_d, 2L, function(m_d_j) {
-      any(m_d_j == 0L | m_d_j == n)
-    })
+    small_sets_M <- apply(M == 0L, 2L, any)
+    small_sets_M_d <- apply(M_d == 0L, 2L, any)
+
+    # Combined set size exceeds max_size
+    large_sets <- apply((M + M_d) > max_size, 2L, any)
 
-    extreme_sets <- which(extreme_sets_M & extreme_sets_M_d)
+    # Too small in both up and down portions or too large overall
+    extreme_sets <- which((small_sets_M & small_sets_M_d) | large_sets)
   } else {
     M_d <- W_d <- NULL
 
-    extreme_sets <- which(extreme_sets_M)
+    extreme_sets <- which(
+      apply(M == 0 | M > max_size, 2L, any)
+    )
   }
 
   if (length(extreme_sets)) {
     # If any sets are extreme, check if they are all extreme.
     if (length(extreme_sets) == ncol(A)) {
       stop(
-        "All sets in `gene_sets` contain fewer than `min_size` genes ",
-        "with nonmissing values or consist of all genes with nonmissing ",
-        "values in at least one sample."
+        "All sets in `gene_sets` contain fewer than `min_size` genes or more ",
+        "than `max_size` genes with nonmissing values."
       )
     }
 

README.Rmd

@@ -107,8 +107,6 @@ system.time({
   )
 })
 
-str(res)
-
 head(res, 10L)
 ```
 

README.md

@@ -130,22 +130,7 @@ system.time({
 ```
 
     ##    user  system elapsed 
-    ##   5.919   0.906   6.403
-
-``` r
-str(res)
-```
-
-    ## 'data.frame':    20000 obs. of  9 variables:
-    ##  $ sample      : Factor w/ 1 level "sample1": 1 1 1 1 1 1 1 1 1 1 ...
-    ##  $ set         : chr  "set18791" "set16136" "set19084" "set2830" ...
-    ##  $ set_size    : int  138 801 841 163 706 749 450 87 161 761 ...
-    ##  $ ES          : num  -1866 709 698 1584 759 ...
-    ##  $ NES         : num  -5.3 4.65 4.68 4.76 4.68 ...
-    ##  $ n_same_sign : int  49042 52788 52782 50951 52785 47193 47813 50722 48979 47243 ...
-    ##  $ n_as_extreme: int  1 8 8 9 11 10 13 14 18 20 ...
-    ##  $ p_value     : num  4.08e-05 1.70e-04 1.71e-04 1.96e-04 2.27e-04 ...
-    ##  $ adj_p_value : num  0.739 0.739 0.739 0.739 0.739 ...
+    ##   4.925   0.901   5.370
 
 ``` r
 head(res, 10L)
@@ -192,16 +177,16 @@ print(sessionInfo(), locale = FALSE, tzone = FALSE)
     ## [1] stats     graphics  grDevices utils     datasets  methods   base     
     ## 
     ## other attached packages:
-    ## [1] fast.ssgsea_0.1.0.9018
+    ## [1] dqrng_0.4.1            fast.ssgsea_0.1.0.9021
     ## 
     ## loaded via a namespace (and not attached):
-    ##  [1] dqrng_0.4.1            digest_0.6.37          RcppArmadillo_15.0.2-2
-    ##  [4] fastmap_1.2.0          xfun_0.53              Matrix_1.7-4          
-    ##  [7] lattice_0.22-7         knitr_1.50             htmltools_0.5.8.1     
-    ## [10] rmarkdown_2.29         cli_3.6.5              grid_4.5.2            
-    ## [13] data.table_1.17.8      compiler_4.5.2         rstudioapi_0.17.1     
-    ## [16] tools_4.5.2            evaluate_1.0.5         Rcpp_1.1.0            
-    ## [19] yaml_2.3.10            rlang_1.1.6
+    ##  [1] digest_0.6.37          RcppArmadillo_15.0.2-2 fastmap_1.2.0         
+    ##  [4] xfun_0.53              Matrix_1.7-4           lattice_0.22-7        
+    ##  [7] knitr_1.50             htmltools_0.5.8.1      rmarkdown_2.29        
+    ## [10] cli_3.6.5              grid_4.5.2             data.table_1.17.8     
+    ## [13] compiler_4.5.2         rstudioapi_0.17.1      tools_4.5.2           
+    ## [16] evaluate_1.0.5         Rcpp_1.1.0             yaml_2.3.10           
+    ## [19] rlang_1.1.6
 
 ## Performance
 

man/fast_ssgsea.Rd

@@ -201,9 +201,8 @@ test_that("extreme sets are removed", {
   expect_identical(
     object = err,
     expected = paste0(
-      "All sets in `gene_sets` contain fewer than `min_size` genes with ",
-      "nonmissing values or consist of all genes with nonmissing values ",
-      "in at least one sample."
+      "All sets in `gene_sets` contain fewer than `min_size` genes or ",
+      "more than `max_size` genes with nonmissing values."
     )
   )
 })
@@ -243,9 +242,8 @@ test_that("extreme directional sets are removed", {
   expect_identical(
     object = err,
     expected = paste0(
-      "All sets in `gene_sets` contain fewer than `min_size` genes with ",
-      "nonmissing values or consist of all genes with nonmissing values ",
-      "in at least one sample."
+      "All sets in `gene_sets` contain fewer than `min_size` genes or ",
+      "more than `max_size` genes with nonmissing values."
     )
   )
 })