{"id":2736,"date":"2026-04-17T17:32:09","date_gmt":"2026-04-17T10:32:09","guid":{"rendered":"https:\/\/bioturing.com\/blog\/?p=2736"},"modified":"2026-04-18T12:23:34","modified_gmt":"2026-04-18T05:23:34","slug":"bioturing-gsea-exact-deterministic-andgpu-accelerated-gene-set-enrichment-analysis","status":"publish","type":"post","link":"https:\/\/bioturing.com\/blog\/bioturing-gsea-exact-deterministic-andgpu-accelerated-gene-set-enrichment-analysis\/","title":{"rendered":"BioTuring-GSEA: Exact, Deterministic, and GPU-Accelerated Gene Set Enrichment Analysis"},"content":{"rendered":"\n

Abstract<\/strong><\/p>\n\n\n\n

Gene Set Enrichment Analysis (GSEA) is a widely adopted method for pathway level interpretation of transcriptomic data. In the canonical formulation of Subra-manian et al. [3], statistical significance is estimated through thousands of phenotype permutations\u2014a procedure that becomes computationally prohibitive for large gene sets and comprehensive pathway databases. Here we present BioTuring-GSEA, a dual-mode framework that addresses this bottleneck through two complementary strategies. The primary mode exploits the formal equivalence between the GSEA statistic at p = 0 and the classical two sample Kolmogorov\u2013Smirnov (KS) statistic [3], enabling fully deterministic, permutation-free p-values via the numerically stable recursive lattice-path algorithm of Viehmann [4], supplemented by the Pelz Good asymptotic series of Simard and L\u2019Ecuyer [2] and Vrbik\u2019s smallsample correction [5]. The second mode retains permutation testing for the general (p>0) case, accelerating it through GPU parallelisation combined with four algorithmic optimisations: massive thread parallelism across independent permutations and gene sets, single shared permutation generation across all gene sets, incremental processing ordered by hit size, and nulldistribution caching for gene sets of equal hit count. Benchmarking on 5,500 gene sets against ranked lists of 50,000 genes demonstrates nearperfect concordance with permutationbased GSEApy (Pearson r \u2265 0.9998) for the deterministic mode, with speedups exceeding 200\u00d7(exact) and 8,000\u00d7(asymptotic); the GPU permutation engine independently achieves 10\u2013100\u00d7speedup over the standard CPU permutation loop.<\/p>\n\n\n\n

Keywords<\/strong>: gene set enrichment analysis; Kolmogorov\u2013Smirnov test; pathway analysis; permutation-free statistics; GPU acceleration; transcriptomics; multiple testing correction.<\/p>\n\n\n\n

1. Introduction<\/strong><\/p>\n\n\n\n

Interpreting genome-wide expression data at the level of biological pathways rather than individual genes is a central challenge in functional genomics. Gene Set Enrichment Analysis (GSEA), introduced by Subramanian et al. [3], determines whether the members of a predefined gene set S tend to cluster towards the top or bottom of a gene list L ranked by phenotype correlation. Unlike threshold-dependent overrepresentation methods, GSEA considers all genes in the experiment, detecting subtle but coordinated expression shifts that would be missed by single-gene analysis [3]. The method has demonstrated broad applicability across cancer biology, metabolic disease, and comparative genomics [3].<\/p>\n\n\n\n

The canonical GSEA algorithm estimates statistical significance by permuting phenotype labels and recomputing the Enrichment Score (ES) thousands of times to construct an empirical null distribution [3]. This permutation strategy preserves gene\u2013gene correlations and is statistically well-motivated; however, it scales poorly. For modern transcriptomic studies profiling >50000 genes against >5000 gene sets simultaneously, the permutation bottleneck renders real-time exploratory analysis impractical.<\/p>\n\n\n\n

A natural optimisation arises for the p = 0 case. As noted explicitly in the original GSEA formulation [3], when p= 0 the ES reduces to the standard two-sample KS statistic, for which an exact combinatorial distribution is known [1]. We exploit this equivalence in a fully analytical implementation. For the general p > 0 case, where permutations remain necessary to preserve gene\u2013gene correlation structure, we describe a complementary GPU-accelerated permutation engine that eliminates the dominant sources of redundancy in the standard implementation. Together, these constitute BioTuring-GSEA, a comprehensive acceleration of GSEA across both its weighted and unweighted formulations<\/p>\n\n\n\n

2. Methods<\/strong><\/p>\n\n\n\n

2.1. The p= 0 Case and KS Equivalence<\/strong><\/p>\n\n\n\n

Let L= {g1,\u2026,gN} be genes ranked in decreasing order of a correlation metric\u2014typically sign(log2 FC) \u00d7(\u2212log10 FDR) for RNA-seq data. Given a gene set S of NH members, the ES running-sum statistic is:<\/p>\n\n\n\n

\n
\"\"<\/figure>\n<\/figure>\n\n\n\n

At p= 0, |rj |p = 1 for all genes and the ES becomes the maximum absolute difference between two empirical CDFs\u2014the ECDF of the n= NH hit positions and the ECDF of the m= N\u2212NH miss positions. This is exactly the two-sample KS statistic [3]:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Hypotheses tested:<\/strong><\/p>\n\n\n\n

\u2022 H0<\/sub>:<\/span> Pathway genes are randomly distributed throughout the ranked list (no association with phenotype)<\/p>\n\n\n\n

\u2022 Ha<\/sub>:<\/span> Pathway genes are enriched at the extremes of the ranked list (associated with phenotype)<\/p>\n\n\n\n

For RNA-seq data ranked by a combined fold-change\u2013significance metric, the p<\/em> = 0 formulation is statistically appropriate: the ranking metric already jointly encodes effect size and confidence, rendering additional per-gene weighting redundant.<\/p>\n\n\n\n

2.2. Exact P-value via Numerically Stable Lattice-Path Counting<\/strong><\/p>\n\n\n\n

Under the null hypothesis H0<\/sub><\/span>, the n hit positions are a uniformly random subset of {1,\u2026,N} and the exact p-value is:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Na\u00efve computation of ( n+m<\/sup>\u2044n<\/sub> ) causes integer overflow for gene-set sizes common in pathway databases (n<\/em> ~ 50\u2013500, m<\/em> ~ 50,000). We therefore adopt Viehmann\u2019s numerically stabilised recursion [4], which evaluates the cumulative lattice-path probability directly as a floating-point quantity bounded in [0, 1]:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

This achieves 13\u201315 significant decimal digits throughout, resolving the subtractive cancellation that afflicts the Durbin matrix method [2].<\/p>\n\n\n\n

2.3. Asymptotic Approximation and Small-Sample Correction<\/strong><\/p>\n\n\n\n

The exact recursion has complexity O<\/em>(n<\/em> \u00d7 m<\/em>) and becomes expensive for n<\/em> \u2273 500. For large gene sets, BioTuring-GSEA employs the Pelz\u2013Good asymptotic series as implemented by Simard and L\u2019Ecuyer [2].<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n
\"\"<\/figure>\n\n\n\n
\"\"<\/figure>\n\n\n\n

This simple argument shift reduces the maximum approximation error from \u2248 7% at n<\/em> = 10 (uncorrected) to < 0.27% [5]. BioTuring-GSEA applies the exact method for n<\/em> < 500 and the Simard\u2013L\u2019Ecuyer asymptotic algorithm otherwise.<\/p>\n\n\n\n

2.4. Deterministic NES and FDR Computation<\/strong><\/p>\n\n\n\n

The NES normalises each ES by the expected null ES to enable cross-pathway comparison. From the moments of the KS distribution [2]:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

\n where 0.8687 = ln(2)\u221a\u03c0\/2<\/span><\/span> \n is the mean of the Kolmogorov limiting distribution. The NES is then \n ES \/ E[ESnull<\/sub>]<\/span>, \n a quantity that is fully deterministic and requires no permutation.\n <\/p>\n\n\n\n

Following [3], the FDR is:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

3. Results<\/strong><\/p>\n\n\n\n

3.1. Concordance with GSEApy: Exact vs. 10000 Permutations<\/strong><\/p>\n\n\n\n

We applied both BioTuring-GSEA (exact method) and GSEApy (10000 permutations as gold standard) to identical ranked gene lists and gene-set databases. Figure 1 shows scatter plots of all four primary output metrics across all evaluated gene sets.<\/p>\n\n\n

\n
\"\"<\/figure>\n<\/div>\n\n\n

Figure 1: Concordance between BioTuring-GSEA (exact method) and GSEApy (10000 permutations) across all primary output metrics<\/strong>. Each panel plots the value reported by GSEApy (y-axis) against the analytically computed value from BioTuring-GSEA (x-axis) for every gene set evaluated. The dashed line indicates the identity (y= x). ES and NES achieve perfect correlation (r= 1.0000); NOM p-value and FDR q-value show near-perfect agreement (r \u22650.9998), with the marginal residual attributable to finite-sample noise inherent to 10000-permutation estimation.<\/p>\n\n\n\n

Table 1 summarises the Pearson correlations. The perfect ES and NES concordance confirms that the analytical normalisation of Eq. 9 is fully equivalent to permutation-derived normalisation. The FDR q-value achieves r= 0.9999, indicating that analytical FDR estimation matches empirical estimation to essentially machine precision.<\/p>\n\n\n\n

Table 1: Pearson correlation between BioTuring-GSEA (exact method) and GSEApy (10000 permutations) for all primary output metrics.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

3.2. Internal Concordance: Exact vs. Asymptotic Method<\/strong><\/p>\n\n\n\n

Figure 2 demonstrates the concordance between the two analytical methods within BioTuring-GSEA\u2014the exact Viehmann recursion and the Simard\u2013L\u2019Ecuyer asymptotic approximation.<\/p>\n\n\n

\n
\"\"<\/figure>\n<\/div>\n\n\n

Figure 2: Internal concordance between BioTuring-GSEA exact and asymptotic methods<\/strong>. All four output metrics achieve r = 1.0000, confirming that the asymptotic Simard\u2013L\u2019Ecuyer algorithm is indistinguishable from the exact Viehmann recursion across the full range of gene-set sizes tested. The NOM p-value and FDR q-value panels, which are most sensitive to numerical differences, show perfect point-on-line agreement, validating the method-switching
boundary at n = 500.<\/p>\n\n\n\n

The perfect concordance across all metrics validates the decision to switch from the exact recursion to the Pelz\u2013Good asymptotic series at n = 500, and confirms that the Vrbik small-sample correction (Eq. 8) successfully bridges the transition region.<\/p>\n\n\n\n

3.3. Effect of Permutation Count on Convergence<\/p>\n\n\n\n

A key advantage of the analytical approach is immunity to the stochastic variance that affects permutation-based methods at finite permutation counts. Figure 3 illustrates this directly by
comparing BioTuring-GSEA NOM values against GSEApy run with 1000 permutations (left) and 10000 permutations (centre), alongside the exact-vs-asymptotic internal comparison (right).<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Figure 3: NOM p-value concordance as a function of permutation count<\/strong>. Left<\/em>: GSEApy with 1000 permutations vs. BioTuring-GSEA exact method (r= 0.9981). The visible scatter reflects the irreducible stochastic variance of finite permutation estimation. Centre: GSEApy with 10000 permutations (r= 0.9998); increasing permutation count reduces variance but does not eliminate it. Right<\/em>: BioTuring-GSEA asymptotic vs. exact method (r = 1.0000); the analytical approaches agree perfectly, confirming that all residual discordance in the left and centre panels originates from permutation sampling error rather than any deficiency of the exact computation.<\/p>\n\n\n\n

This comparison establishes an important asymmetry: permutation-based methods converge towards the analytical result as permutation count increases, but even 10000 permutations retain non-negligible variance (particularly in the p > 0.5 region). The analytical method provides the limiting value of this convergence with zero additional computational overhead for p-value estimation.<\/p>\n\n\n\n

3.4. Computational Performance<\/strong><\/p>\n\n\n\n

Figure 4 shows the runtime ratio (GSEApy \/ BioTuring-GSEA) as a function of permutation count, benchmarked on 50000 genes against 5500 gene sets.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Figure 4: Runtime speedup of BioTuring-GSEA relative to GSEApy as a function of permutation count<\/strong>. The asymptotic method (red) scales linearly with permutation count, reaching >8,000\u00d7speedup at 10000 permutations. The exact method (purple) is independent of permutation count and achieves a stable >200\u00d7speedup throughout. The grey dashed line marks the break-even ratio of 1. Benchmark: 50000 genes \u00d75500 gene sets.<\/p>\n\n\n\n

The runtime advantage of the asymptotic method scales linearly with permutation count because GSEApy\u2019s cost grows proportionally to the number of permutations while BioTuring-GSEA\u2019s cost is fixed. At the standard 1000-permutation setting, the asymptotic method is already \u223c1,000\u00d7faster; at 10000 permutations the ratio exceeds 8,000\u00d7.<\/p>\n\n\n\n

Table 2: Computational speedup of BioTuring-GSEA relative to GSEApy at 10000 permutations. Benchmark: 50000 genes \u00d75500 gene sets.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

4. GPU-Accelerated Permutation Engine for p<\/em> > 0<\/strong><\/p>\n\n\n\n

The results in Section 3 demonstrate that for the p = 0 case, the analytical KS approach provides exact, deterministic results with massive computational advantages over permutation based methods. However, many GSEA applications require the weighted formulation (p > 0) to properly account for gene-level effect sizes or continuous phenotype correlations. In these scenarios, permutation testing remains necessary to preserve the gene\u2013gene correlation structure under the null hypothesis [3].<\/p>\n\n\n\n

The standard implementation of permutation-based GSEA exhibits severe computational redundancies that can be eliminated through careful algorithmic design. BioTuring-GSEA implements a GPU-accelerated permutation engine that preserves the exact mathematics of the canonical GSEA permutation test [3] while achieving substantial speedup over the standard implementation (e.g., GSEApy). The acceleration rests on four complementary strategies:<\/p>\n\n\n\n

    \n
  1. Massive GPU thread parallelism across independent permutations and gene sets<\/li>\n\n\n\n
  2. Single shared permutation generation for all gene sets<\/li>\n\n\n\n
  3. Hit-size sorting with incremental running-sum processing<\/li>\n\n\n\n
  4. Null-distribution caching for gene sets of identical hit count<\/li>\n<\/ol>\n\n\n\n

    These optimisations collectively deliver their greatest performance gains at large permutation counts (P<\/em> \u2265 1,000) and large pathway databases (G<\/em> \u2265 1,000)\u2014precisely the regime where permutation testing becomes computationally prohibitive in standard implementations.<\/p>\n\n\n\n

    4.1. Natural Parallelism in GSEA Permutation Testing<\/p>\n\n\n\n

    The canonical permutation loop [3] possesses two structural properties that map perfectly onto GPU architectures:<\/p>\n\n\n\n

    Independence across permutations.<\/strong> The P phenotype permutations are statistically independent: the outcome of permutation k<\/em> does not depend on permutation k<\/em> \u22121. Each permutation involves reordering the gene list, recomputing correlation scores, and recalculating the ES for all
    gene sets\u2014operations that can be executed concurrently without communication.<\/p>\n\n\n\n

    Independence across gene sets. Given a fixed ranked list L<\/em> and correlation scores rj <\/em>from a particular permutation, gene sets are mutually independent: the ES of gene set A does not affect the ES of gene set B<\/em>. This independence holds both for the observed data and for each null permutation.<\/p>\n\n\n\n

    Consequently, the full P\u00d7G<\/em> permutation matrix (where G<\/em> is the number of gene sets) constitutes an embarrassingly parallel workload that can be dispatched in a single GPU kernel. Standard
    implementations serialise this computation through nested loops, leaving the vast majority of available parallelism unexploited.<\/p>\n\n\n\n

    4.2. Shared Permutation Generation<\/strong><\/p>\n\n\n\n

    Standard implementations generate a fresh set of P random hit indices independently for every gene set. For a typical database (G\u22485,500, P = 10,000), this incurs G\u00d7P = 55 million separate random-number generation calls and memory allocations, representing a severe bottleneck in both computation time and memory bandwidth.<\/p>\n\n\n\n

    \"\"<\/figure>\n\n\n\n