Gene-set inference: calibrating spiGSEA and why competitive is the default

Dharmesh D. Bhuva

What is being tested

spiGSEA() scores a gene set by averaging its member genes’ statistics: each per-gene $t$ is converted to a $z$, the set’s $m$ z-scores are averaged, and the average is standardised by

$$\sqrt{\frac{1 + \rho\,(m - 1)}{m}},$$

the variance of a mean of $m$ equicorrelated unit variates — $\rho$ is the average inter-gene correlation, estimated from the counts (or reused from the fit). Without the $\rho$ term the divisor is $1/\sqrt m$, which overstates the precision of the mean for any positively correlated set and makes every set-level p-value too small.

The function offers two null hypotheses, selected by test:

• self-contained — is the set’s mean statistic different from zero?
• competitive (the default) — is the set’s mean different from the genes outside the set, standardised by the observed spread of the gene-level statistics (the limma::camera construction)?

They answer different questions, and this study measures which is safe to default to. The headline is stark enough to state up front: on a fully null benchmark at realistic inter-gene correlation, the self-contained test calls about 20 of 208 sets per replicate — all false — while the competitive test essentially never calls; and the failure is only partly the correlation term’s fault, which is why no tuning of $\rho$ rescues it.

Methods

Design. The gene-set layer has its own generative grid, separate from the gene-level simulation study: 1,500 genes (150 markers per cell type, 450 housekeeping), $S = 30$ samples of 300 cells, 500 replicates per design point. Candidate sets number 208 in total: at each size $m \in \{10, 25, 50, 100\}$, 25 null sets (genes drawn at random), 25 high sets (drawn from the top expression quintile), and 2 plant sets (which carry signal only in the power arm).

Correlation levels. The base count model draws genes independently, so inter-gene correlation is injected via a latent factor loading and then measured (next section). Three levels are contrasted: none (loading 0, no library-size variation — genuinely independent genes), libonly (loading 0, library-size variation on — which alone leaves $\hat\rho \approx 0.0036$), and real (the calibrated loading, giving the $\hat\rho \approx 0.011$ observed on real spatial data).

Arms. null plants nothing anywhere; power plants an effect in the designated plant sets only, over effect sizes $\beta \in \{0.3, 0.6, 1.2, 2\}$; global shifts a broad share of all genes with no set-specific signal — the scenario that discriminates the two nulls, since a competitive test should ignore it and a self-contained test should not. Every arm runs both tests, and each is also re-scored with the correlation term switched off (rho_mode = "zero") to separate the $\rho$ term’s contribution from everything else.

Scoring. A set’s p-value is the unweighted Cauchy combination over its triplets — level 1 of spiGSEA()’s own cascade — with BH across sets; type-I is read from the per-set p at $\alpha$, call rates from BH at the stated FDR. On the null arm every call is false by construction.

Calibrating the injected correlation

The set statistic’s variance term is driven entirely by $\rho$, so before the study can claim “realistic correlation” the loading has to be mapped to the correlation it actually produces, with the package’s own estimator (spiDE:::.interGeneCor()) doing the measuring.

Realised average inter-gene correlation (as spiDE estimates it) against the injected latent-factor loading. The dashed line is the target measured on real CosMx data; the loading 0.30 reproduces it. At loading 0 the estimate does not return to zero: per-cell library-size variation alone leaves rho at about 0.0036.

The mapping is smooth and monotone; the calibrated loading (0.30) realises $\hat\rho = 0.011$, matching the real-data target, and the loading-0 floor of $\hat\rho = 0.0036$ shows that per-cell library-size variation is itself a correlation source — which is why the “independent genes” contrast in this study switches both the loading and the library-size variation off, not just the loading.

Null calibration: the two tests are not interchangeable

False call rate at a nominal 5% FDR with nothing planted, at realistic inter-gene correlation. Every call shown is false. The self-contained test does not control error, and its errors concentrate in high-expression sets.

The dashed line is the nominal 5 %. The competitive test sits at zero; the self-contained test calls 2.7 % of random sets and 17.5 % of high-expression sets — about 20.6 of the 208 sets per replicate, all false, so its realised false-discovery proportion is 1.00 against a reported 0.05 (the competitive test’s mean is 0.05 calls per replicate). The same ordering holds before any FDR correction: per-set type-I at $\alpha = 0.05$ is 0.082 / 0.36 (self-contained, random / high sets) against 0.013 / 0.067 (competitive).

Null QQ curves of the per-set p-values at realistic correlation, by candidate-set class. Above the diagonal is anti-conservative. The self-contained curve lifts off most strongly for high-expression sets; the competitive curve stays at or below the diagonal for random sets with a mild lift on high-expression sets.

Why the self-contained test fails — two causes, one of them untunable

It is worth being precise about the mechanism, because the obvious suspect — the correlation term — is only half the story.

• The correlation term does real work. Re-scoring the same fits with the $\rho$ term switched off (rho_mode = "zero") raises the self-contained false-call rate on random sets from 2.7 % to 11.8 % — roughly a 4.4-fold increase. Estimating $\rho$ and inflating the variance is necessary.
• It is not sufficient. The self-contained test remains badly anti-conservative with the correct $\rho$, because it also assumes the gene-level z-scores it averages have unit spread. That holds under the null (measured spread 1.01) and fails under signal — the spread reaches 1.82 at the largest planted effect, and 1.4 even with genuinely independent genes, so no choice of $\rho$ repairs it. The competitive test divides by the observed spread, which is why it is immune.
• The control case closes the argument. With genuinely independent genes (loading and library-size variation both off) both tests are calibrated (self-contained false-call rate 0.2 % on random sets at FDR 0.05) — identifying the correlation structure, not the model or the FDR cascade, as the trigger.

The discriminating scenario: a global shift

A coordinated shift across most genes — a global response with no set-specific biology — is exactly where the two null hypotheses part ways: a competitive test asks “is this set different from the rest?” (no), a self-contained test asks “is this set’s mean nonzero?” (yes, like everything else’s).

False call rate on random sets at 5% FDR when a global shift is planted with no set-specific signal, versus the shift size. The self-contained test calls nearly everything at the larger shift; the competitive test correctly ignores it.

At the larger shift the self-contained test calls 96.2 % of random sets against 0.7 % for competitive. On real tissue, where broad transcriptional shifts between conditions are the rule rather than the exception, this failure mode is not hypothetical.

Where the false calls land: expression bias

Because the per-gene statistics track expression level, the self-contained test’s false calls are not spread evenly: under the null its called sets sit at the 86th expression percentile against a 69th-percentile average over all candidate sets — the same “your top pathways are just the abundant ones” pattern seen on real data, here reproduced with known ground truth.

Power of the competitive test

Strict error control is only worth having if the test still detects real set-level signal. On the power arm the effect lives in the plant sets only.

Detection rate of the planted sets at 5% FDR versus effect size, by set size, at realistic correlation. The self-contained curves are shown for completeness but sit on an invalid null (see above); the competitive curves are the usable ones.

Competitive power is usable from around 25 genes per set at moderate effects (detection 0.8 at $m = 25$, $\beta = 1.2$; 0.3 at $m = 50$, $\beta = 0.6$) and saturates by $m = 50$ at the larger effects. The self-contained curves sit above it everywhere, but that advantage is priced in false discoveries: a test that calls 10 % of null sets buys its extra “power” from the same broken null that produced the calibration figures above.

Cost

The gene-set layer is computationally free at this scale: on top of a 261-second model fit, estimating $\rho$ costs about 1.2 s and scoring all 208 sets about 1 s.

Summary

• The self-contained test does not control its error rate on realistic spatial data: FDP 1.00 at a reported 0.05 on the null benchmark, with false calls concentrated in high-expression sets, and near-total failure under a global shift. Two causes: inter-gene correlation (tunable, via the $\rho$ inflation) and the unit-spread assumption on the averaged z-scores (untunable — it breaks under any real signal).
• The competitive test is calibrated at every correlation level and scenario measured, because it standardises by the observed spread and contrasts against the out-of-set genes; its power is usable from about 25 genes per set.
• spiGSEA() therefore defaults to test = "competitive". The self-contained form is retained only to reproduce the original flat script’s fry-based results, and should not be used for inference.

The variance-inflation formula and the estimator of $\rho$ are documented in ?spiGSEA; the gene-level benchmark this layer sits on is in Simulating realistic spatial data and benchmarking spiDE.

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3    spiDE_0.99.17    BiocStyle_2.40.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] SummarizedExperiment_1.42.0 gtable_0.3.6               
#>  [3] rjson_0.2.23                xfun_0.60                  
#>  [5] bslib_0.12.0                Biobase_2.72.0             
#>  [7] lattice_0.22-9              vctrs_0.7.3                
#>  [9] tools_4.6.1                 generics_0.1.4             
#> [11] stats4_4.6.1                parallel_4.6.1             
#> [13] tibble_3.3.1                pkgconfig_2.0.3            
#> [15] Matrix_1.7-5                RColorBrewer_1.1-3         
#> [17] S7_0.2.2                    desc_1.4.3                 
#> [19] S4Vectors_0.50.2            lifecycle_1.0.5            
#> [21] compiler_4.6.1              farver_2.1.2               
#> [23] textshaping_1.0.5           statmod_1.5.2              
#> [25] Seqinfo_1.2.0               codetools_0.2-20           
#> [27] htmltools_0.5.9             sass_0.4.10                
#> [29] yaml_2.3.12                 pkgdown_2.2.1              
#> [31] pillar_1.11.1               jquerylib_0.1.4            
#> [33] BiocParallel_1.46.0         SingleCellExperiment_1.34.0
#> [35] DelayedArray_0.38.2         cachem_1.1.0               
#> [37] limma_3.68.5                magick_2.9.1               
#> [39] abind_1.4-8                 tidyselect_1.2.1           
#> [41] digest_0.6.39               dplyr_1.2.1                
#> [43] bookdown_0.48               labeling_0.4.3             
#> [45] fastmap_1.2.0               grid_4.6.1                 
#> [47] cli_3.6.6                   SparseArray_1.12.2         
#> [49] magrittr_2.0.5              S4Arrays_1.12.0            
#> [51] withr_3.0.3                 scales_1.4.0               
#> [53] rmarkdown_2.32              XVector_0.52.0             
#> [55] matrixStats_1.5.0           otel_0.2.0                 
#> [57] ragg_1.5.2                  SpatialExperiment_1.22.0   
#> [59] evaluate_1.0.5              knitr_1.52                 
#> [61] GenomicRanges_1.64.0        IRanges_2.46.0             
#> [63] rlang_1.3.0                 Rcpp_1.1.2                 
#> [65] glue_1.8.1                  BiocManager_1.30.27        
#> [67] BiocGenerics_0.58.1         jsonlite_2.0.0             
#> [69] R6_2.6.1                    MatrixGenerics_1.24.0      
#> [71] systemfonts_1.3.2           fs_2.1.0