Latent motifs: recurrent cellular neighbourhoods¶

Co-occurrence works pair by pair. But tissue organisation is rarely pairwise — a tumour-immune boundary is a combination of cell types that recurs across the slide.

The recipe here is always the same:

1. Describe each cell's neighbourhood as a vector.
2. Reduce those vectors to a handful of recurring patterns.
3. Map the patterns back onto the tissue.

Three ways to do step 2, with different trade-offs.

In [1]:

from pathlib import Path

import matplotlib.pyplot as plt
import pandas as pd

import scimappro as sp

# The demo data is not shipped with the package (it is 247 MB). Point DATA at
# wherever you unpacked it; this looks in the repository root first, then in the
# path the docs are built from.
DATA = next(
    path for path in (Path("example_data"), Path("../../../example_data"))
    if path.exists()
)
DATA

from pathlib import Path import matplotlib.pyplot as plt import pandas as pd import scimappro as sp # The demo data is not shipped with the package (it is 247 MB). Point DATA at # wherever you unpacked it; this looks in the repository root first, then in the # path the docs are built from. DATA = next( path for path in (Path("example_data"), Path("../../../example_data")) if path.exists() ) DATA

Out[1]:

PosixPath('../../../example_data')

In [2]:

import anndata as ad

adata = ad.read_h5ad(DATA / "adata_scimap.h5ad")
adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False)
adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), verbose=False)
adata.obs["phenotype"].value_counts()

import anndata as ad adata = ad.read_h5ad(DATA / "adata_scimap.h5ad") adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False) adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), verbose=False) adata.obs["phenotype"].value_counts()

Out[2]:

phenotype
ECAD+                  7112
Other myeloid cells    2419
Dendritic cells         863
SMA+                    509
Treg                    216
Unknown                  46
NK cells                 35
Immune                    1
Name: count, dtype: int64

Step 1: the neighbourhood matrix¶

[sp.tl.neighCount][scimappro.tl.neighCount] gives every cell a vector of what surrounds it — the proportion of its neighbours that are each cell type.

In [3]:

adata = sp.tl.neighCount(adata, phenotype="phenotype", method="radius",
                         radius=30, verbose=False)
adata.uns["neighCount"].head()

adata = sp.tl.neighCount(adata, phenotype="phenotype", method="radius", radius=30, verbose=False) adata.uns["neighCount"].head()

Out[3]:

	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
quant_1	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0
quant_2	0.000000	0.600000	0.0	0.0	0.2	0.2	0.0	0.0
quant_3	0.000000	0.800000	0.0	0.0	0.2	0.0	0.0	0.0
quant_4	0.333333	0.666667	0.0	0.0	0.0	0.0	0.0	0.0
quant_5	0.000000	1.000000	0.0	0.0	0.0	0.0	0.0	0.0

Hard assignments: cluster the matrix¶

The most direct approach. mode="spatial" tells [sp.tl.cluster][scimappro.tl.cluster] to read uns[layer] rather than the expression matrix.

!!! note In mode="spatial", layer names a key in uns, not in layers. That is where the neighbourhood functions leave their output.

In [4]:

adata = sp.tl.cluster(adata, mode="spatial", layer="neighCount", method="kmeans",
                      k=6, label="rcn", verbose=False)
adata.obs["rcn"].value_counts()

adata = sp.tl.cluster(adata, mode="spatial", layer="neighCount", method="kmeans", k=6, label="rcn", verbose=False) adata.obs["rcn"].value_counts()

Out[4]:

rcn
2    4098
0    2383
3    1729
5    1247
1     969
4     775
Name: count, dtype: int64

In [5]:

sp.pl.spatialScatterPlot(adata, colorBy="rcn", s=2, figsize=(5, 5), fontSize=7)

sp.pl.spatialScatterPlot(adata, colorBy="rcn", s=2, figsize=(5, 5), fontSize=7)

What is each neighbourhood made of?¶

A cross-tabulation of neighbourhood against cell type is how you read off what a given RCN actually is.

In [6]:

composition = pd.crosstab(adata.obs["rcn"], adata.obs["phenotype"], normalize="index")
composition.round(2)

composition = pd.crosstab(adata.obs["rcn"], adata.obs["phenotype"], normalize="index") composition.round(2)

Out[6]:

phenotype	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
rcn
0	0.01	0.75	0.0	0.00	0.18	0.04	0.01	0.00
1	0.61	0.16	0.0	0.00	0.18	0.02	0.03	0.00
2	0.01	0.89	0.0	0.00	0.07	0.03	0.01	0.00
3	0.04	0.48	0.0	0.00	0.39	0.06	0.03	0.00
4	0.03	0.52	0.0	0.03	0.17	0.18	0.02	0.05
5	0.10	0.22	0.0	0.00	0.59	0.05	0.05	0.00

Soft assignments: LDA¶

Hard clustering forces every cell into one neighbourhood. Real tissue has boundaries and gradients where a cell genuinely belongs to two.

[sp.tl.neighLDA][scimappro.tl.neighLDA] treats each neighbourhood as a document and each cell type as a word, and gives every cell a distribution over motifs.

In [7]:

adata = sp.tl.neighLDA(adata, phenotype="phenotype", neighbor_method="radius",
                       radius=30, numMotifs=6, label="lda", verbose=False)
adata.uns["lda"].head()

adata = sp.tl.neighLDA(adata, phenotype="phenotype", neighbor_method="radius", radius=30, numMotifs=6, label="lda", verbose=False) adata.uns["lda"].head()

Images:   0%|          | 0/1 [00:00<?, ?it/s]

Images: 100%|██████████| 1/1 [00:00<00:00, 61.59it/s]

Out[7]:

	Motif_0	Motif_1	Motif_2	Motif_3	Motif_4	Motif_5
quant_1	0.582967	0.083334	0.083335	0.083520	0.083502	0.083343
quant_2	0.023829	0.558942	0.196270	0.023810	0.023843	0.173306
quant_3	0.023829	0.708414	0.196280	0.023810	0.023842	0.023825
quant_4	0.273638	0.512719	0.020885	0.020833	0.020862	0.151062
quant_5	0.033338	0.833180	0.033416	0.033333	0.033380	0.033353

numMotifs is the knob that matters: too few and distinct neighbourhoods merge, too many and they fragment.

Cluster the motif weights if you do want discrete regions in the end:

In [8]:

adata = sp.tl.cluster(adata, mode="spatial", layer="lda", method="kmeans", k=6,
                      label="lda_cluster", verbose=False)
sp.pl.spatialScatterPlot(adata, colorBy="lda_cluster", s=2, figsize=(5, 5), fontSize=7)

adata = sp.tl.cluster(adata, mode="spatial", layer="lda", method="kmeans", k=6, label="lda_cluster", verbose=False) sp.pl.spatialScatterPlot(adata, colorBy="lda_cluster", s=2, figsize=(5, 5), fontSize=7)

Soft assignments: NMF¶

[sp.tl.neighNMF][scimappro.tl.neighNMF] does the same job through non-negative matrix factorisation: faster than LDA, deterministic given a seed, and it reports its components explicitly so you can read what each one is. There is no scimap equivalent.

In [9]:

adata = sp.tl.neighNMF(adata, phenotype="phenotype", method="radius", radius=30,
                       nmfComponents=6, label="nmf", verbose=False)
adata.uns["nmf_nmf_components"].round(3)

adata = sp.tl.neighNMF(adata, phenotype="phenotype", method="radius", radius=30, nmfComponents=6, label="nmf", verbose=False) adata.uns["nmf_nmf_components"].round(3)

Out[9]:

array([[0.000e+00, 8.820e+00, 0.000e+00, 8.000e-03, 0.000e+00, 0.000e+00,
        0.000e+00, 0.000e+00],
       [0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 6.297e+00, 0.000e+00,
        0.000e+00, 0.000e+00],
       [4.762e+00, 0.000e+00, 0.000e+00, 2.000e-03, 0.000e+00, 0.000e+00,
        0.000e+00, 0.000e+00],
       [0.000e+00, 1.000e-03, 0.000e+00, 0.000e+00, 0.000e+00, 3.682e+00,
        0.000e+00, 1.000e-03],
       [0.000e+00, 0.000e+00, 0.000e+00, 1.800e-02, 0.000e+00, 0.000e+00,
        2.706e+00, 0.000e+00],
       [0.000e+00, 0.000e+00, 0.000e+00, 6.610e-01, 0.000e+00, 0.000e+00,
        0.000e+00, 2.335e+00]])

One row per component, one column per cell type — so each row reads directly as "this neighbourhood is mostly X with some Y".

The per-cell weights:

In [10]:

adata.uns["nmf"].head().round(3)

adata.uns["nmf"].head().round(3)

Out[10]:

	NMF_1	NMF_2	NMF_3	NMF_4	NMF_5	NMF_6
quant_1	0.000	0.000	0.00	0.000	0.0	0.0
quant_2	0.068	0.032	0.00	0.054	0.0	0.0
quant_3	0.091	0.032	0.00	0.000	0.0	0.0
quant_4	0.076	0.000	0.07	0.000	0.0	0.0
quant_5	0.113	0.000	0.00	0.000	0.0	0.0

In [11]:

adata = sp.tl.cluster(adata, mode="spatial", layer="nmf", method="kmeans", k=6,
                      label="nmf_cluster", verbose=False)
sp.pl.spatialScatterPlot(adata, colorBy="nmf_cluster", s=2, figsize=(5, 5), fontSize=7)

adata = sp.tl.cluster(adata, mode="spatial", layer="nmf", method="kmeans", k=6, label="nmf_cluster", verbose=False) sp.pl.spatialScatterPlot(adata, colorBy="nmf_cluster", s=2, figsize=(5, 5), fontSize=7)

Building neighbourhoods without cell types¶

inferenceMode="marker_driven" counts neighbours by the markers they express rather than by a phenotype label, so you can find neighbourhoods before you have called any cell types at all.

In [12]:

markerDriven = sp.tl.neighCount(
    adata.copy(),
    inferenceMode="marker_driven",
    layer=None,
    method="knn",
    knn=10,
    label="neighCount_markers",
    verbose=False,
)
markerDriven.uns["neighCount_markers"].head().round(3)

markerDriven = sp.tl.neighCount( adata.copy(), inferenceMode="marker_driven", layer=None, method="knn", knn=10, label="neighCount_markers", verbose=False, ) markerDriven.uns["neighCount_markers"].head().round(3)

Out[12]:

	CD11B	CD16	CD45	CD57	ECAD	ELANE	FOXP3	NCAM	SMA
quant_1	0.155	0.121	0.069	0.017	0.155	0.155	0.052	0.155	0.121
quant_2	0.158	0.123	0.123	0.000	0.158	0.158	0.000	0.158	0.123
quant_3	0.167	0.111	0.111	0.000	0.167	0.167	0.000	0.167	0.111
quant_4	0.148	0.066	0.131	0.016	0.148	0.148	0.049	0.148	0.148
quant_5	0.170	0.094	0.113	0.000	0.170	0.170	0.000	0.170	0.113

Which to use¶

	Assignment	Speed	Reads back as
neighCount + cluster	hard	fastest	a cross-tabulation
neighLDA	soft	slowest	topic weights
neighNMF	soft	fast	explicit components

Start with neighCount + cluster. Move to NMF when boundaries matter, and to LDA when you want the probabilistic interpretation.

Next¶

• Neighbourhood lag — the same idea over marker expression rather than cell-type composition.
• Search patterns — find one specific neighbourhood rather than all of them.