Neighbourhood-dependent differential expression with spiDE
Dharmesh D. Bhuva
Source:vignettes/spiDE.Rmd
spiDE.RmdIntroduction
spiDE finds neighbourhood-dependent differential
expression in spatial transcriptomics. Within an index
cell type, it tests whether a gene’s response to an experimental
condition changes with the local density (the niche)
of a surrounding cell type. One hypothesis is one (gene, index type,
niche type) triplet.
Three stages:
-
Niches —
buildNiches(): per sample and per bandwidth, a Gaussian kernel density of every cell type evaluated at every cell. -
Fit —
fitSpiDE(): one negative binomial GLM per gene over a design withcell type : condition : nicheinteractions, patient-level random effects, and a per-gene convergence step. -
Test —
testSpiDE(): per-gene Wald statistics, combined across a gene’s correlated niches and across bandwidths, then a hierarchical FDR.
spiDE() runs all three. The model is stated in full in
vignette("spiDE-model"); every benchmark number quoted
about it is in vignette("spiDE-calibration").
Example data
The package ships a small synthetic SpatialExperiment,
toySpiDE: 20 genes, 480 cells, 6 samples (3 Responders, 3
Non-responders), three cell types. Gene G1 is planted
up-regulated in index type A, in Responders, in
proportion to the local density of niche type B.
data(toySpiDE)
toySpiDE
#> class: SpatialExperiment
#> dim: 20 480
#> metadata(0):
#> assays(1): counts
#> rownames(20): G1 G2 ... G19 G20
#> rowData names(0):
#> colnames(480): cell1 cell2 ... cell479 cell480
#> colData names(9): sample_id condition ... Area nCount
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x y
#> imgData names(0):
table(toySpiDE$cell_type, toySpiDE$condition)
#>
#> Non-responder Responder
#> A 85 92
#> B 74 66
#> C 81 82One call
res <- spiDE(toySpiDE, condition = "condition", sigma = c(10, 30, 50),
verbose = FALSE)
res
#> SpiDEResults
#> Bandwidths (sigma): 10, 30, 50
#> Genes: 20
#> Condition: condition
#> Index cell types: A, B, C
#> Niche cell types: A, B, C
#> Tested: TRUE (1 rows in results table)results() returns the tidy table, keyed by (gene, index
type, niche type) with the best bandwidth:
tab <- results(res)
head(tab[order(tab$fdr.niche), ], 10)
#> gene ct_index ct_niche bandwidth.max coef t DirectionGene
#> 1 G1 A B 50 1.588429 6.737216 Up
#> DirectionIndex DirectionNiche fdr.gene fdr.index fdr.niche
#> 1 Up Up 5.879288e-09 2.939642e-10 9.798828e-11The planted effect is recovered, with a positive direction:
tab[tab$gene == "G1" & tab$ct_index == "A" & tab$ct_niche == "B", ]
#> gene ct_index ct_niche bandwidth.max coef t DirectionGene
#> 1 G1 A B 50 1.588429 6.737216 Up
#> DirectionIndex DirectionNiche fdr.gene fdr.index fdr.niche
#> 1 Up Up 5.879288e-09 2.939642e-10 9.798828e-11The result layers, and what each estimates
One fit answers three questions, each a different quantity:
vapply(c("niche", "celltype", "patient"),
function(ty) nrow(results(res, type = ty)), integer(1))
#> niche celltype patient
#> 1 0 0-
"niche"(default): does the neighbourhood modulate the response within an index cell type? The estimate is the gradient of the condition effect along the niche density, within (sample, cell type). Keyed by (gene,ct_index,ct_niche). -
"celltype": does the response differ by cell type, regardless of neighbourhood? The flat, cell-type-wide shift, keyed by (gene,ct_index), with a two-level FDR. -
"patient": does the gene respond at the patient level at all? One test per gene.
ct <- results(res, type = "celltype")
if (nrow(ct)) head(ct[order(ct$fdr.celltype), ], 5) else "no calls on this toy data"
#> [1] "no calls on this toy data"The layers are not redundant: a flat response is not a neighbourhood
effect and is reported separately rather than allowed to load onto the
niche slope. A fourth question — does a patient’s composition,
the mean niche density around their index cells, associate with
expression in those cells — is patient-level and is asked by
compositionTest(), below.
Step by step
The wrapper equals the three stages run explicitly, which is useful to inspect intermediate objects, reuse a fit, or change one stage.
Build niches
buildNiches() adds one reducedDim per
bandwidth. mergeNiches() coarsens cell types into niche
groups; computeSizeFactors() derives a per-sample
library-size offset column.
spe <- buildNiches(toySpiDE, sigma = c(10, 30, 50))
reducedDimNames(spe)
#> [1] "Niche10" "Niche30" "Niche50"Fit
fit <- fitSpiDE(spe, condition = "condition", verbose = FALSE)
fit
#> SpiDEResults
#> Bandwidths (sigma): 10, 30, 50
#> Genes: 20
#> Condition: condition
#> Index cell types: A, B, C
#> Niche cell types: A, B, C
#> Tested: FALSE (0 rows in results table)Each bandwidth’s fit is a SpiDEFit with the coefficients
alpha, dispersions psi, the design
W, a tag per design column, and — for the default mixed,
converged fit — the variance components tau2 and per-gene
convergence diagnostics polish:
f1 <- fits(fit)[[1]]
table(f1$covtype)
#>
#> CellType Niche Response ResponseNiche
#> 3 4 0 6
#> ResponseCellType Other Random
#> 3 3 24
f1$tau2
#> SampleInt SampleCellTypeInt
#> 0.002286559 0.005124837A fit made without the convergence step
(converge = FALSE), or saved by an older version, can be
converged afterwards with polishSpiDE(); it returns the
object with converged coefficients and dispersions and inference
cleared, so testSpiDE() runs again:
fit0 <- fitSpiDE(spe, condition = "condition", converge = FALSE, verbose = FALSE)
fit1 <- polishSpiDE(fit0, spe, verbose = FALSE)The patient-level question
compositionTest() asks the composition question on
patients: per (sample, index type) a pseudobulk profile of the index
cells and the mean niche density around them, then a limma
moderated
across samples. Two terms are reported per (gene, index, niche):
"niche", the association pooled across conditions, and
"condition:niche", its difference between conditions.
ct <- compositionTest(spe, condition = "condition", sigma = 30, verbose = FALSE)
head(ct[order(ct$p), c("gene", "ct_index", "ct_niche", "term", "t", "fdr")])
#> gene ct_index ct_niche term t fdr
#> 205 G5 C B niche 3.012393 0.08958363
#> 165 G5 C A niche -2.866971 0.13161462
#> 131 G11 B C niche 2.703142 0.20066920
#> 144 G4 B C condition:niche -2.679210 0.13990485
#> 35 G15 A B condition:niche -4.241351 0.23139488
#> 72 G12 A C condition:niche 4.306386 0.26230412On a fixture this small nothing is expected to pass. The point is that the within-sample niche question and the between-patient composition question have two tests, and neither is reported as the other.
Niche-only analysis
Not every question has a condition. With
condition = NULL the same framework tests how expression
within an index type changes with the density of a niche type, in every
sample alike:
res_niche <- spiDE(toySpiDE, condition = NULL, sigma = c(10, 30),
verbose = FALSE)
head(results(res_niche))
#> gene ct_index ct_niche bandwidth.max coef t DirectionGene
#> 1 G1 A B 30 0.9247049 3.664604 Up
#> DirectionIndex DirectionNiche fdr.gene fdr.index fdr.niche
#> 1 Up Up 0.03318263 0.001657729 0.0005522458The tested quantity is then the CellType:niche slope.
The "celltype" and "patient" layers are empty
in this mode, and spiGSEA(type = "celltype") errors, there
being no condition to contrast. The mode is mildly anti-conservative
because the niche covariate is spatially autocorrelated; see the model
vignette.
Scaling
The fit sees all genes at once (the dispersion is moderated across genes before the per-gene convergence step). Inference and convergence are independent across genes, processed in gene blocks and dispatched with BiocParallel; the counts may be a DelayedArray, realised one block at a time. On the CPU backend, results are exactly invariant to the block size and the number of workers.
library(BiocParallel)
res <- spiDE(toySpiDE, condition = "condition",
block.size = 500, BPPARAM = MulticoreParam(4))backend = c("auto", "cpu", "gpu") selects the compute
backend for both stages. On the GPU the per-gene Wald covariance and
working weights are batched across each gene block on the device, the
block size is chosen to fit the detected memory (override with
gpu.mem.budget, in bytes), and a serial
BPPARAM is used since execution is already batched. The GPU
path uses single precision on Metal, so it matches the CPU path to a
tolerance rather than exactly.
The defaults
What spiDE() and fitSpiDE() do unless told
otherwise, one line each. The reasons are in
vignette("spiDE-model") and the measurements in
vignette("spiDE-calibration").
| argument | default | what it does |
|---|---|---|
random |
"intercept" |
a random intercept per sample, so response effects are tested against patient-level variation |
re.celltype |
TRUE |
a random intercept per (sample, cell type), so the niche slope is a within-group slope and does not carry a patient’s composition |
converge |
TRUE |
converge each gene to its own optimum after the shared fit, and re-estimate its dispersion there; needs integer counts |
df.method |
"satterthwaite" |
a reference df per tested coefficient |
combine |
"cauchy" |
correlation-agnostic combination of a gene’s niche p-values |
A sample-level covariate (Age, Stage) is collinear with the
per-sample random intercept and is rejected; adjust for cell-level
covariates, including per-cell library size, through
covariates.
Session info
sessionInfo()
#> 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] SpatialExperiment_1.22.0 SingleCellExperiment_1.34.0
#> [3] SummarizedExperiment_1.42.0 Biobase_2.72.0
#> [5] GenomicRanges_1.64.0 Seqinfo_1.2.0
#> [7] IRanges_2.46.0 S4Vectors_0.50.2
#> [9] BiocGenerics_0.58.1 generics_0.1.4
#> [11] MatrixGenerics_1.24.0 matrixStats_1.5.0
#> [13] spiDE_0.99.17 BiocStyle_2.40.0
#>
#> loaded via a namespace (and not attached):
#> [1] fastmap_1.2.0 spatstat.geom_3.8-2 spatstat.explore_3.8-2
#> [4] digest_0.6.39 rsvd_1.0.5 lifecycle_1.0.5
#> [7] statmod_1.5.2 processx_3.9.0 spatstat.data_3.1-9
#> [10] magrittr_2.0.5 compiler_4.6.1 rlang_1.3.0
#> [13] sass_0.4.10 tools_4.6.1 yaml_2.3.12
#> [16] knitr_1.52 S4Arrays_1.12.0 bit_4.6.0
#> [19] DelayedArray_0.38.2 abind_1.4-8 BiocParallel_1.46.0
#> [22] desc_1.4.3 grid_4.6.1 polyclip_1.10-7
#> [25] beachmat_2.28.0 edgeR_4.10.4 spatstat.utils_3.2-4
#> [28] cli_3.6.6 rmarkdown_2.32 ragg_1.5.2
#> [31] otel_0.2.0 rjson_0.2.23 cachem_1.1.0
#> [34] splines_4.6.1 parallel_4.6.1 BiocManager_1.30.27
#> [37] XVector_0.52.0 coro_1.1.0 Matrix_1.7-5
#> [40] jsonlite_2.0.0 bookdown_0.48 BiocSingular_1.28.0
#> [43] callr_3.8.0 bit64_4.8.6 irlba_2.3.7
#> [46] SpaNorm_1.7.9 tensor_1.5.1 systemfonts_1.3.2
#> [49] magick_2.9.1 locfit_1.5-9.12 limma_3.68.5
#> [52] spatstat.univar_3.2-0 jquerylib_0.1.4 goftest_1.2-3
#> [55] spatstat.random_3.5-1 pkgdown_2.2.1 codetools_0.2-20
#> [58] ps_1.9.3 deldir_2.0-4 ScaledMatrix_1.20.0
#> [61] htmltools_0.5.9 torch_0.17.0 R6_2.6.1
#> [64] textshaping_1.0.5 evaluate_1.0.5 lattice_0.22-9
#> [67] bslib_0.12.0 Rcpp_1.1.2 SparseArray_1.12.2
#> [70] nlme_3.1-169 spatstat.sparse_3.2-0 xfun_0.60
#> [73] fs_2.1.0