-
Notifications
You must be signed in to change notification settings - Fork 0
Analysis Templates
Every analysis the viewer runs is a template: a piece of plain scverse source
with a declared contract, kept in src/palms/utils/step_templates/builtin/
as a .tmpl file. When you click a button, the owning tab builds a dictionary of
parameters, the template is rendered once, and that same string is both executed
and recorded in the provenance graph — so the code below is the code that runs,
and the code the exported notebook replays.
This page is the catalogue: what each template computes, its contract, and its full default source. It is generated from the registry itself, so it cannot disagree with what is installed.
- To change any of this, see the Templates tab — overrides are resolved per block, so you can edit one step without freezing the rest.
- To see what a specific run did, see the Notebook tab; the provenance graph records the rendered code, not the template.
- The API behind all of it is in the API Reference.
A template is an ordered set of named blocks, delimited by #--- block <name>
lines. Everything structural is a comment, so the file is valid Python.
| Contract field | Meaning |
|---|---|
params |
The $name placeholders the template expects. A ? marks an optional one. Values are substituted as Python literals. |
requires |
Names the template may use from the executor namespace — sc, sq, sd, pd, np, plt, Path, data_path, sdata, adata, plus anything a declared dependency binds. Referencing anything else is rejected before it can fail as a NameError on replay. |
outputs |
Results the template must bind. The executor raises if one is missing, so a template edit that stops producing a result fails immediately instead of returning stale state. |
assemblies |
The block combinations that are allowed to run. The call site decides which one to use, because the branch structure is what the widgets mean; the registry owns the text. |
| frozen blocks | Blocks that may not be customised, listed where they occur. |
Which assembly runs is decided by the tab, from your widget settings. The Templates tab's preview pane shows the exact assembled and substituted string for the settings currently in the owning tab.
-
normalize— Normalised copy every expression-based analysis reads. -
spatial_neighbors— Spatial neighbour graph, built on the normalised copy. -
clustering.leiden— Leiden community detection on a copy of the normalised data. -
genes.cnv_infercnv— inferCNV over a reference population, run in-process. -
genes.correlation— Two-gene scatter with Pearson/Spearman, one expr.* block per normalisation. -
genes.marker_plot— Scanpy marker-gene figure; one call.* block per plot type. -
genes.rank_genes— Rank marker genes per cluster, keyed by clustering. -
umap.plot— UMAP scatter coloured by gene expression or by a clustering. -
spatial.cooccur— Squidpy co-occurrence across increasing radii. -
spatial.ligrec— Ligand-receptor permutation test (omnipath resources). -
spatial.nhood— Squidpy neighbourhood enrichment z-scores. -
roi.deg— Differential expression between drawn ROI regions. -
roi.export_expression— Write the per-cell ROI expression table to CSV. -
roi.expression— Per-region expression of one gene, with pairwise Welch tests. -
roi.polygons— The ROI polygons drawn in the viewer, inlined as literals.
Run by: every expression-based tab
Builds adata_norm, the normalised copy that every expression-based analysis reads: counts per cell scaled to 10,000, log1p, then PCA. Nothing else normalises, and nothing normalises adata in place.
The copy is the whole point. An earlier version normalised adata itself, which made the provenance graph lie: a step that read adata got different numbers depending on whether some other tab had run first. Every downstream template declares deps=["normalize"] and reads adata_norm, so the dependency is explicit and the notebook reproduces it in the right order. ctx.ensure_normalized() is idempotent — running it twice does not re-normalise.
Contract
No parameters.
-
Requires:
adata,sc -
Outputs:
adata_norm -
Blocks:
main
Default source
# Normalized copy used by expression-based analyses
adata_norm = adata.copy()
sc.pp.normalize_total(adata_norm, target_sum=1e4)
sc.pp.log1p(adata_norm)
sc.pp.pca(adata_norm)Run by: Neighborhood Enrichment and Ligand-Receptor
Builds the spatial neighbour graph squidpy's spatial statistics need, on adata_norm, from each cell's k nearest neighbours. Uses spatial_neighbors_knn, which replaced spatial_neighbors(coord_type='generic') — removed in squidpy 1.9. The two produce an identical graph, so results recorded before the change still match.
This is the one template with no live preview in the Templates tab. k comes from whichever tab called ctx.ensure_spatial_neighbors(), and two tabs have their own slider, so a preview provider would have to pick one of them arbitrarily. It declares a sample-params value instead and the preview says so.
Contract
| Parameter | Type | Required |
|---|---|---|
n_neighs |
int |
yes |
-
Requires:
adata_norm,sq - Outputs: nothing (terminal step)
-
Blocks:
main
Default source
# Spatial neighbor graph (k=$n_neighs)
sq.gr.spatial_neighbors_knn(adata_norm, n_neighs=$n_neighs)Run by: Clustering
Leiden community detection on a copy of adata_norm, writing the labels back to adata.obs[key]. The result is what every cluster-keyed analysis — ranked genes, neighbourhood enrichment, co-occurrence, ligand-receptor — is grouped by.
leiden_labels is a declared output rather than something read back off ctx.adata afterwards. That matters: reading back worked only while the executor namespace and ctx.adata were the same object, and the executor now raises if a template stops binding a declared output, so an edit that breaks the result fails loudly.
Contract
| Parameter | Type | Required |
|---|---|---|
key |
str |
yes |
resolution |
float |
yes |
n_neighbors |
int |
yes |
n_pcs |
int |
yes |
n_top_genes |
int |
no |
flavor |
str |
yes |
n_iterations |
int |
yes |
directed |
bool |
yes |
random_state |
int |
yes |
-
Requires:
adata,adata_norm,sc -
Outputs:
leiden_labels -
Blocks:
head,hvg,scale,pca,tail
Variants — 4 assemblies:
The four assemblies are the four combinations of the two preprocessing checkboxes. Highly-variable-gene selection adds hvg, scaling adds scale, and either one adds pca — because PCA has to be recomputed once the matrix has changed. With neither, the PCA already computed by normalize is reused and only head+tail runs.
head + tailhead + scale + pca + tailhead + hvg + pca + tailhead + hvg + scale + pca + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# Leiden clustering ($key)
adata_leiden = adata_norm.copy()
#--- block hvg
sc.pp.highly_variable_genes(adata_leiden, n_top_genes=$n_top_genes, flavor='seurat')
adata_leiden = adata_leiden[:, adata_leiden.var.highly_variable].copy()
#--- block scale
sc.pp.scale(adata_leiden, max_value=10)
#--- block pca
sc.pp.pca(adata_leiden)
#--- block tail
sc.pp.neighbors(adata_leiden, n_neighbors=$n_neighbors, n_pcs=$n_pcs)
sc.tl.leiden(
adata_leiden, resolution=$resolution, key_added=$key,
flavor=$flavor, n_iterations=$n_iterations, directed=$directed,
random_state=$random_state,
)
adata.obs[$key] = adata_leiden.obs[$key].values
leiden_labels = adata.obs[$key]Run by: CNV
Infers copy-number variation against a reference population with inferCNV, in-process: prepares the input, builds CNV neighbours, clusters at several resolutions and scores each cell.
arrow_shim is a frozen block — it cannot be customised even with overrides enabled. It is a version workaround whose byte identity with a helper elsewhere is pinned by a test, and no validation gate could have told you that editing it breaks on the next pandas. The Templates tab lists frozen blocks in the Contract pane for exactly this reason.
Contract
| Parameter | Type | Required |
|---|---|---|
reference_clustering |
str |
yes |
reference_obs_key |
str |
yes |
reference_categories |
list |
yes |
n_neighbors |
int |
yes |
smoothing_neighbors |
int |
yes |
window_size |
int |
yes |
step |
int |
yes |
lfc_clip |
float |
yes |
resolution |
float |
yes |
min_mapped_fraction |
float |
yes |
include |
list |
no |
-
Requires:
adata,np,pd,sc -
Outputs:
adata_cnv,cnv_cluster_keys,cnv_clusters,cnv_score -
Blocks:
head,subset,prepare,arrow_shim,tail -
Frozen (not customisable):
arrow_shim
Variants — 2 assemblies:
Two assemblies, differing only in whether a subset of clusters was selected to run on.
head + prepare + arrow_shim + tailhead + subset + prepare + arrow_shim + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# CNV inference (inferCNV): reference $reference_clustering
from insitucnv.tl import (
prepare_cnv_input, compute_cnv_neighbors, cluster_cnv_resolutions, run_infercnv,
)
adata_cnv = adata.copy()
#--- block subset
# limit the analysis to the selected cell types plus the reference population
adata_cnv = adata_cnv[
adata_cnv.obs[$reference_clustering].astype(str).isin($include)
].copy()
#--- block prepare
adata_cnv.obs[$reference_obs_key] = adata_cnv.obs[$reference_clustering].values
_n_panel = adata_cnv.n_vars
_raw_counts = adata_cnv.X.copy()
sc.pp.normalize_total(adata_cnv, target_sum=1e4)
sc.pp.log1p(adata_cnv)
sc.pp.pca(adata_cnv)
sc.pp.neighbors(adata_cnv, n_neighbors=$n_neighbors)
adata_cnv.layers['raw_counts'] = _raw_counts
adata_cnv = prepare_cnv_input(
adata_cnv, raw_layer='raw_counts', smoothing_neighbors=$smoothing_neighbors,
add_gene_positions=True, drop_unmapped_genes=True, copy=False,
)
# Refuse a panel that barely maps, before anything reads a copy number from it.
# insitucnv's default gene-position reference is the *human* infercnvpy Maynard
# 2020 table, so a mouse panel matches only the symbols spelled identically in
# both nomenclatures — 8 of 5006 on the dataset that prompted this check. CNV
# inference reads copy number from runs of neighbouring genes along a
# chromosome, so a result built from that many genes is noise rather than a weak
# signal, and nothing downstream says so: the run completes, clusters, and the
# first symptom appears several steps later.
#
# The decision is the mapped *fraction*, not the species: a mouse panel given a
# mouse annotation table should run, and a human panel against a broken
# reference should not. Hand-written because no library API covers it.
if adata_cnv.n_vars < max(1, _n_panel * $min_mapped_fraction):
raise RuntimeError(
f"Only {adata_cnv.n_vars} of {_n_panel} panel genes have genomic "
f"coordinates ({adata_cnv.n_vars / _n_panel:.1%}). CNV inference needs a "
f"gene-position reference matching this panel's species, and the default "
f"one is human — a mouse panel needs an annotation table with gene_name, "
f"chromosome, start and end columns."
)
#--- block arrow_shim # frozen: not customisable
_old_infer = pd.options.future.infer_string
pd.options.future.infer_string = False
try:
for _attr in ('obs', 'var'):
_df = getattr(adata_cnv, _attr).copy()
if pd.api.types.is_string_dtype(_df.index):
_df.index = pd.Index(_df.index.to_numpy(dtype=object))
for _col in _df.columns:
if isinstance(_df[_col].dtype, pd.CategoricalDtype):
_cats = _df[_col].cat.categories
if pd.api.types.is_string_dtype(_cats):
_df[_col] = _df[_col].cat.rename_categories(
dict(zip(_cats, _cats.astype(object)))
)
elif pd.api.types.is_string_dtype(_df[_col]):
_df[_col] = _df[_col].to_numpy(dtype=object)
setattr(adata_cnv, _attr, _df)
finally:
pd.options.future.infer_string = _old_infer
#--- block tail
run_infercnv(
adata_cnv, reference_key=$reference_obs_key,
reference_categories=$reference_categories,
window_size=$window_size, step=$step, lfc_clip=$lfc_clip,
calculate_gene_values=True, copy=False,
)
compute_cnv_neighbors(adata_cnv, copy=False)
cnv_cluster_keys = cluster_cnv_resolutions(
adata_cnv, [$resolution], key_prefix='cnv_leiden_res',
dendrogram=False, copy=False,
)
_cnv_ids = (adata_cnv.obs['cell_id'].values
if 'cell_id' in adata_cnv.obs.columns else adata_cnv.obs_names)
cnv_clusters = pd.Series(
adata_cnv.obs[cnv_cluster_keys[0]].values,
index=_cnv_ids, name=cnv_cluster_keys[0],
)
_X_cnv = adata_cnv.obsm['X_cnv']
# abs().mean(axis=1) works natively on a CSR matrix, so the whole
# n_cells x n_bins CNV matrix never has to be densified for a row mean.
cnv_score = pd.Series(
np.asarray(np.abs(_X_cnv).mean(axis=1)).ravel(),
index=_cnv_ids, name='cnv_score',
)Run by: Gene Correlation
Scatters two genes cell by cell, annotates Pearson and Spearman coefficients with p-values, and saves the figure.
The expr.fraction block calls sc.pp.calculate_qc_metrics with inplace=False on purpose. With inplace=True scanpy overwrites obs['total_counts'] — one of Xenium's own columns — so the template would mutate adata on its way to reading a total. A template must not change the data it is describing.
Contract
| Parameter | Type | Required |
|---|---|---|
gene_a |
str |
yes |
gene_b |
str |
yes |
norm_label |
str |
yes |
xlabel |
str |
yes |
ylabel |
str |
yes |
title_prefix |
str |
yes |
paths |
list |
yes |
clustering |
str |
no |
selected |
list |
no |
-
Requires:
Path,adata,adata_norm,plt,sc -
Outputs:
fig,x,pr,pp,sr,sp -
Blocks:
head,expr.raw,expr.fraction,expr.log1p_cpm,filter,stats,title.plain,title.filtered,tail
Variants — 6 assemblies:
Six assemblies: three normalisations (raw counts, fraction of the cell's total, log1p CPM from adata_norm) times whether a cluster subset is selected, which adds the filter block and changes the title to say the count is filtered.
head + expr.raw + stats + title.plain + tailhead + expr.raw + filter + stats + title.filtered + tailhead + expr.fraction + stats + title.plain + tailhead + expr.fraction + filter + stats + title.filtered + tailhead + expr.log1p_cpm + stats + title.plain + tailhead + expr.log1p_cpm + filter + stats + title.filtered + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# Gene correlation ($norm_label): $gene_a vs $gene_b
from scipy.stats import pearsonr, spearmanr
#--- block expr.raw
_expr = sc.get.obs_df(adata, keys=[$gene_a, $gene_b])
x = _expr[$gene_a].to_numpy(dtype='float32')
y = _expr[$gene_b].to_numpy(dtype='float32')
#--- block expr.fraction
_expr = sc.get.obs_df(adata, keys=[$gene_a, $gene_b])
# Totals from scanpy's QC helper rather than a hand-rolled row sum.
# inplace=False deliberately: inplace=True writes obs['total_counts'], which is
# one of Xenium's own columns, so it would overwrite the raw value on adata.
_qc_obs, _ = sc.pp.calculate_qc_metrics(adata, percent_top=None, log1p=False,
inplace=False)
_totals = _qc_obs['total_counts'].to_numpy(dtype='float64')
_totals[_totals == 0] = 1
x = (_expr[$gene_a].to_numpy() / _totals).astype('float32')
y = (_expr[$gene_b].to_numpy() / _totals).astype('float32')
#--- block expr.log1p_cpm
_expr = sc.get.obs_df(adata_norm, keys=[$gene_a, $gene_b])
x = _expr[$gene_a].to_numpy(dtype='float32')
y = _expr[$gene_b].to_numpy(dtype='float32')
#--- block filter
_sel = adata.obs[$clustering].astype(str).isin($selected).to_numpy()
x, y = x[_sel], y[_sel]
#--- block stats
pr, pp = pearsonr(x, y)
sr, sp = spearmanr(x, y)
_p = lambda p: f'{p:.2e}' if p < 0.001 else f'{p:.4f}'
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(x, y, s=1, alpha=0.3, rasterized=True, color='#1f77b4')
ax.set_xlabel($xlabel)
ax.set_ylabel($ylabel)
#--- block title.plain
ax.set_title($title_prefix + f' [n={len(x):,}]')
#--- block title.filtered
ax.set_title($title_prefix + f' [n={len(x):,} (filtered)]')
#--- block tail
ax.text(
0.03, 0.97,
f'Pearson r = {pr:.3f}, p = {_p(pp)}\nSpearman \u03c1 = {sr:.3f}, p = {_p(sp)}',
transform=ax.transAxes, va='top', ha='left', fontsize=9,
bbox={'boxstyle': 'round,pad=0.3', 'fc': 'white', 'alpha': 0.7},
)
fig.tight_layout()
for _path in $paths:
Path(_path).parent.mkdir(parents=True, exist_ok=True)
fig.savefig(_path, dpi=300, bbox_inches='tight')Run by: Markers
Draws one scanpy marker-gene figure — dotplot, heatmap, matrixplot, tracksplot or correlation matrix — for a set of genes grouped by a clustering, and saves it.
Contract
| Parameter | Type | Required |
|---|---|---|
plot_name |
str |
yes |
groupby |
str |
yes |
markers |
dict |
no |
paths |
list |
yes |
categories |
dict |
no |
-
Requires:
Path,adata,adata_norm,pd,plt,sc -
Outputs:
fig -
Blocks:
head,relabel,call.dotplot,call.heatmap,call.matrixplot,call.tracksplot,call.correlation_matrix,save
Variants — 10 assemblies:
Ten assemblies: which of the five call.* plot types you picked, times whether the clusters have been given names (relabel). There used to be twenty — a second axis for whether the save carried dpi=150 — but the output format is no longer a per-plot choice: the save block loops over a paths list holding one file per format Preferences asks for.
head + call.dotplot + savehead + relabel + call.dotplot + savehead + call.heatmap + savehead + relabel + call.heatmap + savehead + call.matrixplot + savehead + relabel + call.matrixplot + savehead + call.tracksplot + savehead + relabel + call.tracksplot + savehead + call.correlation_matrix + savehead + relabel + call.correlation_matrix + save
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# $plot_name: $groupby
adata_norm.obs[$groupby] = adata.obs[$groupby].values
#--- block relabel
# .map() rather than .cat.rename_categories(): naming two clusters the same
# thing is a request to merge them, and rename_categories refuses it outright
# ("Categorical categories must be unique"). dict.fromkeys dedupes the names
# while keeping cluster order.
_display = $categories
adata_norm.obs[$groupby] = pd.Categorical(
adata_norm.obs[$groupby].astype(str).map(_display),
categories=list(dict.fromkeys(_display.values())),
)
#--- block call.dotplot
sc.pl.dotplot(adata_norm, var_names=$markers, groupby=$groupby, show=False)
#--- block call.heatmap
sc.pl.heatmap(adata_norm, var_names=$markers, groupby=$groupby, show=False)
#--- block call.matrixplot
sc.pl.matrixplot(adata_norm, var_names=$markers, groupby=$groupby, show=False)
#--- block call.tracksplot
sc.pl.tracksplot(adata_norm, var_names=$markers, groupby=$groupby, show=False)
#--- block call.correlation_matrix
sc.tl.dendrogram(adata_norm, $groupby)
sc.pl.correlation_matrix(adata_norm, $groupby, show=False)
#--- block save
# One figure, every configured format — the viewer writes a PNG to look at and
# a PDF to publish, and the notebook writes the same files to the same places.
fig = plt.gcf()
for _path in $paths:
Path(_path).parent.mkdir(parents=True, exist_ok=True)
fig.savefig(_path, bbox_inches='tight', dpi=300)Run by: Rank Genes
Ranks marker genes per cluster with sc.tl.rank_genes_groups, and returns the result as a tidy DataFrame.
Note key_added=$rank_key. scanpy writes uns['rank_genes_groups'] in place, so ranking a second clustering used to overwrite the first — a session that ranked two clusterings ended holding markers for whichever ran last. Keying per clustering is also what lets the sc.pl.rank_genes_groups* plots reach a specific result, since they all take key=.
Contract
| Parameter | Type | Required |
|---|---|---|
groupby |
str |
yes |
method |
str |
yes |
n_genes |
int |
yes |
rank_key |
str |
yes |
-
Requires:
adata,adata_norm,sc -
Outputs:
rank_df -
Blocks:
main
Default source
# Rank genes: groupby=$groupby, method=$method, n_genes=$n_genes
# key_added, so a second ranking does not overwrite this one: scanpy writes
# uns['rank_genes_groups'] in place, and a session that ranked two clusterings
# ended holding markers for whichever ran last. Keying is also what lets the
# sc.pl.rank_genes_groups* plots below reach a specific result, since they all
# take key=.
adata_norm.obs[$groupby] = adata.obs[$groupby].values
sc.tl.rank_genes_groups(
adata_norm, groupby=$groupby, method=$method, n_genes=$n_genes,
key_added=$rank_key,
)
rank_df = sc.get.rank_genes_groups_df(adata_norm, group=None, key=$rank_key)Run by: UMAP
Draws the UMAP embedding as a publication figure — one panel per selected gene, each with its own colour scale, or a single panel coloured by a clustering with the labels drawn on the points.
embed.xenium reads analysis/umap/gene_expression_2_components/projection.csv, which is where the viewer's own UMAP window gets its coordinates. That matters more than it looks: recomputing the embedding with sc.tl.umap gives an equally valid layout that is not the one that was on screen, so the notebook would disagree with the session it claims to reproduce.
Contract
| Parameter | Type | Required |
|---|---|---|
color |
list |
yes |
cmap |
str |
no |
ncols |
int |
no |
paths |
list |
yes |
groupby |
str |
no |
categories |
dict |
no |
-
Requires:
Path,adata,adata_norm,data_path,pd,plt,sc -
Outputs:
fig -
Blocks:
embed.xenium,embed.recompute,relabel,color.genes,color.clusters,save
Variants — 6 assemblies:
Six assemblies from three choices: whether the dataset ships Xenium's UMAP or the embedding has to be recomputed (a Crop Dataset export has no analysis/ folder), whether the colour is a list of genes or a clustering, and whether that clustering's clusters have been given names.
embed.xenium + color.genes + saveembed.recompute + color.genes + saveembed.xenium + color.clusters + saveembed.xenium + relabel + color.clusters + saveembed.recompute + color.clusters + saveembed.recompute + relabel + color.clusters + save
Default source — by block; an assembly above picks which of these run, in that order
#--- block embed.xenium
# UMAP: $color
#
# Xenium ships its own UMAP, and it is the one the viewer draws. Reading those
# coordinates rather than recomputing an embedding is what makes this figure the
# figure that was on screen — sc.tl.umap would produce a different, equally
# valid layout, and the notebook would quietly disagree with the session.
_umap = pd.read_csv(
data_path / 'analysis' / 'umap' / 'gene_expression_2_components' / 'projection.csv',
index_col=0)
# Joined on obs['cell_id'], not on obs_names: spatialdata_io indexes the table
# positionally ('0', '1', '2', ...) and keeps the barcode in a column, while the
# projection is indexed by barcode. Reindexing on obs_names matches nothing at
# all — every coordinate comes back NaN and the figure is empty.
adata_norm.obsm['X_umap'] = (
_umap.reindex(adata_norm.obs['cell_id']).to_numpy(dtype='float32'))
#--- block embed.recompute
# UMAP: $color
#
# No analysis/ folder — a Crop Dataset export has none — so compute one.
sc.pp.neighbors(adata_norm)
sc.tl.umap(adata_norm, random_state=0)
#--- block relabel
# .map() rather than .cat.rename_categories(): two clusters may deliberately
# carry the same display name, which is simply a request to merge them, but
# rename_categories refuses it outright ("Categorical categories must be
# unique"). dict.fromkeys dedupes the names while keeping cluster order, so the
# legend still reads in the order the clusters are numbered.
_display = $categories
adata_norm.obs[$groupby] = pd.Categorical(
adata.obs[$groupby].astype(str).map(_display),
categories=list(dict.fromkeys(_display.values())),
)
#--- block color.genes
# One panel per gene, each with its own colour bar; scanpy lays out the grid.
fig = sc.pl.umap(adata_norm, color=$color, cmap=$cmap, ncols=$ncols,
show=False, return_fig=True)
#--- block color.clusters
fig = sc.pl.umap(adata_norm, color=$color, legend_loc='on data',
show=False, return_fig=True)
#--- block save
for _path in $paths:
Path(_path).parent.mkdir(parents=True, exist_ok=True)
fig.savefig(_path, dpi=300, bbox_inches='tight')Run by: Co-occurrence
Squidpy's co-occurrence score across increasing radii — how the probability of finding one cluster near another changes with distance.
Contract
| Parameter | Type | Required |
|---|---|---|
cluster_key |
str |
yes |
interval |
int |
yes |
-
Requires:
adata,adata_norm,sq - Outputs: nothing (terminal step)
-
Blocks:
main
Default source
# Co-occurrence: $cluster_key (interval=$interval)
adata_norm.obs[$cluster_key] = adata.obs[$cluster_key].values
sq.gr.co_occurrence(adata_norm, cluster_key=$cluster_key, interval=$interval)Run by: Ligand-Receptor
Squidpy's ligand-receptor permutation test over the omnipath interaction database, between the clusters of a chosen clustering.
Contract
| Parameter | Type | Required |
|---|---|---|
cluster_key |
str |
yes |
n_perms |
int |
yes |
threshold |
float |
yes |
seed |
int |
yes |
include |
list |
no |
resources |
list |
no |
-
Requires:
adata,adata_norm,sq -
Outputs:
ligrec_res -
Blocks:
head,include,resources,tail
Variants — 4 assemblies:
Four assemblies: whether specific interaction datasets are included, and whether resources are restricted — each adds its own block that fills in interactions_params.
head + tailhead + resources + tailhead + include + tailhead + include + resources + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# Ligand-receptor: $cluster_key (n_perms=$n_perms)
from omnipath.constants import InteractionDataset
adata_norm.obs[$cluster_key] = adata.obs[$cluster_key].values
interactions_params = {}
#--- block include
interactions_params['include'] = tuple(
InteractionDataset[_n] for _n in $include
)
#--- block resources
interactions_params['resources'] = $resources
#--- block tail
ligrec_res = sq.gr.ligrec(
adata_norm, cluster_key=$cluster_key, n_perms=$n_perms,
threshold=$threshold, seed=$seed, use_raw=False, copy=True,
transmitter_params={'categories': 'ligand'},
receiver_params={'categories': 'receptor'},
interactions_params=interactions_params,
)Run by: Neighborhood Enrichment
Squidpy's neighbourhood-enrichment permutation test: which cluster pairs are adjacent more or less often than chance, as z-scores.
Needs the graph from spatial_neighbors first; the tab calls ctx.ensure_spatial_neighbors(k) before running this, and the provenance graph records the dependency.
Contract
| Parameter | Type | Required |
|---|---|---|
cluster_key |
str |
yes |
uns_key |
str |
yes |
n_perms |
int |
yes |
seed |
int |
yes |
-
Requires:
adata,adata_norm,sq - Outputs: nothing (terminal step)
-
Blocks:
main
Default source
# Neighborhood enrichment: $cluster_key (n_perms=$n_perms)
adata_norm.obs[$cluster_key] = adata.obs[$cluster_key].values
sq.gr.nhood_enrichment(
adata_norm, cluster_key=$cluster_key, n_perms=$n_perms, seed=$seed,
)
nhood_zscore = adata_norm.uns[$uns_key]['zscore']Run by: ROI Analysis
Differential expression between the regions you drew: assigns cells to ROIs, then ranks genes across them.
The centroid-to-ROI test goes through sd.models.PointsModel.parse and sd.polygon_query rather than a hand-rolled point-in-polygon loop. adata.obsm['spatial'] is in microns and the ROIs are in pixels; the scale between them is declared as a spatialdata transformation instead of applied by hand, which is what keeps the notebook's coordinate conventions honest.
Contract
| Parameter | Type | Required |
|---|---|---|
method |
str |
yes |
pixel_size |
float |
yes |
clustering |
str |
no |
selected |
list |
no |
-
Requires:
adata,np,pd,roi_polygons,sc -
Outputs:
roi_deg_df,roi_adata -
Blocks:
head,filter,loop_head,loop_filter,tail
Variants — 2 assemblies:
Two assemblies, differing in whether the analysis is restricted to selected clusters.
head + loop_head + tailhead + filter + loop_head + loop_filter + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# ROI differential expression (method=$method)
# Cell centroids as a spatialdata points element. adata.obsm['spatial'] is in
# microns and the ROIs are in pixels; the scale between them is *declared* as a
# transformation rather than applied by hand, so sd.polygon_query does the
# point-in-polygon test in the ROIs' own frame.
roi_cells = sd.models.PointsModel.parse(
pd.DataFrame({
'x': adata.obsm['spatial'][:, 0],
'y': adata.obsm['spatial'][:, 1],
'cell_index': np.arange(adata.n_obs),
}),
coordinates={'x': 'x', 'y': 'y'},
transformations={'global': sd.transformations.Scale(
[1 / $pixel_size, 1 / $pixel_size], axes=('x', 'y'))},
)
roi_region = np.full(adata.n_obs, '', dtype=object)
#--- block filter
# Cluster filter: cells must be inside an ROI *and* in the selected clusters
cluster_mask = adata.obs[$clustering].astype(str).isin($selected).to_numpy()
#--- block loop_head
for _i, _poly in enumerate(roi_polygons):
# polygon_query returns None, not an empty frame, for an ROI with no cells.
_hit = sd.polygon_query(roi_cells, _poly, target_coordinate_system='global')
_idx = (_hit.compute()['cell_index'].to_numpy() if _hit is not None
else np.empty(0, dtype=int))
#--- block loop_filter
_idx = _idx[cluster_mask[_idx]]
#--- block tail
# Last ROI wins where two overlap: a cell carries one region label.
roi_region[_idx] = f'Region {_i + 1}'
roi_adata = adata[roi_region != ''].copy()
roi_adata.obs['roi_region'] = pd.Categorical(roi_region[roi_region != ''])
# Normalised here rather than reusing adata_norm on purpose: this step is
# self-contained, and normalising the ROI subset is not the same computation as
# subsetting a globally normalised matrix.
sc.pp.normalize_total(roi_adata, target_sum=1e4)
sc.pp.log1p(roi_adata)
sc.tl.rank_genes_groups(
roi_adata, 'roi_region', method=$method, reference='rest', key_added=$method,
)
roi_deg_df = sc.get.rank_genes_groups_df(roi_adata, group=None, key=$method)Run by: ROI Analysis
Writes the per-cell ROI expression table produced by roi.expression to CSV.
A terminal step: it produces a file, not a value, so nothing depends on it. Its preview renders the filename the save dialog would propose, and the preview header says the path is a sample — a value that cannot come from a widget is named rather than faked.
Contract
| Parameter | Type | Required |
|---|---|---|
gene |
str |
yes |
path |
str |
yes |
-
Requires:
roi_expr_cells - Outputs: nothing (terminal step)
-
Blocks:
main
Default source
# Export ROI per-cell expression of $gene
roi_expr_cells.to_csv($path, index=False)Run by: ROI Analysis
Per-region expression of a single gene, with summary statistics per ROI and pairwise Welch t-tests between every pair of regions, Benjamini-Hochberg corrected.
The Welch + BH block is hand-rolled numpy/scipy rather than a library call, and says so in a comment. That is allowed where no API covers the case — the rule is to prefer the library where one exists, not to avoid writing statistics.
Contract
| Parameter | Type | Required |
|---|---|---|
gene |
str |
yes |
pixel_size |
float |
yes |
clustering |
str |
no |
selected |
list |
no |
-
Requires:
adata,np,pd,roi_polygons,sc -
Outputs:
roi_expr_cells,roi_expr_stats,roi_expr_tests -
Blocks:
head,filter,loop_head,loop_filter,tail
Variants — 2 assemblies:
Two assemblies, as for roi.deg: with or without a cluster filter.
head + loop_head + tailhead + filter + loop_head + loop_filter + tail
Default source — by block; an assembly above picks which of these run, in that order
#--- block head
# ROI expression of $gene, per drawn region
from itertools import combinations
from scipy import stats
# Cell centroids as a spatialdata points element. adata.obsm['spatial'] is in
# microns and the ROIs are in pixels; the scale between them is *declared* as a
# transformation rather than applied by hand, so sd.polygon_query does the
# point-in-polygon test in the ROIs' own frame.
roi_cells = sd.models.PointsModel.parse(
pd.DataFrame({
'x': adata.obsm['spatial'][:, 0],
'y': adata.obsm['spatial'][:, 1],
'cell_index': np.arange(adata.n_obs),
}),
coordinates={'x': 'x', 'y': 'y'},
transformations={'global': sd.transformations.Scale(
[1 / $pixel_size, 1 / $pixel_size], axes=('x', 'y'))},
)
roi_region = np.zeros(adata.n_obs, dtype=int) # 0 = outside every ROI
#--- block filter
# Cluster filter: cells must be inside an ROI *and* in the selected clusters
cluster_mask = adata.obs[$clustering].astype(str).isin($selected).to_numpy()
#--- block loop_head
for _i, _poly in enumerate(roi_polygons):
# polygon_query returns None, not an empty frame, for an ROI with no cells.
_hit = sd.polygon_query(roi_cells, _poly, target_coordinate_system='global')
_idx = (_hit.compute()['cell_index'].to_numpy() if _hit is not None
else np.empty(0, dtype=int))
#--- block loop_filter
_idx = _idx[cluster_mask[_idx]]
#--- block tail
# Last ROI wins where two overlap: a cell belongs to one region.
roi_region[_idx] = _i + 1
_cells = sc.get.obs_df(
adata, keys=[$gene], obsm_keys=[('spatial', 0), ('spatial', 1)],
)
roi_expr_cells = pd.DataFrame({
'region_id': roi_region,
'cell_id': (adata.obs['cell_id'].to_numpy() if 'cell_id' in adata.obs
else adata.obs_names.to_numpy()),
'x_centroid_um': _cells['spatial-0'].to_numpy(),
'y_centroid_um': _cells['spatial-1'].to_numpy(),
'expression': _cells[$gene].to_numpy(),
})
# Region-grouped, as the exported CSV has always been; stable, so cells keep
# their order within a region.
roi_expr_cells = (roi_expr_cells[roi_expr_cells['region_id'] > 0]
.sort_values('region_id', kind='stable')
.reset_index(drop=True))
roi_expr_stats = (
roi_expr_cells.groupby('region_id')['expression']
.agg(['count', 'mean', 'median', 'std', 'min', 'max'])
.reindex(range(1, len(roi_polygons) + 1))
)
roi_expr_stats['count'] = roi_expr_stats['count'].fillna(0).astype(int)
# Pairwise Welch's t-tests between regions, Benjamini-Hochberg corrected.
# Hand-rolled deliberately: scanpy has no per-region two-sample test, and
# scipy's false_discovery_control is already the right call for the correction.
_groups = [(_r, _g['expression'].to_numpy())
for _r, _g in roi_expr_cells.groupby('region_id') if len(_g) >= 2]
_tests = []
for (_r1, _e1), (_r2, _e2) in combinations(_groups, 2):
_t, _p = stats.ttest_ind(_e1, _e2, equal_var=False)
_tests.append({'region_1': _r1, 'region_2': _r2, 't': _t, 'p': _p})
roi_expr_tests = pd.DataFrame(_tests, columns=['region_1', 'region_2', 't', 'p'])
roi_expr_tests['p_adj'] = (
stats.false_discovery_control(roi_expr_tests['p'], method='bh')
if len(roi_expr_tests) > 1 else roi_expr_tests['p']
)Run by: ROI Analysis
Inlines the polygons you drew in the viewer as literal coordinates, so the notebook reproduces the exact regions without needing the napari shapes layer.
Two deliberate details. napari stores shapes as N×2 (y, x) pixel arrays, and the flip to (x, y) happens once, here — an earlier version flipped and unflipped on consecutive lines, so the comment asserted a convention the code reversed. And make_valid, not buffer(0): on a self-intersecting ROI buffer(0) silently deletes a lobe rather than repairing it.
Contract
| Parameter | Type | Required |
|---|---|---|
polygons |
list |
yes |
-
Requires:
np -
Outputs:
roi_polygons -
Blocks:
main
Default source
# ROI polygons drawn in the viewer. napari stores shapes as Nx2 (y, x) pixel
# arrays; they are flipped to (x, y) once here, so every consumer works in the
# same frame as sdata.shapes['rois'] instead of flipping back and forth.
# make_valid, not buffer(0): buffer(0) silently *deletes* a lobe of a
# self-intersecting ROI rather than repairing it.
from shapely import make_valid
from shapely.geometry import Polygon
roi_polygons = [make_valid(Polygon(np.asarray(_p)[:, ::-1])) for _p in $polygons]Reference
Cells
Genes
Spatial
- ROI Analysis
- Ligand-Receptor
- Neighborhood Enrichment
- Co-occurrence
- Spatial Domains
- Annot Nhood
- Annot Distance
Images
Tools
Tutorials
- Getting Started
- Clustering and DEG
- H&E Registration
- ARMS Overlay
- ROI Analysis
- Annotations
- Recovering a Cache