Concepts¶
Steps and pipelines¶
A step smooths an (n_obs, n_genes) expression matrix over one embedding stored in
adata.obsm. A pipeline is an ordered list of steps: each consumes the previous step’s
output. [KompotGP(), KnnGaussian()] therefore denoises along the expression manifold first,
then smooths the already-denoised expression over physical space.
step |
default basis |
engine |
|---|---|---|
|
Gaussian kernel over |
|
|
FFT Nadaraya-Watson on a fine grid (KDEpy) |
|
|
Gaussian-process regression (kompot/mellon) |
Steps are frozen dataclasses – specifications, not fitted objects. They carry no data, so the
same pipeline can be reused across datasets and recorded verbatim in adata.uns.
Shorthands cover the common pipelines, so the common cases need no imports:
shorthand |
pipeline |
meaning |
|---|---|---|
|
|
spatial smoothing only |
|
|
cell-state smoothing only |
|
|
both, cell state first |
|
|
both, spatial first |
|
|
spatial, KDE engine |
|
|
spatial, GP engine |
|
|
no smoothing; raw score only |
Doing just one of the two is the ordinary case, not a special one: steps="spatial" and
steps="dm" are single-element pipelines.
Blending: a symmetric alternative to composition¶
"dm+spatial" composes the two smoothers: the spatial step consumes the cell-state-denoised
expression, so the result inherits the spatial parent’s footprint and collapses toward it.
steps="blend" (a Blend) does something different – it runs the
two views independently on the raw expression, standardises each branch’s score to a common
footing, averages them, and finally rescales the average back onto the raw score’s scale. Because
neither branch consumes the other, the blend stays a genuinely distinct field, roughly equidistant
from both parents rather than collapsing onto either:
ss.smooth(adata, genes, "sig", steps="blend") # default: spatial + dm
ss.smooth(adata, genes, "sig", steps=ss.Blend(calibrate="iqr")) # robust calibration
The range calibration is what keeps a blended field usable alongside the others. The raw
standardised average lives in z-units, which do not share a colour bar with raw, spatial
or dm. The final affine step matches the blend’s centre and spread to the raw score’s – by
mean and standard deviation (calibrate="std", the default) or by median and inter-quartile
range (calibrate="iqr", robust to outlier cells). The map is monotone, so it never reorders
cells; the realised scale and shift are recorded in
provenance().
Which composition to reach for. "spatial" is the fast default and the right first thing to
try – a truncated Gaussian over tissue coordinates, about a second on a full slide. Use "dm"
when position is irrelevant and you only want to denoise along the expression manifold. Between the
two combinations, prefer "blend" when you want space and cell state to contribute on equal
footing – cell-state denoising and spatial coherence, with neither view dominating; reach for
"dm+spatial" when you specifically want the composed field that leans on the spatial footprint.
The tutorial features "blend" for exactly this reason.
The storage contract¶
Everything smooth() produces is written into the AnnData:
key |
contents |
|---|---|
|
smoothed signature score, |
|
unsmoothed score from the same genes and combiner |
|
|
|
provenance: genes, pipeline, resolved bandwidths, package version |
Nothing else is needed to render the result later. spatial_smooth.plot reads those keys
and never recomputes, so a smoothed object can be written to .h5ad, shipped, reloaded, and
plotted in an environment without kompot, KDEpy or palantir installed. The package’s
test suite asserts this adversarially: it blocks those imports and replaces every compute entry
point with a function that raises, then renders a reloaded file.
provenance() reads the record back with the pipeline decoded, including
the bandwidth each step actually resolved to:
>>> ss.provenance(adata, "sig")["steps"]
[{'kind': 'knn_gaussian', 'basis': 'spatial', 'k': 400, 'sigma': None,
'sigma_factor': 6.0, 'workers': -1,
'resolved': {'sigma_nominal': 78.3, 'sigma_used': 78.3, 'k_used': 400,
'kernel_mass_retained': 0.955,
'sigma_effective': 71.6,
'sigma_effective_p1': 52.2, 'sigma_effective_p99': 84.1}}]
sigma_nominal is the Gaussian’s width before truncation – a bandwidth no cell actually
experiences. (sigma_used is the same number under its older name.) sigma_effective is what
the k-truncated kernel behaves like. Quote the latter – see Truncation: what bandwidth the data actually sees below.
Smooth every gene once, then derive any signature for free¶
The Gaussian process is the one expensive step – a diffusion-map GP over a whole panel takes
minutes, where the kNN kernel takes about a second. So do it once.
smooth_all() smooths every gene through a view and stores the result; every
later smooth() call with all_genes=True – for a signature, a single
gene, or a "blend" – reads those pre-smoothed layers instead of recomputing:
ss.compute_diffusion_map(adata) # Palantir -> obsm["DM_EigenVectors"]
ss.smooth_all(adata, steps="spatial") # every gene, once, over space
ss.smooth_all(adata, steps="dm") # every gene, once, over the diffusion-map GP
# all of these are now cache hits -- no smoother runs again:
ss.smooth(adata, genes, "spatial_only", steps="spatial", all_genes=True)
ss.smooth(adata, genes, "dm_only", steps="dm", all_genes=True)
ss.smooth(adata, genes, "composed", steps="dm+spatial", all_genes=True)
ss.smooth(adata, genes, "blended", steps="blend", all_genes=True)
The smoothers are per-gene operations, so a signature’s smoothed columns gathered from the
all-genes layer are identical to smoothing that signature alone – bit-for-bit for the neighbour
smoothers, to floating-point precision for the GP. A four-mode comparison therefore pays for the
diffusion-map GP exactly once; "composed" and "blended" both re-use that single solve.
Underneath, each smoother’s output is memoized on the AnnData, keyed by a stable hash of its
input matrix, parameters and basis (spatial_smooth.CACHE_KEY). The cache is bounded
(cache_max_entries=N) and opt-out (cache=False);
clear_smooth_cache() drops the cache artifacts while leaving the results
intact.
Scoring, and why gene-level smoothing is free¶
The multi-gene score is mean_z by default: standardise each gene, then average. The mean and
standard deviation always come from the raw matrix, for both the raw and the smoothed score.
Two consequences, both intended.
First, raw and smoothed scores share one scale, so they belong on a common colour bar.
Second, for a row-stochastic smoother – one whose weights sum to one, so it maps a constant field to itself – smoothing the genes and then scoring is exactly scoring and then smoothing the score. Writing \(W\) for the operator and \(x_j\) for gene \(j\),
because \(W\mathbf{1} = \mathbf{1}\) carries the constants through untouched.
KnnGaussian and Kde are both
row-stochastic. The pipeline therefore smooths per gene at no cost in correctness, which is what
keeps a Gaussian-process step – linear, but not row-stochastic, so it does not commute –
meaningful inside the same framework.
Choosing a smoother¶
step |
engine |
full slide (~1.6e5 cells) |
gives you |
|---|---|---|---|
|
truncated Gaussian over |
~1 s |
the default; fast, sharp |
|
FFT Nadaraya-Watson on a fine grid |
~1 s |
a rendered field; resolution-bound |
|
Gaussian-process regression |
minutes |
a length scale, a posterior, fit-on-one-condition |
KnnGaussian and KompotGP produce visually equivalent spatial fields, with the GP
marginally sharper when given enough landmarks; the kNN kernel is roughly two orders of magnitude
faster per gene. Reach for the GP when you want its extras – an explicit length scale,
uncertainty, or groupby/condition (fit on one condition, evaluate everywhere) – and for
smoothing over a diffusion map, where it is the established choice.
Truncation: what bandwidth the data actually sees¶
KnnGaussian evaluates its Gaussian only over each cell’s k
nearest neighbours. That truncation is not free, and it is the number most likely to end up
misquoted in a methods section, so the package measures it.
Write \(d_{ij}\) for the distance from cell \(i\) to its neighbour \(j\), and
\(\mathcal{N}_k(i)\) for the k nearest neighbours of \(i\) (itself included). The
weights are
so the operator is row-stochastic by construction. Let \(r_k(i)\) be the distance to the
k-th neighbour. Two diagnostics follow, both recorded by
provenance().
Retained kernel mass. For a continuous isotropic 2-D Gaussian the mass inside radius
\(r\) is \(1 - e^{-r^{2}/2\sigma^{2}}\), so an estimate of the fraction of the untruncated
kernel that survives the k-neighbour cutoff is
Effective bandwidth. A 2-D Gaussian satisfies \(\mathbb{E}[d^{2}] = 2\sigma^{2}\). Inverting that on the realised weights gives the bandwidth the truncated kernel behaves like:
This \(m\) is an estimate, not an exact accounting: it equals the discrete retained weight
fraction \(\sum_{j \in \mathcal{N}_k(i)} \tilde{w}_{ij} \big/ \sum_{j} \tilde{w}_{ij}\)
(unnormalised \(\tilde{w}\)) only where the local point density is uniform. Measured against
the discrete truth on a real section it is biased low by 0.006 at k=400 and 0.022 at
k=100 – i.e. it under-reports retention, so the warning below fires slightly early rather
than slightly late.
Because \(r_k(i)\) is fixed by a neighbour count, it contracts where cells are dense and
expands where they are sparse: \(\sigma_{\mathrm{eff}}\) varies from cell to cell, and the
smoother is truncated-Gaussian and implicitly density-adaptive rather than fixed-bandwidth.
provenance() therefore stores kernel_mass_retained (\(m\)), sigma_effective
(\(\overline{\sigma_{\mathrm{eff}}}\)) and its 1st/99th percentiles, and a UserWarning
fires when \(m < 0.9\).
Measured on the tutorial’s 10x Xenium mouse-brain section (\(n = 36{,}602\), median nearest-neighbour distance 13.05 µm, hence nominal \(\sigma = 6 \times 13.05 = 78.3\) µm):
quantity |
|
|
|---|---|---|
retained mass \(m\) |
0.58 |
0.96 |
\(\overline{\sigma_{\mathrm{eff}}}\) |
47.3 µm |
71.6 µm |
\(\sigma_{\mathrm{eff}}\) p1 → p99 |
24.9 → 68.6 µm |
52.2 → 84.1 µm |
At k = 100 the nominal \(\sigma\) overstates the realised bandwidth by roughly 40%, which
is why the default is k = 400. Quote sigma_effective, never sigma_nominal.
Implementation: spatial_smooth.smoothers.knn_gaussian_operator() (return_info=True).
Bandwidths are scale-invariant¶
Every default bandwidth is a multiple of the median nearest-neighbour distance of the coordinates,
so the same factor smooths the same amount whether the coordinates are in microns or millimetres.
KnnGaussian(sigma_factor=6.0) is about six cell spacings (~50 um on a section with 8 um
spacing).
KompotGP inherits mellon’s empirical length scale,
\(\ell = \texttt{ls\_factor} \cdot \operatorname{geomean}(d_\mathrm{NN}) \cdot e^{3}\)
(mellon.parameters.compute_ls), which is scale-invariant for the same reason. Passing an
explicit ls (or an explicit sigma to KnnGaussian) bypasses the empirical path and puts
you back in absolute coordinate units.
Warning
Over a diffusion map, kompot’s native ls_factor=10 is right. Over physical coordinates
it is roughly 200x the cell spacing and collapses the field into a single global gradient. Use
ls_factor near 0.3 there – which is precisely what the "spatial-gp" shorthand does.
A small effective length scale also needs enough landmarks to resolve it: keep n_landmarks
large enough that the landmark spacing stays below the length scale.
Plot control: kwargs pass straight through¶
spatial_smooth.plot.signature() is a thin wrapper, not a reimplementation. Everything after
name is forwarded verbatim to the backend, and color is set for you from the stored
provenance; package defaults (cmap, percentile colour limits, a grey na_color) are injected
only for keys you did not pass.
|
underlying call |
|---|---|
|
|
|
|
|
|
|
squidpy if installed, else scanpy |
ss.pl.signature(adata, "sig", backend="scanpy", cmap="magma", vmax="p99.5", frameon=False)
ss.pl.signature(adata, "sig", backend="squidpy", figsize=(6, 6))
spatial_smooth.plot.compare() lays several stored results side by side (optionally with the
raw score as the first panel), which is how the mode figures in this documentation are built.
Restricting, conditioning, and other engines¶
A few further knobs, each one argument away:
Subset of cells.
ss.smooth(adata, genes, name, subset_key="annot", include=["CA1"])smooths only the named cells and returns a new, smallerAnnData– cells filtered out neither train the smoother nor receive the field, so use the return value.Fit on one condition, evaluate everywhere.
KompotGP(groupby=..., condition=...)trains the GP on one group and imputes the field for all cells – useful when one arm of an experiment is the reference.Other engines.
Kde(a fine-grid FFT Nadaraya-Watson estimator via KDEpy) renders a field rather than a neighbour average; the"spatial-gp"shorthand puts the Gaussian process on tissue coordinates with a sensiblels_factor:ss.smooth(adata, genes, "kde", steps="spatial-kde") ss.smooth(adata, genes, "spatial_gp", steps="spatial-gp")