Gene-set enrichment on a fitted spiDE model

Tests whether the genes of a set shift coherently in their response-niche (or cell-type response) statistics. Gene-set statistics average over tens to hundreds of genes, which recovers the power the per-gene niche tests lose on sparse spatial data, with the inter-gene correlation carried explicitly so that co-regulation is not mistaken for evidence.

Usage

spiGSEA(object, spe = NULL, genesets, ...)

# S4 method for class 'SpiDEResults'
spiGSEA(
  object,
  spe = NULL,
  genesets,
  type = c("niche", "celltype"),
  test = c("competitive", "self-contained"),
  fdr = 0.05,
  min.size = 5L,
  max.size = 500L,
  rho = NULL,
  rho.genes = 2000L,
  block.size = NULL,
  backend = c("auto", "cpu", "gpu"),
  gpu.mem.budget = NULL,
  BPPARAM = BiocParallel::SerialParam(),
  verbose = TRUE
)

Arguments

object: a SpiDEResults from testSpiDE() or spiDE().
spe: the SpatialExperiment::SpatialExperiment the model was fitted on. Only needed when the inter-gene correlation has to be estimated – testSpiDE() normally stores it on each fit, in which case this can be omitted and no pass over the counts happens at all.
genesets: a named list of gene identifiers (character, matched against rownames) or integer row indices.
...: ignored; present for compatibility with the generic.
type: on a condition-free fit (condition = NULL) "niche" tests the two-way celltype:niche statistics, which are the tested effects in that mode, and "celltype" has nothing to test and errors. Otherwise: "niche" (default) tests the three-way celltype:condition:niche statistics; "celltype" tests the CellType:condition statistics, which are empty unless the design carries that block.
test: "competitive" (default) compares the set's mean statistic with the genes outside the set, as limma::camera does; "self-contained" compares it with zero. See the section below: the benchmark shows the self-contained form does not control error on correlated data, which is why it is not the default.
fdr: the FDR threshold applied at every level of the cascade.
min.size, max.size: sets outside this size range (after intersecting with the fitted genes) are dropped. The lower bound keeps the mean from being dominated by one gene; the upper bound drops sets so broad that the competitive null they are tested against is largely themselves.
rho: the average inter-gene correlation. NULL (default) uses the value testSpiDE() stored on each fit, and falls back to estimating it from Pearson residuals – separately per bandwidth, since each is a different design and leaves different residuals – only if the fit does not carry one. Supply a scalar to use one value throughout, or one value per bandwidth, to override. Pass 0 only if you are certain the genes are independent – they are not, and the test is severely anti-conservative without this term.
rho.genes: number of genes to subsample when estimating rho (NULL uses all). A random subset of genes gives a random subset of gene PAIRS, so the estimate is unbiased; the default trades a little precision on a nuisance scalar for a large saving on the counts pass.
block.size: genes per block when streaming the counts. NULL auto-sizes against the device memory budget on the GPU backend.
backend: the compute backend for the inter-gene correlation pass ("auto", "cpu" or "gpu"), which is the only step that touches the counts and so the only one worth accelerating – the gene-set arithmetic downstream is a sets x contrasts matrix and negligible by comparison. Ignored when rho is supplied, since no pass is then needed.
gpu.mem.budget: device memory budget in bytes for the GPU block sizer.
BPPARAM: a BiocParallelParam. Gene blocks contribute additive partial sums, so they parallelise exactly; a multi-worker param is downgraded to serial on the GPU backend to stop workers contending for one device.
verbose: report progress.

Value

a data.frame of significant sets, keyed by (geneset, ct_index) and additionally ct_niche when type = "niche", with the set size, the z-statistic at its most informative bandwidth, the direction, and the q-value from each level of the cascade. Empty if nothing survives.

Details

Evidence is combined across bandwidths with the same log-likelihood-weighted Cauchy combination used for the gene-level results, then gated by a nested Benjamini-Hochberg cascade: set level, then index cell type within surviving sets, then (for type = "niche") niche cell type within those.

Which test, and what it licenses you to say

test = "competitive" (the default) asks whether the set's mean differs from the mean of the genes outside it, as limma::camera does – null: the set responds no more than the assayed background.

test = "self-contained" asks whether the set's mean differs from zero – null: no gene in the set responds. It reproduces the flat-script pipeline this replaces, and is offered for that comparison, but it does not control error on correlated data and should not be used for inference.

The simulation benchmark (research/, scenario gsea) measured this on ground truth: with nothing planted, at realistic inter-gene correlation and a nominal FDR of 0.05, the self-contained test called 20.6 of 208 sets per replicate – every one of them false, a realised FDP of 1.00 – while the competitive test called 0.05 sets. The cause is not the correlation term: it is that the self-contained form assumes the gene-level statistics being averaged have unit spread, which holds under the null but not under signal (observed spread 1.8 at the largest effect tested). The competitive form divides by the observed spread and so is immune.

Two consequences worth knowing. Under a global shift affecting most genes, the self-contained test called 96.3% of random sets and the competitive test 0.7%. And because spiDE's per-gene statistics track expression, the self-contained test's false calls concentrate in abundant sets: its called sets sat at the 86th expression percentile against the 51st for random sets. Where that distinction is load bearing, also compare each set against expression-matched random sets of the same size rather than reading the q-value alone.

References

Wu D, Smyth GK (2012). "Camera: a competitive gene set test accounting for inter-gene correlation." Nucleic Acids Research 40(17):e133. (Source of the variance inflation factor, and of the two-sample form used by test = "competitive".)

Wu D, Lim E, Vaillant F, Asselin-Labat ML, Visvader JE, Smyth GK (2010). "ROAST: rotation gene set tests for complex microarray experiments." Bioinformatics 26(17):2176-2182. (Self-contained gene-set null.)

Examples

data(toySpiDE)
spe <- buildNiches(toySpiDE, sigma = 20)
res <- spiDE(spe, condition = "condition", sigma = 20, random = "none", verbose = FALSE)
gs <- list(setA = rownames(spe)[1:5], setB = rownames(spe)[6:10])
spiGSEA(res, spe, gs, min.size = 3, fdr = 1)
#> spiGSEA: 2 of 2 sets within [3, 500] genes
#> spiGSEA: using the rho stored on the fit by testSpiDE()
#> spiGSEA: sigma 20  rho = 0.0046
#> spiGSEA: 12 significant rows over 2 sets
#>    geneset collection size ct_index ct_niche bandwidth.max          z Direction
#> 1     setA        all    5        C        A            20  1.6687975        Up
#> 2     setB        all    5        A        C            20  1.6551465        Up
#> 3     setA        all    5        A        B            20  1.3134296        Up
#> 4     setA        all    5        A        C            20  0.8807061        Up
#> 5     setA        all    5        C        B            20 -0.6861079      Down
#> 6     setA        all    5        B        C            20 -1.0701246      Down
#> 7     setB        all    5        B        A            20 -0.8886714      Down
#> 8     setB        all    5        B        C            20 -0.4207446      Down
#> 9     setB        all    5        A        B            20  0.3694734        Up
#> 10    setA        all    5        B        A            20 -0.1601932      Down
#> 11    setB        all    5        C        B            20 -0.5561613      Down
#> 12    setB        all    5        C        A            20 -0.1169303      Down
#>    fdr.geneset fdr.index fdr.niche
#> 1    0.5729268 0.3894411 0.1903150
#> 2    0.5729268 0.6705865 0.1957896
#> 3    0.5729268 0.3894411 0.3780764
#> 4    0.5729268 0.3894411 0.3784769
#> 5    0.5729268 0.3894411 0.4926451
#> 6    0.5729268 0.7111448 0.5691265
#> 7    0.5729268 0.7954601 0.6739416
#> 8    0.5729268 0.7954601 0.6739416
#> 9    0.5729268 0.6705865 0.7117749
#> 10   0.5729268 0.7111448 0.8727289
#> 11   0.5729268 0.8375260 0.9069153
#> 12   0.5729268 0.8375260 0.9069153