Simulating realistic spatial data and benchmarking spiDE
Dharmesh D. Bhuva
Introduction
spiDE tests for neighbourhood-dependent
differential expression. Within an index cell type, it
asks whether a gene’s response to an experimental condition
depends on the local density of surrounding cell types — the
niche. This is a three-way interaction:
cell type × condition × niche density. Detecting it
reliably is hard, because spatial transcriptomics data are sparse, the
biological signal is subtle, and cells from the same patient are not
independent observations.
To trust a method that looks for such a subtle effect, we need two things: simulated data whose statistical structure genuinely resembles real spatial transcriptomics, and a benchmark that measures the method against known ground truth. This report provides both.
It is organised in two halves:
- The simulator. A detailed description of the generative model, followed by a guided tour of one simulated dataset — what the counts look like, how genes are categorised, how the differential-expression effect is planted, and how cells are arranged in space.
- The benchmark. A large study that fits
spiDEto thousands of simulated datasets and measures its calibration (false-positive control), statistical power, false-discovery-rate control, robustness to spatial arrangement, and computational cost. Every summary figure is accompanied by an explanation of what it plots and how to read it.
The benchmark figures are read from the canonical summary tables in
the research repository (reports/benchmarks/tables/). They
are produced by a self-contained study designed to run on a compute
cluster; the code and instructions for regenerating the tables live in
the package’s research/ directory.
Part 1 — The simulation framework
The generative model
The simulator builds a SpatialExperiment one patient
sample at a time and is inspired by the negative-binomial model of the
splatter
package (Zappia, Phipson & Oshlack, Genome Biology 2017).
Each ingredient targets a specific property of real data.
Gene-level abundance and overdispersion. Every gene \(g\) has a baseline mean expression \(\mu_g\). Counts are drawn from a negative binomial whose dispersion is not constant but follows a mean–variance trend: the biological coefficient of variation (BCV) falls with abundance, \[\text{BCV}_g = b_0 + b_1 / \sqrt{\mu_g}, \qquad \text{size}_g = 1/\text{BCV}_g^2,\] so lowly expressed genes are relatively far more variable than highly expressed ones — exactly as observed empirically. Defaults are \(b_0 = 0.3\), \(b_1 = 1\).
Three gene categories. Real panels mix genes with very different roles, and the simulator mirrors this with three categories (summarised below). This matters because the genes that carry a niche effect in real tissue are expressed genes — cell-type markers and metabolic/housekeeping programs — not the silent background.
| category | role | expression | typical detection |
|---|---|---|---|
| cell-type marker | defines a cell type; expressed in that type, near-silent elsewhere | low mean in its own type, very low background elsewhere | ~10–20 % of cells of its type |
| housekeeping | broadly active programs (incl. metabolism) | moderate mean in all cells | ~40–55 % of all cells |
| background | silent / lowly expressed | very low everywhere | ~1–4 % of cells |
Cell types and space. Each cell is assigned a type
and a position in a square field. One type is the index
(A); one is the niche (B) whose local
density drives the effect; the rest are bystanders (C, and
D in the multi-niche layout). Niche cells are placed so
their density varies in space according to the chosen layout
(Part 1, “The four spatial layouts”).
The planted effect (the ground truth). A small
fraction of genes are made niche-responsive. For these genes,
in index cells of Responder samples,
the log-mean is raised in proportion to the local niche density: \[\log \mu_{g,i} \;\mathrel{+}=\; \beta \cdot
\mathbb{1}[\text{index cell}]\cdot
\mathbb{1}[\text{Responder}] \cdot \rho_i,\] where \(\rho_i \in [0,1]\) is the niche potential
at cell \(i\) and \(\beta\) is the effect size. The responsive
genes are drawn from the expressed categories — a mix of
index-cell-type markers and housekeeping genes — so the effect lands
where it plausibly would in tissue. This term, and nothing else, is what
spiDE must recover.
Patient structure. A per-(gene, sample) random intercept adds between-patient variability, so that cells from one patient are correlated — the pseudo-replication that makes naive cell-level inference anti-conservative. Each patient is labelled Responder or Non-responder; the condition is shared by all of a patient’s cells.
Library size and filtering. Each cell has a log-normal size factor scaling its total counts, and genes below a low total-count floor are dropped (as real pipelines do), except planted true positives which are always kept.
The chunk below implements this model compactly, drawing on the
package’s own count engine (spiDE:::.simGeneParams,
spiDE:::.simCounts), and generates one reference dataset of
24 patient samples that the rest of Part 1
dissects.
#> reference dataset: 4000 genes x 3840 cells across 24 samplesWhat the counts look like
Individual gene–cell measurements are dominated by zeros and ones, yet each cell still accumulates a realistic total depth because the (few) housekeeping genes carry most of the reads.
Overdispersion decreases with abundance
Estimating each gene’s mean and BCV from the counts recovers the characteristic downward mean–variance trend the model was built from: rare genes are the most overdispersed.
The three gene categories
The categories differ in both how much they are expressed and where. The left panel shows their mean-expression distributions; the right panel shows the fraction of cells in which each category is detected — the property that makes markers usable as cell-type labels and keeps the background quiet.
Anatomy of one dataset: gene categories and the planted effect
The two heatmaps below dissect the reference dataset. The first shows what makes a gene a marker; the second shows the planted differential-expression effect across all 24 patients.
Cell-type specificity. For a representative sample of genes from each category, we compute the fraction of each cell type’s cells in which the gene is detected (a non-zero count), pooled over all patients. Cell-type markers are detected mainly in a single cell type; housekeeping genes are detected across all types; background genes are rarely detected anywhere.
The planted effect. Restricting to the index
cell type A, we compute each gene’s mean
expression in every patient. Rows are genes split by category
and by whether they carry the planted effect (DE vs non-DE) —
including a row of DE background genes, since the full
study spikes the effect into 300 background genes alongside 300 markers
and 300 housekeeping genes; columns are the 24 patients split into
Responders and Non-responders. Only the DE genes — and only in Responder
patients — are systematically brighter: this is the condition-dependent,
niche-driven signal spiDE is designed to find. Non-DE genes
and the entire Non-responder column are flat by construction. Note the
DE background row is much dimmer than the DE marker/housekeeping rows:
the same planted log-fold change is far harder to see against a
near-silent baseline, which foreshadows the category breakdown in Part
2.
The four spatial layouts
The niche effect is only meaningful relative to where the niche cell type sits. The study uses four spatial layouts, each modelling a different kind of tissue organisation and each producing a differently shaped ground-truth signal.
- gradient — the niche cell type
Bincreases smoothly in density from one side of the tissue to the other. This models a monotonic tissue axis (for example distance from a boundary, a vessel, or the tissue edge). The planted signal in index cells is a smooth ramp. - clustered —
Bcells concentrate in a handful of Gaussian hotspots, leaving the rest of the tissue sparse. This models focal structures such as immune aggregates or tumour nests, where the niche is patchy. The planted signal appears as bright islands. - random — cell types are placed with no spatial structure at all; the niche “density” is essentially constant everywhere. This is a negative control: there is no real spatial niche, so a well-behaved method should find no niche effect. It is the reference against which false spatial discoveries are measured.
- multiniche — the
Bniche (patchy) coexists with a second, independent structured cell typeDoccupying different hotspots. This is a specificity test: the method must attribute the effect toBand not be misled by the competing structure ofD.
The chunk below draws one sample of each layout and shows three
views, stacked as rows: where each cell type sits, the resulting
true planted signal (deterministic, before any sampling noise)
in the index cell type, and — the step that actually reaches
spiDE — the noisy simulated log-counts of one example
marker gene once that signal has been passed through the negative
binomial count model.
Notice that under the random layout the true signal
is spatially flat: there is nothing for buildNiches() to
lock onto, which is exactly why it serves as the negative control. The
bottom row makes the detection problem concrete: even where a real ramp
or hotspot exists (gradient, clustered, multiniche), the per-cell counts
are dominated by Poisson-like sampling noise on top of the signal — a
single cell’s log-count is a very weak clue on its own, which is why
spiDE pools evidence across many cells and many patients
rather than testing cell-by-cell.
Part 2 — Benchmarking spiDE
The remainder of this report summarises a study that fits
spiDE to many datasets drawn from the model above, sweeping
the number of patient samples, the number of cells per sample, the
effect size and the spatial layout, and repeating each design point over
many random replicates. Each dataset is fit once and every metric is
derived from that fit.
The benchmark experiments
Every result in Part 2 comes from a single study, defined in
research/R/scenarios.R and run on a compute cluster. It is
organised as four experiment families. Each family expands into (design
point × replicate) jobs; each simulated dataset is fit
once, and every metric is re-derived from that one fit
(re-gating discoveries at each FDR threshold rather than re-fitting).
The full study uses 40 replicates per design point, and
every summary metric below is reported as a mean over those replicates
(medians, standard deviations and quartiles are also computed, for error
bars).
The niche effect is spiked into a fixed number of genes from
each gene class — 300 cell-type markers, 300 housekeeping genes
and 300 background genes — so that detectability can be compared across
classes. Most other settings are held fixed while one axis at a time is
swept, so each figure isolates a single variable. The table below is
generated from the shipped summary tables wherever a quantity can be
read off them; a handful of rows describe simulator inputs that never
appear as a column in the output (gene-category counts, the niche
bandwidth, the between-patient SD) and are typed from
research/config.R directly.
| quantity | fixed value | swept range (in which family) |
|---|---|---|
| total genes | ~10,000 | — |
| — cell-type markers (per type) | 500 | — |
| — housekeeping genes | 1,500 | — |
| — background genes | remainder (~7,000) | — |
| true niche-responsive genes | 900: 300 markers + 300 housekeeping + 300 background | — |
patient samples, S |
24 | 4, 10, 16, 24, 30 (null, power) |
| cells per sample | 300 | ** 25, 100, 250, 500, 1,000** (power) |
effect size, β |
2 | 1, 1.5, 2, 3 (power) |
| spatial layout | gradient | clustered, gradient, multiniche, random (power) |
niche bandwidth, σ |
50 | — |
| inference mode | intercept | fixed, intercept, slope, twostage (null, power) |
| p-value combiner | Cauchy | brown, cauchy (combiner) |
PQL cell fraction, re.prop |
0.2 | 0.05, 0.1, 0.2, 0.25, 1 (pql) |
| between-patient SD | 0.5 | 0.6 (pql) |
| nominal FDR thresholds | — | 0.01, 0.05, 0.1, 0.2 |
| replicates per design point | 40 | — |
The four families and what each varies:
| family | question it answers | swept axis | held fixed | modes |
|---|---|---|---|---|
| null | type-I error / calibration | S × {gradient, clustered}, β = 0 |
300 cells | fixed / intercept / slope / twostage |
| power | power, FDR, recall (overall and per gene class), spatial robustness and fit/inference timing | four sweeps: S; cells; β;
layout |
gradient layout (except the layout sweep); the other axes at their fixed value | fixed / intercept / slope / twostage |
| combiner | Cauchy vs Brown | null (β = 0) and alternative arms |
S = 24, 300 cells, β = 2 |
intercept |
| pql | subsampling speed/accuracy | re.prop |
S = 24, 1,000 cells |
intercept |
The power family now also carries the layout sweep, so
its results supply the spatial-robustness metrics, and
its per-fit wall-clock times supply the scalability
metrics — dedicated spatial and timing families are no longer required.
The combiner and pql families are reported in
the two companion vignettes.
Every summary table carries, for each metric m, the mean
(bare column m) plus m_sd and m_n
across replicates. The standard error of the mean is
m_sd / sqrt(m_n); the se_of() helper defined
above returns it, and errbars() adds a matching error bar
to any figure whose y aesthetic is the mean of metric m.
Every figure below therefore shows replicate-to-replicate uncertainty,
not just a point estimate.
Throughout, four arms are compared — three inference modes of the one-stage GLM, plus a separate estimator:
- fixed (
random = "none") — every cell is treated as an independent observation; - intercept (
random = "intercept") — a per-patient random intercept absorbs between-patient shifts; - slope (
random = "slope") — additionally allows the niche effect itself to vary between patients; - twostage (
twoStageSpiDE()) — a different estimator entirely: each patient’s niche slope is estimated first, and the condition contrast is then taken over the patient slopes with a moderated t-test. Its rows are computed on the same simulated datasets (the two-stage harness seeds from this study’s generator), with the arm shown beingstage1 = "ols"; the dedicated vignette Two-stage estimation: benchmarked against the published simulation study documents that choice, the scoring differences (each arm is gated by the FDR machinery it ships), and the results in full. Here it is drawn alongside the GLM modes for context.
The mixed-effects modes exist because a patient contributes many cells but only one condition label, so cells are pseudo-replicates. The benchmark quantifies what that pseudo-replication does to calibration, power and false discoveries; the two-stage estimator instead removes the pseudo-replication by construction, at a cost in power measured below.
Calibration
Each default of the estimator is measured here on the null family, in
the order the model vignette introduces them: the random mode and the
reference degrees of freedom, the CellType:condition design
term, and the 0.99.17 switches.
The random modes and the reference degrees of freedom
What is plotted. Datasets are simulated with no planted effect. A well-calibrated test should then reject the null hypothesis at exactly the chosen significance level \(\alpha = 0.05\). The figure plots the realised rejection rate against the number of patients, for each mode and two layouts. The dashed line is the target, \(0.05\).
How to read it. A curve on the dashed line is perfectly calibrated. Above the line means anti-conservative — the test declares significance too readily, which would inflate false discoveries. Below the line means conservative — the test is cautious and will sacrifice some power.
How it is computed. From the null family
(no planted effect, β = 0): for each of the 40 replicates
we record the fraction of all ~10,000 genes whose two-sided
ResponseNiche Wald p-value — the single index-A × condition
× B-niche coefficient — falls below 0.05, then average over
replicates. Two layouts (gradient, clustered) and five sample sizes are
shown.
Interpretation. The two reference-df choices behave very differently, and the gap is widest exactly where studies are smallest.
Under the legacy between-patient reference (solid lines), the mixed modes are severely over-conservative with only a handful of patients: at \(S=4\) the \(S-2=2\) reference df pushes type-I error to \(\approx 0.001\), a 50-fold undershoot of the nominal \(0.05\), and it climbs back toward nominal only as patients accrue. That is not free caution — the power figures below show it costs essentially all recall at small \(S\).
The Satterthwaite reference (dashed) targets
precisely this, giving each coefficient the df its own information
content earns rather than one flat \(S-2\). It holds close to nominal across the
entire sampled range. The honest caveat is that it does
overshoot mildly at larger \(S\) — worst measured \(\approx 0.065\), a \(\approx 30\%\) relative inflation — so it
is slightly liberal where between is drastically
conservative. That residual points at the covariance, and a
Kenward–Roger correction is the indicated next step (see The spiDE
model). Trading a bounded \(30\%\)
overshoot for a 50-fold undershoot, plus the power it unlocks, is why
Satterthwaite is now the package default.
The fixed and twostage arms have no
df.method to vary (neither uses a between-patient reference
df — two-stage inherits the moderated-t df of its patient-level
limma contrast), so both appear solid throughout. They sit at opposite
ends of the calibration picture: fixed is drastically
conservative at every design point, while twostage is
the best-calibrated arm in this simulation study — essentially on the
nominal line at \(S = 4\) and drifting
mildly conservative as patients accrue, the mirror image of the
mixed modes’ mild anti-conservative drift. The dedicated two-stage
benchmark vignette unpacks that trade (calibration bought at a cost in
power); note the ordering reverses on the real cohort, where
the niche-shuffle report shows the mixed GLM better calibrated across
all index types.
A quantile–quantile (QQ) plot makes the same calibration check across the whole distribution of p-values rather than at one threshold. Because each point is an order statistic pooled over many genes and replicates rather than an independent replicate mean, there is no per-point standard error to draw here — the type-I error and FDP figures above and below are where replicate-to-replicate uncertainty is quantified.
What is plotted. For the largest sample size, the sorted null p-values are plotted (on a \(-\log_{10}\) scale) against their expected values under perfect calibration, split by mode (colour) and layout (facet — both layouts available in the null family). How to read it. Points on the diagonal are calibrated; points above the diagonal are anti-conservative (too many small p-values — excess false positives); points below are conservative.
The Satterthwaite reference does not just fix small-\(S\) calibration — it also delivers real power the between reference cannot, most dramatically where between’s severe small-\(S\) conservatism leaves it with almost nothing to detect. The figures below (raw power, then recall after FDR control) show this directly.
The design term: CellType:condition
The design carries one CellType:condition column per
cell type — the flat, cell-type-wide response — so that the tested
three-way slope estimates the gradient of the response along
the niche and not its level. The alternative, which has no such column
and lets a global condition effect stand for every type, is
run as a second design arm (niche-only) over the families
that can discriminate between them. Both are calibrated:
The arm without the term reports more niche calls, and the difference is large:
That is not power. The two arms estimate different quantities.
Without the term, the flat part of the planted effect — \(\beta\) times the mean niche potential, a
shift in every index cell of every Responder — has nowhere to go but the
uncentred niche slope, and one power dataset fit under builds that
differ only by the term shows the niche-only slope exceeding
the other by \(0.25\times\) the
flat-response coefficient, gene by gene (\(r =
0.88\) on planted genes, \(0.92\) on null genes). That inflates the
estimate by about 10%. The \(t\) is
\(1.5\)–\(1.6\times\) larger because the standard
error is \(0.67\)–\(0.72\times\): with no intercept column in
that subspace, the slope is a regression through the origin, whose
precision is bought by assuming there is no flat shift to estimate.
Two-sided Cauchy combination gives more discoveries than one-sided Brown
in every build, so sidedness is not involved; the recall saturates at
about 0.5 from \(\beta = 2\); and on
the real cohort the design with the term recovered 41% more
triplets. The ctresp family makes the same point from the
other side by planting a flat shift with no niche
effect, so that every niche call on those genes is false by
construction:
#> *(No `ctresp` results in the shipped summaries yet.)*The design without the term produces spurious niche calls at every effect size, with a rate that climbs steeply as the FDR gate loosens — roughly \(9.5\times\) the rate with the term at \(\beta_{ct} = 1\) — because the leaked shift sits just below the gate. Null calibration is the same either way (median type-I 0.0450 against 0.0449 across the design points), so the term is a question of the estimand, not of validity, and the estimand the package ships is the within-cell-type gradient.
The 0.99.17 defaults: nested intercept and per-gene convergence
fitSpiDE(re.celltype = TRUE, converge = TRUE) — a
ridge-penalised (sample × cell type) intercept block and per-gene
convergence — became the defaults on real-cohort evidence: without the
nested block the tested slopes carried a between-patient composition
association that the niche-shuffle null preserved (The real
cohort report). They are run here as three switch arms on one
frozen 0.99.17 library, nested-converged (both),
nested-only and converged-only, in intercept
mode. One caveat governs the reading: the simulator places niche cells
by the same potential in every sample, so a sample’s mean niche density
varies only by sampling noise and is independent of its expression
intercept. It plants no composition effect, so on this
study the nested block can only cost, never pay.
| S | design | tpr | fdp |
|---|---|---|---|
| 4 | celltype-response | 0.040 | 0.346 |
| 4 | nested-converged | 0.007 | 0.590 |
| 10 | celltype-response | 0.068 | 0.072 |
| 10 | nested-converged | 0.072 | 0.059 |
| 16 | celltype-response | 0.193 | 0.041 |
| 16 | nested-converged | 0.205 | 0.033 |
| 24 | celltype-response | 0.315 | 0.034 |
| 24 | nested-converged | 0.350 | 0.028 |
| 30 | celltype-response | 0.361 | 0.033 |
| 30 | nested-converged | 0.402 | 0.031 |
Averaged over \(S \ge 10\) and both
layouts the shipped design term alone is at 0.058,
converged-only at 0.066, nested-only at 0.068
and the pair at 0.071, with a per-cell standard error of about 0.0016.
The two switches add roughly and act at different sample sizes. At \(S = 4\) the nested block is harmless
(0.037–0.043) and convergence carries the whole excess (0.086–0.095):
per-gene profile-ML dispersion gives up edgeR’s cross-gene moderation,
which matters most when data are few. From \(S
= 16\) convergence fades onto the design-term arm while the
nested block holds a constant \(+0.010\), which points at its reference df
or the estimation of its variance component. The pair gains recall
(0.402 against 0.361 at \(S = 30\))
with slightly better FDP from \(S =
10\) and a worse one at \(S =
4\) (0.59 against 0.35). None of this is a reason to revert: the
defaults were chosen for a confound that real cohorts carry and this
simulator does not contain. It is the measurement the two candidate
refinements — a moderated dispersion at the converged mean, and the
nested block’s reference df — are judged against; the first is measured
next, and it is not the dispersion that matters.
The dispersion rule and the standard-error scale (0.99.18)
Four arms on the same frozen build cross the convergence stage’s
dispersion rule (per-gene profile-ML, or fitNB’s cross-gene
moderated value kept at the converged mean) with the scale of the
standard errors (the working Pearson dispersion, or the edgeR-v4
quasi-likelihood dispersion moderated across genes), all with the nested
block and the convergence step, in intercept mode with Satterthwaite
df.
| S | design | tpr | fdp |
|---|---|---|---|
| 4 | celltype-response | 0.040 | 0.346 |
| 4 | nested-converged | 0.007 | 0.590 |
| 4 | psi-moderated-ql | 0.001 | 0.000 |
| 10 | celltype-response | 0.068 | 0.072 |
| 10 | nested-converged | 0.072 | 0.059 |
| 10 | psi-moderated-ql | 0.044 | 0.005 |
| 16 | celltype-response | 0.193 | 0.041 |
| 16 | nested-converged | 0.205 | 0.033 |
| 16 | psi-moderated-ql | 0.196 | 0.013 |
| 24 | celltype-response | 0.315 | 0.034 |
| 24 | nested-converged | 0.350 | 0.028 |
| 30 | celltype-response | 0.361 | 0.033 |
| 30 | nested-converged | 0.402 | 0.031 |
| 30 | psi-moderated-ql | 0.360 | 0.014 |
The dispersion rule is irrelevant: the two Pearson arms agree to the third decimal at every \(S\) (0.091 at \(S = 4\), 0.068 at \(S = 30\)), and so do the two QL arms, so the small-\(S\) convergence cost above was never the loss of moderation. The scale is the whole difference. The QL scale is the first configuration measured on this study that holds the nominal level at every \(S\) including \(S = 4\) (0.050–0.055 throughout). Its recall at FDR 0.05 is that of the design term alone (0.360 against 0.361 at \(S = 30\); the 0.99.17 pair’s 0.402 is bought with its inflated null — raw power at 0.05 is 0.666, 0.660 and 0.694 respectively) and its realised FDP sits far below nominal (0.014 at \(S = 30\), 0.000 at \(S = 4\) against 0.59), which says the hierarchical cascade is conservative on well-calibrated p-values — a separate lever, unmeasured. Per task, the moderated dispersion costs 0.6× the profile one from \(S = 16\) (the profile step is the expensive part of the convergence stage) and the QL pre-pass adds about 5% there and more only on tiny designs; on the real cohort both variants call the same triplets as the 0.99.17 defaults (The real cohort report). From 0.99.18 the moderated dispersion and the QL scale are the package defaults, with the machinery in SpaNorm.
Statistical power: what buys detection?
What is plotted. With a real effect planted, raw power is the fraction of truly niche-responsive genes the test flags before multiple-testing correction — a measure of the per-gene test’s sensitivity, pooled across all three gene categories (a category-resolved version follows in “Which genes are detected?” below). The three panels sweep, respectively, the number of patient samples, the number of cells per sample, and the effect size \(\beta\). Higher is better; there is no reference line.
How it is computed. From the power family,
in which ~9 % of genes (900 of ~10,000: 300 markers + 300 housekeeping +
300 background) carry the planted effect. Raw power is the fraction of
those 900 true genes whose uncorrected ResponseNiche p-value is below
0.05, averaged over 40 replicates. Each panel sweeps one axis while the
other two are held at their fixed values (S = 24, 300
cells, β = 2).
Interpretation. All three levers increase power, and
the mixed-effects modes are uniformly more powerful than the
fixed-effects mode — because modelling the between-patient variance
correctly leaves a cleaner estimate of the effect. The
twostage arm sits between them: above
fixed from \(S = 10\)
onward but below the mixed modes throughout, part of which is
calibration rather than information (the null figure above shows the
mixed modes buying detections with an inflated type-I error where
two-stage sits at or under nominal). Three points are worth drawing out
for the GLM modes:
The reference df is itself a power lever. Compare the solid (
between) and dashed (Satterthwaite) lines for the same mode: the gap is dramatic at small \(S\), where the flat \(S-2\) reference is so conservative that the mixed modes detect almost nothing, and narrows as patients accrue andbetweenbecomes less punitive. This is the power counterpart of the calibration figure above, and the reason Satterthwaite is the default — the conservatismbetweenbuys at small \(S\) is paid for almost entirely in lost detections. Note the quoted figures in the two bullets below are for thebetweenarm.More patients help the most. Power climbs steeply with the number of samples (from ~0.02 at 4 patients to ~0.64 at 30 for the mixed modes), because the condition contrast is fundamentally a between-patient comparison.
More cells also help here — unlike a simple between-patient mean contrast, the niche effect is a within-patient spatial gradient, so additional cells genuinely sharpen each patient’s estimate of that gradient (power rises from ~0.37 at 100 cells to ~0.78 at 1,000). Cells are informative for the niche interaction even though they cannot manufacture independent evidence about the condition label itself.
From raw power to controlled discoveries
Raw power measures a single gene’s test sensitivity. What a user actually receives is the set of genes that survive genome-wide multiplicity control, and that set is much smaller. This section makes the distinction explicit — it is the reason raw power can be high while recall stays low — and shows across conditions where discoveries are made and where the false-discovery rate is controlled.
How discoveries and their errors are computed. From
the power family: at each design point all ~10,000 genes are
fit, and a gene is discovered when it clears
spiDE’s three-level hierarchical Benjamini–Hochberg gate —
gene level, then index cell type, then niche cell type — at the nominal
threshold. The realised false-discovery proportion
(FDP) is the fraction of discoveries that are not among the 900
true genes; recall is the fraction of the 900 true
genes discovered. Both are averaged over 40 replicates.
Why raw power is high but recall is low
What is plotted. For the intercept mode across the samples sweep, raw (uncorrected) power and recall at a 5 % FDR on the same axes. Higher is better; the vertical distance between the curves is the cost of multiplicity control.
Interpretation. Raw power rises smoothly with sample size, reaching ~0.33 by 10 patients and ~0.64 by 30. Recall behaves completely differently: it stays near zero until the study is large enough — here about 24 patients — and then climbs steeply. This threshold is characteristic of genome-wide multiplicity control. A gene is only discovered if it clears the hierarchical Benjamini–Hochberg gate, whose gene-level step must separate the 900 true genes from ~9,100 nulls. A raw p-value just under 0.05 is nowhere near small enough to survive that correction; only once the effect is strong enough to push the true genes’ p-values deep into the tail do they start to clear the bar — at which point many clear it together. High raw power is therefore necessary but not sufficient for recall: raw power says the per-gene test can see the effect, while recall says the effect is strong enough to be reported genome-wide at a controlled error rate.
The same threshold governs cells and effect size
What is plotted. Recall at a fixed 5 % FDR as each experimental axis is swept in turn, for all three modes, pooled across all three gene categories (see “Which genes are detected?” below for the category-resolved version). Higher is better; there is no reference line.
Interpretation. More cells per sample and larger effects cross the same threshold: recall is ~0 below it (≤100 cells, or \(\beta \le\) 1.5) and rises once the signal clears genome-wide correction (to ~0.45 at 1,000 cells, ~0.38 at \(\beta =\) 3). So the practical message is not “cells don’t help” — they clearly do here — but that some minimum combination of patients, cells and effect size is needed before any discoveries are made at all.
Where discoveries are made, is the FDR controlled?
What is plotted. The realised FDP at a nominal 5 % FDR for each mode, shown only where a mode makes an appreciable number of discoveries (recall > 5 %) — below the detection threshold the FDP is the ratio of a handful of calls and is not meaningful. As with power and recall above, the three panels split by sample size, cells per sample and effect size. The dashed line is the nominal 0.05. How to read it. On or below the line means control holds; above means it fails.
Interpretation. Wherever a mode actually makes discoveries, the mixed-effects modes hold FDP close to the nominal 0.05 on all three axes (~0.004–0.028) — with one point worth flagging honestly rather than glossing over: at the largest effect size studied (\(\beta =\) 2), FDP creeps mildly above nominal for the mixed modes too (~0.03), though still well below the fixed mode’s anti-conservatism at the same point (~0.04). The fixed-effects mode is anti-conservative whenever discoveries are scarce or the signal is strong — its worst points on each axis are ~0.24 at 4 patients, ~0.52 at 25 cells, and ~0.33 at \(\beta =\) 1 — and while it settles closer to nominal as more patients are added (~0.03 by 30 patients), it does not reliably improve with more cells or a bigger effect: those levers buy the fixed mode more discoveries, not better-calibrated ones. The panels start later on the x-axis where a mode makes essentially no discoveries (recall ≤ 5 %), so there is nothing meaningful to plot until it starts discovering. The twostage arm, where it discovers at all, holds FDP at ~0.047–0.233 — near nominal, with the lower recall the previous figures show. Combined with its lack of an honest error guarantee even where it is well within the nominal rate, this is why the mixed-effects modes are the recommended default despite the fixed mode’s superficially higher raw power.
Which genes are detected? Performance by gene category
The planted effect is spread across three biologically distinct kinds of gene (see Part 1): cell-type markers (sparse, expressed only in the index type), housekeeping genes (abundant, expressed everywhere) and background genes (near-silent — the “random” category). Because detection power depends on how much a gene is expressed, splitting performance by category shows which genes a study of a given size can actually recover. This section covers all three continuous knobs the power family sweeps — patient samples, cells per sample, and effect size — each split by category and inference mode.
Detection versus sample size
What is plotted. For the samples sweep, raw (uncorrected) power (top row) and recall at a 5 % FDR (bottom row) versus the number of patients, split into one column per gene category and coloured by inference mode. Higher is better; error bars are \(\pm 1\) standard error over the 40 replicates.
How it is computed. From the power family, but with the true genes partitioned by category: 300 markers, 300 housekeeping and 300 background genes each carry the planted effect. Raw power and recall are computed exactly as before (uncorrected \(p < 0.05\); survival of the hierarchical gate at 5 % FDR) but averaged within each category, over the 40 replicates.
Interpretation. Detection tracks expression level. Housekeeping genes, which are abundant and broadly detected, are recovered with the fewest patients — their raw power is already high (~0.68) by 10 patients, and recall clears a substantial detection level by ~24 patients. Cell-type markers, expressed in only the index cells and so sampled at far lower depth, need more patients to reach the same raw power, only crossing that same level at ~NA patients. Background genes, planted on a near-silent baseline, are the hardest of all and stay barely recovered even at the largest sample size studied here (recall ~0 at 30 patients). The ordering is the same for all three inference modes — and, as elsewhere, the mixed-effects modes are uniformly more powerful than the fixed mode. The practical lesson is that a spatial study’s effective sensitivity is gene-specific: an effect of the same true size is far easier to detect on an abundant gene than on a sparse marker, so power and sample-size planning should be framed for the expression regime of the genes of interest.
Does buying more cells or a bigger effect change the ranking?
What is plotted. Recall at a 5 % FDR, split by gene
category (columns) and inference mode (colour), against cells per sample
(left) and effect size (right), at the fixed S = 24 used
for those two sweeps. Error bars are \(\pm
1\) standard error.
Interpretation. At the fixed S = 24
used for these two sweeps, background genes are
recovered only marginally even at the favourable end of each axis
(intercept-mode recall reaches just ~0.01 at 1,000 cells, and ~0.01 at
\(\beta =\) 3).
Housekeeping genes are detected readily across almost
the whole range of both axes. Marker genes — sparse,
but no longer starved of patients at this S — now show
real, climbing recall on both axes too, rather than staying flat at zero
the way they do at a smaller S: intercept-mode marker
recall reaches ~0.35 at 1,000 cells and ~0.19 at \(\beta =\) 3. The ordering — housekeeping
easiest, then markers, then background — holds throughout; what a larger
S changes is how much of each axis’s range clears the
detection threshold at all, not the ordering itself. As elsewhere, treat
any point where the fixed mode’s raw recall exceeds the
mixed-effects modes’ with caution: the FDR-control figure above shows
the fixed mode can be anti-conservative, so a higher recall there is not
necessarily a real advantage — the realised FDP, not raw recall, is what
determines whether a set of discoveries can be trusted. The
takeaway is unchanged: more patients is the more reliable lever than
more cells or a bigger effect for recovering sparse gene classes without
sacrificing FDR control.
Robustness to the spatial arrangement of cells
The niche effect is recovered under structured layouts and must vanish under the spatially random control.
Overview: pooled across gene categories, intercept mode
What is plotted. Two metrics per layout.
Recall @ 5 % FDR is the fraction of responsive genes
recovered at a 5 % false-discovery threshold. Niche
specificity is the fraction of the time the correct
niche cell type (B) is identified as the strongest
association for a responsive gene — a check that the method attributes
the effect to the right neighbour, not a bystander. Both range from 0 to
1; higher is better, except that under the random
negative control recall should be near zero.
How it is computed. From the power family’s
layout sweep: one dataset per layout (S =
24, 300 cells, β = 2, intercept mode) is fit. Recall @ 5 %
FDR uses the hierarchical gate; niche specificity asks, for each true
gene, whether the B-niche coefficient has the largest \(|t|\) among all candidate niche cell types.
Both are averaged over true genes and 40 replicates.
Interpretation. The gradient layout is the easiest
(a smooth, global signal) and gives the highest recall and specificity;
the patchy clustered and multiniche layouts are harder but still recover
the effect and, crucially, still point at the correct niche type. The
random control behaves as it must: recall collapses to
essentially zero, confirming that the method does not manufacture a
niche effect where there is no spatial structure. Because this overview
pools across gene categories, and housekeeping genes are by far the most
detectable at this study’s fixed S = 24, the figure above
is effectively telling the housekeeping story — the breakdown below
shows why.
Breakdown by gene category and inference mode
What is plotted. The same two metrics, now split by gene category (columns) and inference mode (fill/dodge), for every layout. Error bars are \(\pm 1\) standard error.
Interpretation. At this study’s fixed
S = 24, housekeeping genes are detected
often enough under every structured layout for the breakdown to be
informative, and — unlike at a smaller S —
marker genes now clear the bar too, most clearly under
the gradient layout (intercept-mode recall ~0.03). Recall collapses
under the random control for both, as it must. Niche specificity for
housekeeping and marker genes is high (~NA–NA) wherever there is real
spatial structure to attribute the effect to, versus close to chance
(~NA–NA) under the random control. Background genes
remain essentially undiscovered at every layout and mode (recall ~0.01
at best), so their niche-specificity bars, computed only over the tiny
handful of genes that do clear the gate, are noisy and sit close to
chance regardless of layout — they should be read as “not enough
discoveries to tell” rather than as evidence the method mis-attributes
the niche for this category. This is a direct consequence of the
sample-size threshold established earlier, not a layout-specific
weakness: markers only joined housekeeping in being reliably discovered
once S reached the level used here, and background genes
are expected to follow the same pattern at a large enough
S.
Computational scalability
What is plotted. Wall-clock time for the two stages of the pipeline — fitting the model and running inference — as the number of cells grows, at a realistic ~10,000-gene scale. Lower is better.
How it is computed. From the power family’s
cells sweep: at ~10,000 genes and S = 24
(intercept mode), the fit and inference stages are timed separately as
the number of cells per sample grows, averaged over replicates.
Interpretation. Fitting dominates the cost and grows sub-linearly with the number of cells; inference is cheaper and is the only stage that is chunked and parallelised across genes. At tens of thousands of cells and ten thousand genes the whole pipeline completes in minutes on a single core, and inference parallelises further across cores.
Summary
- The
spiDEsimulator reproduces the defining features of real spatial data: many genes, counts dominated by zeros and ones, a downward mean–variance trend, three biologically distinct gene categories, four spatial niche layouts, and a planted, condition-dependent niche effect that serves as ground truth. - Under this realistic model the mixed-effects modes control both type-I error and the false-discovery rate, where the fixed-effects mode is anti-conservative — most visibly when patients are few.
- Raw per-gene power rises smoothly with patients, cells and effect size, but recall shows a threshold: it stays near zero until the signal is strong enough to clear genome-wide multiplicity control, then climbs sharply. A high raw power is necessary but not sufficient for a discovery, which is why power can look strong while recall is still low.
- Detection is gene-category-dependent: abundant
housekeeping genes are recovered with far fewer patients, cells or a
smaller effect than sparse cell-type markers, and near-silent background
genes need the most of all — at this study’s design point
(
S =24 for the cells/effect/layout sweeps), marker genes are now meaningfully discovered while background genes remain essentially undetected, underscoring that patient count, more than cells or effect size, is what unlocks detection for sparse gene classes. - The niche effect is recovered across spatial layouts with correct niche attribution, and correctly vanishes under a spatially random control — though, consistent with the point above, this is clearly visible for housekeeping and, at this study’s sample size, marker genes; background genes still need more patients before their layout-specific recall becomes interpretable.
- The
CellType:conditiondesign term fixes what the niche layer estimates: without it a flat cell-type response loads onto the niche slope and its standard error shrinks to a regression through the origin, which reports more niche calls of a broader quantity; with it, the estimand is the within-cell-type gradient. The 0.99.17 defaults cost a mild liberal shift on this null, which plants no composition effect, and are justified on the real cohort. - The two-stage estimator, run paired on the same
simulated datasets, is the best-calibrated arm in this simulation study
(at or under nominal at every null design point for its
ols/nbstage-1 arms) but delivers less power than the mixed modes at every design point measured — and on the real cohort the niche-shuffle null reverses the calibration ordering; the dedicated vignette Two-stage estimation: benchmarked against the published simulation study reports its full results, including the adversarial scenarios the published simulator cannot express. - The pipeline scales to ten thousand genes and tens of thousands of cells in minutes.
The other reports on this site take up the parts this study cannot: The real cohort (the null tail, its cause and the fix, measured on real data), The two-stage estimator, Combining correlated niche p-values, Gene-set inference, and What was tried and rejected.
The complete study — the simulator, the scenario definitions, and
SLURM scripts for a compute cluster — is in the package’s
research/ directory. To regenerate the summary tables used
above at full scale, run the study (see research/README.md)
and then Rscript research/R/install_results.R. This
vignette reads those tables directly, so rebuilding it (no
precompilation step) picks up the refreshed results.
#> R version 4.5.2 (2025-10-31)
#> Platform: x86_64-conda-linux-gnu
#> Running under: Red Hat Enterprise Linux 9.8 (Plow)
#>
#> Matrix products: default
#> BLAS/LAPACK: /home/uqdbhuva/miniconda3/envs/latest-r/lib/libopenblasp-r0.3.30.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Australia/Brisbane
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] patchwork_1.3.2 ggplot2_4.0.1
#> [3] SpatialExperiment_1.20.0 SingleCellExperiment_1.32.0
#> [5] SummarizedExperiment_1.40.0 Biobase_2.70.0
#> [7] GenomicRanges_1.62.0 Seqinfo_1.0.0
#> [9] IRanges_2.44.0 S4Vectors_0.48.0
#> [11] BiocGenerics_0.56.0 generics_0.1.4
#> [13] MatrixGenerics_1.22.0 matrixStats_1.5.0
#> [15] spiDE_0.99.18
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 SparseArray_1.10.1 lattice_0.22-9
#> [4] digest_0.6.38 magrittr_2.0.4 evaluate_1.0.5
#> [7] grid_4.5.2 RColorBrewer_1.1-3 fastmap_1.2.0
#> [10] jsonlite_2.0.0 Matrix_1.7-4 limma_3.66.0
#> [13] viridisLite_0.4.2 scales_1.4.0 codetools_0.2-20
#> [16] jquerylib_0.1.4 abind_1.4-8 cli_3.6.5
#> [19] rlang_1.1.6 XVector_0.50.0 withr_3.0.2
#> [22] cachem_1.1.0 DelayedArray_0.36.0 yaml_2.3.10
#> [25] S4Arrays_1.10.0 tools_4.5.2 parallel_4.5.2
#> [28] BiocParallel_1.44.0 dplyr_1.1.4 vctrs_0.6.5
#> [31] R6_2.6.1 lifecycle_1.0.4 magick_2.9.0
#> [34] pkgconfig_2.0.3 pillar_1.11.1 bslib_0.9.0
#> [37] gtable_0.3.6 glue_1.8.0 Rcpp_1.1.2
#> [40] statmod_1.5.1 tidyselect_1.2.1 tibble_3.3.0
#> [43] xfun_0.54 knitr_1.50 farver_2.1.2
#> [46] rjson_0.2.23 htmltools_0.5.8.1 labeling_0.4.3
#> [49] rmarkdown_2.30 compiler_4.5.2 S7_0.2.1