| Title: | Context-specific, neighbourhood-dependent differential expression |
|---|---|
| Description: | spiDE identifies context-specific, neighbourhood-dependent differential expression in spatial transcriptomics data. Within an index cell type, it asks whether a gene's expression changes with the local density (niche) of another cell type within patients, and whether that dependence differs between two conditions. Niche covariates are Gaussian kernel density estimates of each cell type. For every index cell type a per-gene negative binomial GLM with a patient intercept is fitted with SpaNorm's Newton solver; its niche slopes are then tested across patients, either from each patient's own slopes with weighted limma (the slopes engine) or with a patient-clustered CR2 sandwich (the sandwich engine). The pooled test asks whether a slope is consistently non-zero; the condition-specific test whether it differs between conditions. |
| Authors: | Dharmesh D. Bhuva [aut, cre] (ORCID: <https://orcid.org/0000-0002-6398-9157>), Chin Wee Tan [aut] (ORCID: <https://orcid.org/0000-0001-9695-7218>) |
| Maintainer: | Dharmesh D. Bhuva <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.99.38 |
| Built: | 2026-11-11 01:10:15 UTC |
| Source: | https://github.com/BiocStaging/spiDE |
Computes, for each biological sample and each bandwidth sigma, a
Gaussian kernel density estimate of every cell type evaluated at each cell's
location. The resulting cells x cellTypes matrices are stored as
reducedDim(spe, "Niche<sigma>"), one per bandwidth, with a consistent
set of cell-type columns (cell types absent from a sample are set to 0).
buildNiches(spe, ...) ## S4 method for signature 'ANY' buildNiches( spe, sigma = c(10, 30, 50, 70), cell_type = "cell_type", sample_id = "sample_id", edge = TRUE, diggle = TRUE, name = "Niche", BPPARAM = BiocParallel::SerialParam(), ... )buildNiches(spe, ...) ## S4 method for signature 'ANY' buildNiches( spe, sigma = c(10, 30, 50, 70), cell_type = "cell_type", sample_id = "sample_id", edge = TRUE, diggle = TRUE, name = "Niche", BPPARAM = BiocParallel::SerialParam(), ... )
spe |
a SpatialExperiment with spatial coordinates, cell type labels, and sample identifiers. |
... |
ignored. |
sigma |
a numeric vector of kernel bandwidths (in the units of
|
cell_type |
a character, the colData column holding cell type labels. |
sample_id |
a character, the colData column identifying samples. |
edge, diggle
|
logicals passed to |
name |
a character, the prefix for the stored reducedDims (default "Niche"). |
BPPARAM |
a BiocParallelParam for parallelising over samples. |
the input spe with one reducedDim per bandwidth added.
data(toySpiDE) spe <- toySpiDE spe <- buildNiches(spe, sigma = c(10, 20)) SingleCellExperiment::reducedDimNames(spe)data(toySpiDE) spe <- toySpiDE spe <- buildNiches(spe, sigma = c(10, 20)) SingleCellExperiment::reducedDimNames(spe)
Fits, for every index cell type, one negative binomial GLM per gene of its
expression on the local density of every other (niche) cell type, with an
intercept per patient so that each niche effect is a within-patient slope.
The slopes are then tested across patients by testSpiDE(), pooled or
between conditions.
fitSpiDE(spe, condition = NULL, ...) ## S4 method for signature 'SpatialExperiment' fitSpiDE( spe, condition = NULL, engine = c("slopes", "sandwich"), index = NULL, niche = NULL, sigma = NULL, cell_type = "cell_type", sample_id = "sample_id", section = NULL, covariates = character(), depth = c("loglib", "nonlinear", "spatial_spline", "none"), offset = NULL, strata = NULL, assay = "counts", name = "Niche", min.cells = 10L, min.patients = 6L, min.detect = 0.1, genes = NULL, tile = NULL, min.tiles = 5L, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE, ... )fitSpiDE(spe, condition = NULL, ...) ## S4 method for signature 'SpatialExperiment' fitSpiDE( spe, condition = NULL, engine = c("slopes", "sandwich"), index = NULL, niche = NULL, sigma = NULL, cell_type = "cell_type", sample_id = "sample_id", section = NULL, covariates = character(), depth = c("loglib", "nonlinear", "spatial_spline", "none"), offset = NULL, strata = NULL, assay = "counts", name = "Niche", min.cells = 10L, min.patients = 6L, min.detect = 0.1, genes = NULL, tile = NULL, min.tiles = 5L, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE, ... )
spe |
a SpatialExperiment with raw counts and a niche reducedDim from
|
condition |
|
... |
unused. |
engine |
|
index |
|
niche |
|
sigma |
the niche bandwidth (needed only if |
cell_type, sample_id
|
the colData columns of the cell type and of the patient (the unit that is replicated: one sample per patient). |
section |
|
covariates |
a character vector of cell-level colData columns to adjust for (centred; patient-level ones are absorbed by the patient intercepts). |
depth |
how sequencing depth enters: |
offset |
|
strata |
|
assay |
the counts assay (raw integer counts). |
name |
the niche reducedDim prefix. |
min.cells |
minimum cells of an index type a patient must contribute. |
min.patients |
minimum patients an index type needs. |
min.detect |
minimum fraction of an index type's cells in which a gene is detected for it to be tested in that type. |
genes |
|
tile |
side of the square tiles of the spatial variance estimate, in coordinate units (default 3 x the bandwidth). |
min.tiles |
minimum tiles a patient needs for its own spatial variance. |
BPPARAM |
a BiocParallelParam; genes are fitted in parallel. |
verbose |
report progress. |
Two engines estimate the between-patient error the tests need:
"slopes" (default)fits the shared model once, without the
condition, then gives every patient a one-step NB estimate of its own
niche slopes with a within-patient spatial sandwich variance.
testSpiDE() combines them across patients with weighted limma. It is
the more powerful engine when effects vary between patients, and its
fit never sees the condition, so any condition can be tested later.
"sandwich"fits the condition-specific model directly (the condition x niche columns) and tests its coefficients with a patient-clustered CR2 sandwich and Bell-McCaffrey degrees of freedom. It is simpler, and strongest for effects shared by all patients, but needs the condition at fit time.
The slopes engine passes every calibration null of five comparisons in four spatial cohorts, the sandwich engine all but one borderline tail; see the calibration vignette.
a SpiDEFit.
testSpiDE(), spiDE(), testNicheAbundance()
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) fitdata(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) fit
Sums the density columns of grouped cell types in each niche
reducedDim, producing coarser niche covariates (e.g. merging
Macrophage subtypes into a single "Macrophage" niche). Cell types absent from
groups are retained as their own (singleton) niche.
mergeNiches(spe, groups, ...) ## S4 method for signature 'ANY' mergeNiches(spe, groups, sigma = NULL, name = "Niche", ...)mergeNiches(spe, groups, ...) ## S4 method for signature 'ANY' mergeNiches(spe, groups, sigma = NULL, name = "Niche", ...)
spe |
a SpatialExperiment with niche reducedDims (see |
groups |
a named list mapping each merged niche name to a character vector of cell-type column names to sum. |
... |
ignored. |
sigma |
a numeric vector of bandwidths to update; |
name |
a character, the reducedDim name prefix (default "Niche"). |
The group membership is recorded in metadata(spe)$spiDE_niche_groups
(per niche reducedDim) so that a subsequent fitSpiDE() never tests an
index cell type against a merged niche it is a member of (its own density
would otherwise contaminate the niche).
the input spe with merged niche reducedDims.
data(toySpiDE) spe <- toySpiDE spe <- buildNiches(spe, sigma = 30) spe <- mergeNiches(spe, groups = list(AC = c("A", "C")), sigma = 30) colnames(SingleCellExperiment::reducedDim(spe, "Niche30"))data(toySpiDE) spe <- toySpiDE spe <- buildNiches(spe, sigma = 30) spe <- mergeNiches(spe, groups = list(AC = c("A", "C")), sigma = 30) colnames(SingleCellExperiment::reducedDim(spe, "Niche30"))
Every index type's shared fit has an intercept per patient and gene: the
patient's expression level of the gene in that type, at zero niche density
and a reference depth: the index type's mean log library size under
depth = "loglib" or "nonlinear", each patient's own mean
under "spatial_spline". The intercepts absorb
every difference between patients – the condition's main effect,
composition, batch – which is why the niche slopes are within-patient
slopes. plotPatientEffects() shows what they capture. Kept by
fitSpiDE() from spiDE 0.99.34 on, for both engines.
patientIntercepts(object, ...) ## S4 method for signature 'SpiDEFit' patientIntercepts(object, gene = NULL, index = NULL, ...) ## S4 method for signature 'SpiDEResults' patientIntercepts(object, ...)patientIntercepts(object, ...) ## S4 method for signature 'SpiDEFit' patientIntercepts(object, gene = NULL, index = NULL, ...) ## S4 method for signature 'SpiDEResults' patientIntercepts(object, ...)
object |
a SpiDEFit or SpiDEResults. |
... |
unused. |
gene, index
|
optional character vectors to subset. |
a data.frame: gene, index, patient,
intercept (log scale), ncells (the patient's cells of the
index type), plus the patient-level colData columns recorded by
fitSpiDE().
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30, verbose = FALSE) head(patientIntercepts(fit, gene = "G1"))data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30, verbose = FALSE) head(patientIntercepts(fit, gene = "G1"))
The slopes engine estimates every patient's own niche slopes (a one-step negative binomial estimate from the shared fit); this returns them in long form, e.g. to plot how a gene's slope varies between patients and conditions.
patientSlopes(object, ...) ## S4 method for signature 'SpiDEFit' patientSlopes(object, gene = NULL, index = NULL, niche = NULL, ...) ## S4 method for signature 'SpiDEResults' patientSlopes(object, ...)patientSlopes(object, ...) ## S4 method for signature 'SpiDEFit' patientSlopes(object, gene = NULL, index = NULL, niche = NULL, ...) ## S4 method for signature 'SpiDEResults' patientSlopes(object, ...)
object |
a SpiDEFit or SpiDEResults from the slopes engine. |
... |
unused. |
gene, index, niche
|
optional character vectors to subset. |
a data.frame: gene, index, niche,
patient, slope, var_spatial (the within-patient
spatial sandwich variance), var_model, ncells, plus the
patient-level colData columns recorded by fitSpiDE().
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) head(patientSlopes(fit, gene = "G1"))data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) head(patientSlopes(fit, gene = "G1"))
The number of genes called for every (index, niche) pair, split by direction: for the pooled test, a slope that rises or falls with the niche's density; for the condition-specific test, a slope higher in one condition than in the other. The heatmap gives exact counts and marks untested pairs (an index type against its own niche) in grey; the graph draws cell types on a fixed circle with an arrow from niche to index type, its width the number of genes – it reads best when calls are few (the condition test), the heatmap when every pair has some.
plotCallMap( x, test = NULL, fdr = x@fdr, style = c("heatmap", "graph"), min.calls = 1L )plotCallMap( x, test = NULL, fdr = x@fdr, style = c("heatmap", "graph"), min.calls = 1L )
x |
a SpiDEResults. |
test |
|
fdr |
the q-value threshold of a call. |
style |
|
min.calls |
the graph's smallest drawn edge. |
Read it for which cell-type pairs carry the signal, and for pairs whose
calls all go one way toward one niche type (the shape segmentation spillover
gives; see plotSpillover()). An arrow is an association, not signalling.
a ggplot. For the heatmap, p$data has one row per (index,
niche, direction): calls and tested.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, verbose = FALSE) plotCallMap(res, test = "pooled") plotCallMap(res, style = "graph")data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, verbose = FALSE) plotCallMap(res, test = "pooled") plotCallMap(res, style = "graph")
One sample's cells, coloured by the log1p density of a niche type: the
covariate spiDE regresses on. With index, only the index type's
cells are coloured – the covariate as the index type's model reads it –
and the others are dim. The panel is dark (the sticker's ink) so that the
lowest densities stay visible; the legend, titles and axes keep the page's
style.
plotNicheMap( spe, niche, sample, index = NULL, sigma = NULL, cell_type = "cell_type", sample_id = "sample_id", section = NULL, name = "Niche", point.size = NULL, scale.bar = TRUE, unit = NULL )plotNicheMap( spe, niche, sample, index = NULL, sigma = NULL, cell_type = "cell_type", sample_id = "sample_id", section = NULL, name = "Niche", point.size = NULL, scale.bar = TRUE, unit = NULL )
spe |
a SpatialExperiment with niche reducedDims ( |
niche |
a niche column (a cell type or a |
sample |
the sample (patient) to show. |
index |
|
sigma |
the bandwidth (needed only if |
cell_type, sample_id
|
the colData columns of cell type and sample. |
section |
|
name |
the niche reducedDim prefix. |
point.size |
the point size (default: from the number of cells). |
scale.bar |
draw a scale bar of about a fifth of the width. |
unit |
the coordinates' unit, for the scale bar's label (e.g.
|
Read it to judge the bandwidth: a speckled field that follows single cells is too small; a field nearly flat within a section is too large (the patient intercept then absorbs it). It also shows where the index cells sit relative to the niche.
a ggplot; p$data has one row per cell of the sample:
x, y (from each section's own origin), cell_type,
density (log1p),
section and focal (coloured).
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) plotNicheMap(spe, niche = "B", sample = "S1", index = "A")data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) plotNicheMap(spe, niche = "B", sample = "S1", index = "A")
The index type's cells' expression of each gene (counts per 10,000, log axis) against the log1p density of the niche type, binned, one panel per gene.
plotNicheResponse( spe, x, gene, index, niche, adjust = TRUE, bins = 8L, assay = "counts", name = "Niche" )plotNicheResponse( spe, x, gene, index, niche, adjust = TRUE, bins = 8L, assay = "counts", name = "Niche" )
spe |
the SpatialExperiment the fit was made from. |
x |
a SpiDEResults or SpiDEFit; |
gene |
one or more genes, one panel each. |
index, niche
|
the index and niche cell types. |
adjust |
divide out each patient's intercept (see Details). |
bins |
the number of density bins ( |
assay |
the counts assay. |
name |
the niche reducedDim prefix. |
With adjust = TRUE (the default) each cell's expression is divided by
its patient's fitted intercept, so the patients share one reference level
and overlay: faint points are each patient's mean in a density bin, points
and bars the mean across patients and its 95\
the dashed line each group's fitted slope. This is what the tests compare.
With adjust = FALSE the data are raw: each patient's own slope is a
segment through its mean (open circle), and the dashed line through the
patients' means is the between-patient association that the intercepts
absorb – testNicheAbundance()'s question, not a niche slope.
Read it for the data behind a slope, whether the log1p scale is roughly
linear, and whether a few extreme cells drive it. The display normalises
depth by library size (the model fits a slope on it). Every layer carries
patient, so + facet_wrap(~ patient) splits it by patient.
a ggplot; p$data has one row per (gene, patient, bin):
gene, patient, bin, expr (mean CP10k),
density (mean log1p density), condition and
gene_lab. Bins without counts are left out of that layer (not out
of the means).
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotNicheResponse(spe, res, gene = "G1", index = "A", niche = "B") plotNicheResponse(spe, res, gene = "G1", index = "A", niche = "B", adjust = FALSE)data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotNicheResponse(spe, res, gene = "G1", index = "A", niche = "B") plotNicheResponse(spe, res, gene = "G1", index = "A", niche = "B", adjust = FALSE)
Every index type's shared fit gives each patient an intercept per gene
(see patientIntercepts()): they absorb every difference between patients
– the condition's main effect, composition, batch – which is what makes
the niche slopes within-patient slopes. This relates the principal
components of the patients' gene-centred intercepts to patient covariates:
every patient-level colData column, plus the patients' index cells (log10),
mean log library size, and mean log1p density of each tested niche type.
Each tile is the adjusted of a component on a covariate.
plotPatientEffects( x, index, niche = NULL, covariates = NULL, type = c("association", "pca"), colour.by = NULL, n.pcs = 5L, ntop = 500L )plotPatientEffects( x, index, niche = NULL, covariates = NULL, type = c("association", "pca"), colour.by = NULL, n.pcs = 5L, ntop = 500L )
x |
a SpiDEResults or SpiDEFit (spiDE >= 0.99.34). |
index |
the index cell type. |
niche |
|
covariates |
|
type |
|
colour.by |
for |
n.pcs |
the number of components shown. |
ntop |
the number of genes, most expressed first, whose intercepts enter the components (a gene with any non-finite value is left out). |
Structure in the intercepts is harmless to the tests: it is absorbed by
design. With niche, the same analysis runs on the patients' slopes on
that niche (each gene's slopes scaled by its between-patient spread, as in
spiGSEA()): a slide or batch that organises the slopes can reach the
condition test, and belongs in strata.
a ggplot. For "association", p$data has pc,
covariate, r2 and group (patient or derived); for
"pca", patient, PC1, PC2 and the covariate.
patientIntercepts(), testNicheAbundance() for the composition
association gene by gene.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotPatientEffects(res, index = "A") plotPatientEffects(res, index = "A", niche = "B")data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotPatientEffects(res, index = "A") plotPatientEffects(res, index = "A", niche = "B")
The evidence the tests weigh: every patient's own slope of the gene's expression in the index type on the niche type's density (log1p), with an interval from its within-patient variance – the variance the tests use – and a point sized by the patient's weight in them. The pooled estimate and its confidence interval are the grey band. With a condition the patients are split by it, each group's weighted mean is a coloured line, and each strip gives the condition-specific slope difference and its q-value.
plotPatientSlopes(x, gene, index, niche, ci = 0.95)plotPatientSlopes(x, gene, index, niche, ci = 0.95)
x |
a SpiDEResults from the slopes engine. |
gene |
one or more genes, one column each. |
index, niche
|
the index and niche cell types. |
ci |
the confidence level of the intervals. |
Read it for whether a call rests on patients that agree or on one or two heavily weighted ones, and for why a condition test is missing (each group needs three effective patients). Patients are ordered by group, then by their mean slope over the genes shown, so a row is one patient across genes.
a ggplot. p$data has one row per (gene, patient) with a
slope: gene, patient, slope, se (within
patient), weight (the patient's share of the pooled test's
weights), condition and gene_lab (the strip label, plotmath).
Change it like any ggplot, e.g. + theme(axis.text.y = element_blank()).
patientSlopes() for the numbers, plotNicheResponse() for the cells.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotPatientSlopes(res, gene = c("G1", "G2"), index = "A", niche = "B")data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", verbose = FALSE) plotPatientSlopes(res, gene = c("G1", "G2"), index = "A", niche = "B")
Histograms of every triplet's p-value, one row per test and one column per fifth of each index type's genes by mean expression, with the level of a uniform distribution dashed. A spike at 0 on a flat floor is signal. A floor that rises toward 0 only in the brightest genes, or a hump at 1 (low power, conservative p-values), deserves a look at the calibration vignette. This is a sanity check, not a calibration test: real signal also moves p-values toward 0, and calibration needs data where there is nothing to find.
plotPValues(x, bands = 5L)plotPValues(x, bands = 5L)
x |
a SpiDEResults. |
bands |
the number of expression bands. |
a ggplot; p$data has one row per tested triplet with
p, test and band.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, procedure = "all", verbose = FALSE) plotPValues(res)data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, procedure = "all", verbose = FALSE) plotPValues(res)
Every tested triplet, placed by how much higher the gene is in the niche type than in the index type (log2, from each cell type's mean counts per 10,000) against the signed z of its test; calls are coloured by direction. A gene several-fold higher in the niche type whose measured level in the index cells rises with that niche's density is what transcripts leaking across segmentation boundaries would produce: the association is real in the data, but it needs orthogonal evidence (nuclear counts, a spillover-corrected segmentation) before it is read as biology.
plotSpillover( x, spe, test = "pooled", fold = 4, label = 6L, fdr = x@fdr, assay = "counts", name = "Niche" )plotSpillover( x, spe, test = "pooled", fold = 4, label = 6L, fdr = x@fdr, assay = "counts", name = "Niche" )
x |
a SpiDEResults. |
spe |
the SpatialExperiment the results came from. |
test |
|
fold |
the niche-marker threshold (shaded). |
label |
the number of markers labelled (strongest first). |
fdr |
the q-value threshold of a call. |
assay, name
|
the counts assay and the niche reducedDim prefix. |
a ggplot; p$data is the test's rows of the results table
with fold (log2), z and status.
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, verbose = FALSE) plotSpillover(res, spe)data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, verbose = FALSE) plotSpillover(res, spe)
A forest plot: each triplet's estimate (the pooled slope, or the slope
difference between conditions) with its 95\
where it is not called. gene, index and niche filter
the table before the n smallest p-values are taken, so
gene = c("CD74", "HLA-DRA") shows those genes' triplets. With
spe, a hollow triangle marks a gene at least fold-fold higher
in the niche type than in the index type (see plotSpillover()).
plotTopTriplets( x, test = NULL, n = 30L, gene = NULL, index = NULL, niche = NULL, spe = NULL, fold = 4, fdr = x@fdr, assay = "counts", name = "Niche" )plotTopTriplets( x, test = NULL, n = 30L, gene = NULL, index = NULL, niche = NULL, spe = NULL, fold = 4, fdr = x@fdr, assay = "counts", name = "Niche" )
x |
a SpiDEResults. |
test |
|
n |
the number of triplets. |
gene, index, niche
|
optional filters. |
spe |
|
fold |
the niche-marker threshold. |
fdr |
the q-value threshold of a call. |
assay, name
|
the counts assay and the niche reducedDim prefix. |
a ggplot; p$data is the plotted rows of the results table
with lo, hi, called, direction,
marker and label (plotmath).
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, verbose = FALSE) plotTopTriplets(res, test = "pooled", n = 15, spe = spe) plotTopTriplets(res, gene = "G1")data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) res <- spiDE(spe, condition = "condition", sigma = 30, verbose = FALSE) plotTopTriplets(res, test = "pooled", n = 15, spe = spe) plotTopTriplets(res, gene = "G1")
Features (genes, or gene sets from spiGSEA()) against every (index,
niche) column: the fill is the test's t (gold: the slope rises with the
niche's density, or is higher in the second condition; violet: the
reverse), a dot marks a call, grey a feature not tested in that index type.
Rows are clustered on their t profiles. Read it for a gene's or pathway's
response across all cell-type pairs, and for pairs that share responses.
plotTripletHeatmap( x, features = NULL, test = NULL, n = 30L, fdr = NULL, limit = NULL, cluster = TRUE )plotTripletHeatmap( x, features = NULL, test = NULL, n = 30L, fdr = NULL, limit = NULL, cluster = TRUE )
x |
a SpiDEResults (genes) or the table |
features |
|
test |
|
n |
the number of features when |
fdr |
the q-value threshold of a call (default: the results', or 0.05 for sets). |
limit |
the |t| at which the fill saturates (default: the 98th percentile). |
cluster |
order rows by hierarchical clustering of their t profiles. |
a ggplot; p$data has the plotted rows with feature
and called.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, verbose = FALSE) plotTripletHeatmap(res, test = "pooled", n = 10)data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, verbose = FALSE) plotTripletHeatmap(res, test = "pooled", n = 10)
Extract spiDE's test results
results(object, ...) ## S4 method for signature 'SpiDEResults' results(object, test = NULL, fdr = NULL, ...)results(object, ...) ## S4 method for signature 'SpiDEResults' results(object, test = NULL, fdr = NULL, ...)
object |
a SpiDEResults. |
... |
unused. |
test |
|
fdr |
|
a data.frame, one row per (gene, index, niche, test), ordered by
p-value: estimate (the pooled slope, or the difference in slope
between conditions, on the log scale per unit log1p niche density),
se, t, df, p, q (BH within the test's
family), in_family and, for the slopes engine,
p.heterogeneity (the label-free test that the triplet's slopes
vary between patients, the filter of procedure = "heterogeneity").
Condition-specific rows outside the family have q = NA.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") results(res, test = "pooled")data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") results(res, test = "pooled")
Builds the niche covariates at bandwidth sigma (unless spe
already carries them), fits the per-patient niche models (fitSpiDE()) and
tests them (testSpiDE()): the pooled niche test always, and the
condition-specific test when a condition is given.
spiDE(spe, condition = NULL, ...) ## S4 method for signature 'SpatialExperiment' spiDE( spe, condition = NULL, sigma, engine = c("slopes", "sandwich"), procedure = c("pooled_or_heterogeneity", "heterogeneity", "filtered", "all"), fdr = 0.05, cell_type = "cell_type", sample_id = "sample_id", name = "Niche", BPPARAM = BiocParallel::SerialParam(), verbose = TRUE, ... )spiDE(spe, condition = NULL, ...) ## S4 method for signature 'SpatialExperiment' spiDE( spe, condition = NULL, sigma, engine = c("slopes", "sandwich"), procedure = c("pooled_or_heterogeneity", "heterogeneity", "filtered", "all"), fdr = 0.05, cell_type = "cell_type", sample_id = "sample_id", name = "Niche", BPPARAM = BiocParallel::SerialParam(), verbose = TRUE, ... )
spe |
a SpatialExperiment with raw counts. |
condition |
|
... |
further arguments to |
sigma |
the niche bandwidth. |
engine |
|
procedure, fdr
|
passed to |
cell_type, sample_id
|
the colData columns of cell type and patient. |
name |
the niche reducedDim prefix: the niches are read from, or built
into, |
BPPARAM |
a BiocParallelParam. |
verbose |
report progress. |
a SpiDEResults.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") head(results(res))data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") head(results(res))
spiDE gives each role one colour and uses it in every plot: gold and violet for a slope that rises or falls with niche density (and for the diverging t scale), dark violet and pink for the first and second level of the condition, the viridis "magma" ramp for magnitudes (light to dark on the page), a brighter part of the same ramp (dark to bright) for plots in tissue coordinates, drawn on a dark panel, and the Okabe-Ito colours for a nuisance factor such as slide. The colours were checked under simulated deuteranopia, protanopia and tritanopia. Use them to make your own figures match, or pass others to a ggplot2 scale to override a spiDE plot's.
spiDEColours(role = NULL, n = 9L)spiDEColours(role = NULL, n = 9L)
role |
|
n |
the number of colours of a ramp, or of nuisance levels. |
a character vector of colours.
spiDEColours("condition") spiDEColours("magnitude", 5) spiDEColours("spatial", 5)spiDEColours("condition") spiDEColours("magnitude", 5) spiDEColours("spatial", 5)
Returned by fitSpiDE(); tested with testSpiDE().
## S4 method for signature 'SpiDEFit' x$name## S4 method for signature 'SpiDEFit' x$name
x |
a SpiDEFit or SpiDEResults. |
name |
a slot name. |
For the slopes engine each element of index holds
genes, niches (tested niche columns), patients,
beta (genes x patients x niches: each patient's slope),
v_tile (the within-patient spatial sandwich variance of each slope),
v_model (the model variance), factor (patients x niches: the
per-patient variance factor), pooled (genes x niches: the shared
fit's slope), mean_expr and ncells. For the sandwich engine
it holds genes, niches and a data.frame coef of the
tested coefficients (test = "pooled" or "condition") with their
CR2 standard errors and Bell-McCaffrey degrees of freedom. Both hold
status (genes fitted; genes left out because the solver could not
fit them or their coefficients ran off, the latter also counted as
runaway; and genes whose dispersion is on its search bound) and, under
depth = "spatial_spline", depth_r2 (patients x niches: the
of each niche column on the patient's library-size spline, a
diagnostic of how much of the niche covariate the spline could absorb).
From spiDE 0.99.34 both also hold intercept (genes x patients:
each patient's intercept in the shared, condition-free fit; see
patientIntercepts()), niche_mean (patients x niches: the mean
log1p density of each tested niche over the patient's index cells),
loglib_mean (the patients' mean log library size) and
ncells (the patients' index cells).
An object of class SpiDEFit, created by fitSpiDE(); its
$ accessor returns a slot, and show() prints a summary.
enginea character, "slopes" or "sandwich".
sigmaa numeric, the niche bandwidth the fit used.
conditiona character (length 0 or 1), the condition the fit was given (the sandwich engine fits condition-specific slopes only if it was given one; the slopes engine only records it).
patientsa data.frame, one row per patient (sample): its id, the
number of cells it contributes, and every colData column that is constant
within patients (so a condition can be named at testSpiDE() time).
indexa named list, one element per index cell type, holding the engine's estimates (see Details).
paramsa list of the settings the fit used.
Returned by testSpiDE() and spiDE(); read it with results().
## S4 method for signature 'SpiDEResults' x$name## S4 method for signature 'SpiDEResults' x$name
x |
a SpiDEResults. |
name |
a slot name. |
An object of class SpiDEResults, created by testSpiDE();
its $ accessor returns a slot, and show() prints a summary.
tablea data.frame, one row per (gene, index, niche, test).
conditiona character (length 0 or 1), the condition tested.
contrasta character, e.g. "Responder - Non-responder".
procedurea character, the condition test's family:
"pooled_or_heterogeneity", "heterogeneity",
"filtered" or "all".
fdra numeric, the FDR level of the filter.
fitthe SpiDEFit tested.
Asks, for every gene set and every (index cell type, niche cell type) column, whether the set's genes as a group change with the niche (the pooled test) or change with it differently between the two conditions (the condition-specific test).
spiGSEA(object, genesets, ...) ## S4 method for signature 'SpiDEResults' spiGSEA( object, genesets, test = NULL, type = c("competitive", "self-contained"), min.size = 5L, max.size = 500L, BPPARAM = BiocParallel::SerialParam(), ... )spiGSEA(object, genesets, ...) ## S4 method for signature 'SpiDEResults' spiGSEA( object, genesets, test = NULL, type = c("competitive", "self-contained"), min.size = 5L, max.size = 500L, BPPARAM = BiocParallel::SerialParam(), ... )
object |
a SpiDEResults from the slopes engine. |
genesets |
a named list of character vectors of gene names. |
... |
unused. |
test |
|
type |
|
min.size, max.size
|
set size limits after intersecting with the genes tested in each index type. |
BPPARAM |
a BiocParallelParam; the condition-specific test is run in parallel over sets. |
A gene set is tested exactly like a gene. Every patient gets a slope for the
set: the average of its genes' niche slopes in that patient, each gene's
slope first put on a common scale (divided by the gene's typical spread of
slopes between patients, so a highly expressed gene does not dominate). By
default the set's slope is competitive: the same average over the
column's other tested genes is subtracted, so the test asks whether the set
moves more than the typical gene. This guards against a shift that
moves every gene, e.g. from segmentation spillover. With type =
"self-contained" nothing is subtracted.
The patients' set slopes then go through the same pooled and condition
tests as the genes' slopes (see testSpiDE() and the model vignette).
Because the patient is the unit of replication, the correlation between the
genes of a set is accounted for without being estimated: it is already in
how much the set's slope varies from patient to patient.
Needs the slopes engine, whose fit holds each patient's slopes. For gene sets
in the between-patient test, use testNicheAbundance() with genesets.
a data.frame, one row per (set, index, niche, test):
size (genes of the set tested there), estimate (the mean set
slope, pooled, or its difference between conditions, in units of the
genes' between-patient spread), se, t, df, p,
q (Benjamini-Hochberg within the (index, niche, test) column),
q.global (over every row of that test), n_patients and
direction.
data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") sets <- list(first = paste0("G", 1:5), second = paste0("G", 6:12)) spiGSEA(res, sets, min.size = 3)data(toySpiDE) res <- spiDE(toySpiDE, condition = "condition", sigma = 30, index = "A", procedure = "all") sets <- list(first = paste0("G", 1:5), second = paste0("G", 6:12)) spiGSEA(res, sets, min.size = 3)
A between-patient question that the within-patient engines cannot see:
does a gene's expression in index cell type differ between patients
whose type- cells sit, on average, among more type- cells?
(fitSpiDE() and testSpiDE() instead ask whether expression changes with
niche density WITHIN each patient.)
testNicheAbundance(spe, ...) ## S4 method for signature 'ANY' testNicheAbundance( spe, condition = NULL, sigma, index = NULL, niche = NULL, covariates = character(), assay = "counts", cell_type = "cell_type", sample_id = "sample_id", name = "Niche", min.cells = 10L, prior.count = 1, max.leverage = 3, genesets = NULL, type = c("competitive", "self-contained"), min.size = 5L, max.size = 500L, verbose = TRUE )testNicheAbundance(spe, ...) ## S4 method for signature 'ANY' testNicheAbundance( spe, condition = NULL, sigma, index = NULL, niche = NULL, covariates = character(), assay = "counts", cell_type = "cell_type", sample_id = "sample_id", name = "Niche", min.cells = 10L, prior.count = 1, max.leverage = 3, genesets = NULL, type = c("competitive", "self-contained"), min.size = 5L, max.size = 500L, verbose = TRUE )
spe |
a SpatialExperiment with niche reducedDims (see |
... |
further arguments passed to the method. |
condition |
a character, the colData column of the condition (constant
within sample), or |
sigma |
a numeric, the bandwidth (one value). |
index, niche
|
character vectors restricting the index / niche cell types (NULL = all). An index type is never tested against its own niche. |
covariates |
a character vector of sample-level colData columns to adjust for (constant within sample; the first non-missing value per sample is used). |
assay |
a character, the counts assay. |
cell_type, sample_id
|
the colData columns of cell type and sample. |
name |
the niche reducedDim prefix. |
min.cells |
minimum cells of the index type a sample must contribute. |
prior.count |
the pseudocount in the log2-CPM. |
max.leverage |
a number above 1: samples whose leverage in a design
exceeds this multiple of the average leverage are down-weighted until none
does ( |
genesets |
a named list of character vectors (gene identifiers as in
|
type |
with |
min.size, max.size
|
with |
verbose |
report progress. |
For each index cell type and niche cell type , regresses the
per-sample pseudobulk log2-CPM of every gene in the type- cells on
the per-sample mean log1p niche density of type around those cells,
across samples, with a limma moderated -test. Without a
condition the design is ~ niche + covariates and one term is
reported, "niche": the association across all samples. With a
condition the design is ~ niche * condition + covariates
with the condition coded , and two terms are reported:
"niche", the association averaged over the two conditions, and
"condition:niche", its difference between them (second level minus
first) – the patient-level counterpart of the condition-specific niche test
that testSpiDE() runs within patients. Where the conditions leave fewer
than two samples each, only "niche" is reported, across all samples.
Influential samples. With tens of samples, one sample at the edge
of a niche type's abundances can carry the whole slope, and a gene whose
expression is extreme in that same sample then gets a far smaller
-value than the moderated allows. Samples whose leverage
(hat value) in the design exceeds max.leverage times the average are
therefore down-weighted until none does. The weights come from the design
alone, never from expression, so no gene's own values decide how much a
sample counts and the estimates stay unbiased; the permutation nulls set
the default.
Each row reports the largest leverage before the cap and how many samples
were down-weighted.
Gene sets. With genesets, sets are tested instead of genes,
each like a gene. Every gene of the index type is standardised across
samples, and a set's score in a sample is the mean of its genes' standardised
pseudobulk expression there. With type = "competitive" (the default),
the mean over all other genes is subtracted, so a shift that moves every gene
together does not make a set call. The scores then go through the same
regression, leverage cap and multiple testing as genes. This is the
between-patient counterpart of spiGSEA(), which tests sets on the
within-patient slopes. The set version has not yet been checked against
permutation nulls.
Multiple testing. q.global is a Benjamini-Hochberg
adjustment over every row of the same term, so the two terms are separate
families: a strong "niche" signal no longer loosens the threshold
for "condition:niche".
This is the association that fitSpiDE()'s per-patient intercepts
deliberately absorb. It is a between-patient effect with experimental units; it is not
neighbourhood-dependent differential expression, and the two must not be
conflated (see the model vignette).
It was compositionTest() up to spiDE 0.99.22.
Samples contributing fewer than
min.cells cells of the index type are dropped for that index type,
and an index type with fewer than three remaining samples is skipped.
a data.frame with one row per (gene, index type, niche type, term):
gene, index, niche, term, estimate
(log2-CPM per unit log1p density), t, p,
n_patients, leverage (the largest sample leverage before the
cap, as a multiple of the average), downweighted (the samples
down-weighted), q (BH within each (index, niche, term) over genes)
and q.global (BH over every row of the same term). With
genesets, one row per (set, index type, niche type, term): the
column gene becomes set, followed by size, and a
direction column is added; estimate is in standard
deviations of the set score per unit log1p density.
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) na <- testNicheAbundance(spe, condition = "condition", sigma = 30) head(na[order(na$p), ]) sets <- list(first = paste0("G", 1:6), second = paste0("G", 7:12)) nas <- testNicheAbundance(spe, condition = "condition", sigma = 30, genesets = sets, verbose = FALSE) head(nas)data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) na <- testNicheAbundance(spe, condition = "condition", sigma = 30) head(na[order(na$p), ]) sets <- list(first = paste0("G", 1:6), second = paste0("G", 7:12)) nas <- testNicheAbundance(spe, condition = "condition", sigma = 30, genesets = sets, verbose = FALSE) head(nas)
Always runs the pooled test: for every (gene, index cell type, niche
cell type), whether the gene's expression in the index type changes with the
density of the niche type, within patients, consistently across patients.
If a condition is given (here, or recorded by fitSpiDE()), it also
runs the condition-specific test: whether that slope differs
between the two conditions – a between-patient contrast of within-patient
slopes.
testSpiDE(object, condition = NULL, ...) ## S4 method for signature 'SpiDEFit' testSpiDE( object, condition = NULL, procedure = c("pooled_or_heterogeneity", "heterogeneity", "filtered", "all"), fdr = 0.05, strata = NULL, BPPARAM = BiocParallel::SerialParam(), ... )testSpiDE(object, condition = NULL, ...) ## S4 method for signature 'SpiDEFit' testSpiDE( object, condition = NULL, procedure = c("pooled_or_heterogeneity", "heterogeneity", "filtered", "all"), fdr = 0.05, strata = NULL, BPPARAM = BiocParallel::SerialParam(), ... )
object |
a SpiDEFit from |
condition |
|
... |
unused. |
procedure |
|
fdr |
the FDR level of the filter. |
strata |
|
BPPARAM |
a BiocParallelParam; the slopes engine's condition-specific test is run in parallel over genes. |
The condition-specific test is carried out within a family of triplets chosen without looking at the condition, so the choice is exact under any relabelling of patients and the multiplicity budget is spent on triplets with something for the condition to explain:
procedure = "pooled_or_heterogeneity" (the slopes engine's
default): the triplets that pass either filter below.
procedure = "heterogeneity": the triplets whose slopes vary
between patients more than their sampling variances allow (Cochran's Q
across patients, the statistic behind the between-patient variance;
Benjamini-Hochberg over every triplet at fdr). A condition effect
is such variation, including one whose two conditions' slopes cancel in
the pooled slope. With few patients the pooled test passes little, and
this filter carries the condition test.
procedure = "filtered" (the sandwich engine's default, and
the slopes engine's to spiDE 0.99.36): the triplets whose pooled test
passes at fdr. It keeps a modest condition difference on a strong
pooled slope, which the heterogeneity test, spread over every patient,
can miss.
procedure = "all": every triplet.
The heterogeneity filter needs each patient's slopes, so it is the slopes engine's only. The calibration vignette has the measurements behind the default.
The slopes engine's pooled test refers its moderated t statistic to limma's degrees of freedom scaled by the Kish effective share of the patients (the proportional rule), so a column whose weights rest on a few patients is not credited with the others. To spiDE 0.99.36 the df were the smaller of limma's df and the Kish effective number of patients less one, which is conservative with few patients.
a SpiDEResults; read it with results().
data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) res <- testSpiDE(fit, condition = "condition", procedure = "all") res head(results(res))data(toySpiDE) spe <- buildNiches(toySpiDE, sigma = 30) fit <- fitSpiDE(spe, index = "A", sigma = 30) res <- testSpiDE(fit, condition = "condition", procedure = "all") res head(results(res))
The theme of every spiDE plot: a black panel border, no grid lines, black
axis ticks, boxed facet strips and italic legend titles, with text scaled
by rl. Add it to your own ggplots so that they match spiDE's, and
change any element with a further + theme(...).
theme_spiDE(rl = 1.1)theme_spiDE(rl = 1.1)
rl |
a positive number, the text size relative to ggplot2's default. |
a ggplot2 theme.
library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_spiDE()library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_spiDE()
A seeded toy dataset of 20 genes across 2,400 cells from 16 patients (eight
Responders, eight Non-responders; one sample of 150 cells each) with three
cell types (A, B, C) laid out on a 500 x 500 field. Cell type B is concentrated in the right of the field, and
a neighbourhood-dependent signal is planted: gene G1 is up-regulated in the
index cell type A, in Responders, in proportion to the local density of the
niche cell type B. Used throughout the examples and vignette.
data(toySpiDE)data(toySpiDE)
A SpatialExperiment
with a counts assay and colData
columns sample_id, condition, cell_type, Age, Area, and nCount.
Generated by data-raw/make_toySpiDE.R.
data(toySpiDE) toySpiDEdata(toySpiDE) toySpiDE