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:

1. 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.
2. The benchmark. A large study that fits spiDE to 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 samples

What 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.

Distribution of individual gene-cell counts (left) and per-cell sequencing depth (right). Almost every measurement is 0, 1 or 2, as in imaging-based spatial data, yet total depth per cell is realistic.

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.

Empirical mean-variance trend. Each point is a gene; its biological coefficient of variation (BCV) falls as mean expression rises. The line is the model’s true trend, BCV = 0.3 + 1/sqrt(mean).

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.

Left: per-gene mean expression by category, spanning orders of magnitude. Right: detection rate (fraction of cells with a non-zero count) by category. Housekeeping genes are broadly detected; markers and background are sparse.

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.

Cell-type specificity of a sample of genes. Rows are genes grouped by category (marker rows are further split by the cell type they mark); columns are cell types; colour is the fraction of that type’s cells in which the gene is detected. Markers are detected in one type; housekeeping is detected everywhere; background is detected nowhere.

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 planted differential-expression effect. Each tile is a gene’s mean expression in the index cell type of one patient. Rows are grouped by gene category and DE status – including DE background genes – and the 24 patient columns are split into Responders and Non-responders. DE genes are elevated specifically in Responders; everything else is flat.

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 B increases 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 — B cells 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 B niche (patchy) coexists with a second, independent structured cell type D occupying different hotspots. This is a specificity test: the method must attribute the effect to B and not be misled by the competing structure of D.

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.

Three views of one sample per spatial scenario, top to bottom: cell-type layout; the true (noise-free) planted signal in index (A) cells; and the noisy simulated log-counts of one example marker gene after passing that signal through the negative-binomial count model. The bottom row is what spiDE actually sees.

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 being stage1 = "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.

Realised type-I error at alpha = 0.05 under the null, versus the number of patient samples, for two layouts. Colour is the inference mode; linetype is the reference-df method (between-patient S-2 vs the per-coefficient Satterthwaite approximation). On the dashed line is calibrated; above is anti-conservative; below is conservative.

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.

QQ plot of null p-values at the largest sample size, split by layout. On the diagonal is calibrated; above it is anti-conservative; below it is 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:

Null type-I error of the response-niche test by sample size, for the niche-only design and the design with CellType:condition restored (gradient layout, satterthwaite df). Restoring the term costs ~k degrees of freedom out of several hundred, so calibration should be unchanged; a visible gap would mean the added block is absorbing signal rather than nuisance.

The arm without the term reports more niche calls, and the difference is large:

True-positive rate at FDR 0.05 versus effect size, by model design. The added term spends degrees of freedom, so a small power cost is expected; a large one would argue against restoring it.

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.

Null type-I error of the response-niche test by sample size for the shipped arm and the three 0.99.17 switch arms (intercept mode, satterthwaite df; 40 replicates, ribbon = 95% interval of the mean). Both switches cost, with different S-dependence: convergence at small S, the nested block as a constant offset from S = 10.

Recall (TPR) and realised FDP at nominal FDR 0.05 in the samples sweep, shipped arm against the 0.99.17 defaults.

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 — will be judged against.

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).

Raw (uncorrected) power to detect the planted niche effect versus the number of patient samples, the number of cells per sample, and the effect size. Higher is better. Colour is the inference mode; linetype is the reference-df method (between vs Satterthwaite) – the fixed and two-stage modes have no df.method and always appear solid.

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 and between becomes less punitive. This is the power counterpart of the calibration figure above, and the reason Satterthwaite is the default — the conservatism between buys at small $S$ is paid for almost entirely in lost detections. Note the quoted figures in the two bullets below are for the between arm.

• 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.

Raw (uncorrected) power rises smoothly with sample size, but recall after genome-wide FDR control stays near zero until the study is large enough, then rises sharply – a detection threshold imposed by multiplicity correction.

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.

Recall at 5% FDR versus samples, cells per sample, and effect size. The same detection threshold appears on every axis: recall is near zero until the signal is strong enough, then rises. Colour is the inference mode; linetype is the reference-df method (between vs Satterthwaite) – the fixed and two-stage modes have no df.method and always appear solid.

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.

Realised FDP at a nominal 5% FDR versus samples, cells per sample and effect size, shown only where a mode makes appreciable discoveries. The fixed mode is anti-conservative when discoveries are scarce; the mixed modes control FDP wherever they discover.

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.

Detection split by gene category (columns) and inference mode (colour), for raw power (top) and recall at 5% FDR (bottom), versus sample size. Abundant housekeeping genes are recovered with far fewer patients than sparse markers; near-silent background genes are hardest.

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.

Recall at 5% FDR by gene category and inference mode, versus cells per sample (left) and effect size (right), at S = 24. Housekeeping genes are recovered first on both axes; markers lag; background genes lag furthest of all.

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.

Recall at 5% FDR and niche specificity across the four layouts. Structured layouts (gradient, clustered, multiniche) recover the effect and identify the correct niche; the random control correctly yields near-zero recall.

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.

Recall at 5% FDR and niche specificity across the four layouts, split by gene category (columns) and arm (bars), at S = 24. The twostage arm appears only in the recall row: its results table carries no per-gene strongest-niche ranking, so it has no niche-specificity metric.

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.

Fit and inference time versus the number of cells, at ~10,000 genes. Both stages scale gracefully to tens of thousands of cells.

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 spiDE simulator 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:condition design 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/nb stage-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.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] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] patchwork_1.3.2             ggplot2_4.0.3              
#>  [3] SpatialExperiment_1.22.0    SingleCellExperiment_1.34.0
#>  [5] SummarizedExperiment_1.42.0 Biobase_2.72.0             
#>  [7] GenomicRanges_1.64.0        Seqinfo_1.2.0              
#>  [9] IRanges_2.46.0              S4Vectors_0.50.2           
#> [11] BiocGenerics_0.58.1         generics_0.1.4             
#> [13] MatrixGenerics_1.24.0       matrixStats_1.5.0          
#> [15] spiDE_0.99.17               BiocStyle_2.40.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        rjson_0.2.23        xfun_0.60          
#>  [4] bslib_0.12.0        lattice_0.22-9      vctrs_0.7.3        
#>  [7] tools_4.6.1         parallel_4.6.1      tibble_3.3.1       
#> [10] pkgconfig_2.0.3     Matrix_1.7-5        RColorBrewer_1.1-3 
#> [13] S7_0.2.2            desc_1.4.3          lifecycle_1.0.5    
#> [16] compiler_4.6.1      farver_2.1.2        textshaping_1.0.5  
#> [19] statmod_1.5.2       codetools_0.2-20    htmltools_0.5.9    
#> [22] sass_0.4.10         yaml_2.3.12         pkgdown_2.2.1      
#> [25] pillar_1.11.1       jquerylib_0.1.4     BiocParallel_1.46.0
#> [28] DelayedArray_0.38.2 cachem_1.1.0        limma_3.68.5       
#> [31] magick_2.9.1        abind_1.4-8         tidyselect_1.2.1   
#> [34] digest_0.6.39       dplyr_1.2.1         bookdown_0.48      
#> [37] labeling_0.4.3      fastmap_1.2.0       grid_4.6.1         
#> [40] cli_3.6.6           SparseArray_1.12.2  magrittr_2.0.5     
#> [43] S4Arrays_1.12.0     withr_3.0.3         scales_1.4.0       
#> [46] rmarkdown_2.32      XVector_0.52.0      otel_0.2.0         
#> [49] ragg_1.5.2          evaluate_1.0.5      knitr_1.52         
#> [52] viridisLite_0.4.3   rlang_1.3.0         Rcpp_1.1.2         
#> [55] glue_1.8.1          BiocManager_1.30.27 jsonlite_2.0.0     
#> [58] R6_2.6.1            systemfonts_1.3.2   fs_2.1.0