Search patterns¶

You spot an interesting region on one slide. Where else does it occur?

[sp.tl.spatialSimilarityLookup][scimappro.tl.spatialSimilarityLookup] turns that into a query: mark a reference region, and every cell in the dataset is scored by how closely its neighbourhood resembles it.

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: mark a reference region¶

Any obs column naming regions will do. In a real analysis these come from [sp.helpers.addROI_omero][scimappro.helpers.addROI_omero] — regions you drew in OMERO — see Add ROIs.

Here we define one by thresholding a marker, which is what [sp.pl.addRoiScatter][scimappro.pl.addRoiScatter] does.

In [3]:

adata = sp.pl.addRoiScatter(
    adata, marker="ECAD", threshold=0.5, roiName="tumour", label="roi"
)
adata.obs["roi"].value_counts()

adata = sp.pl.addRoiScatter( adata, marker="ECAD", threshold=0.5, roiName="tumour", label="roi" ) adata.obs["roi"].value_counts()

Out[3]:

roi
tumour    11135
Other        66
Name: count, dtype: int64

Cells labelled 'Other' are treated as outside every ROI and are not used as queries.

In [4]:

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

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

Step 2: search¶

The neighbourhood-weighted expression is computed for every cell, the median lag vector of the reference region becomes the query, and every cell is scored by similarity to it.

In [5]:

adata = sp.tl.spatialSimilarityLookup(
    adata,
    roiColumn="roi",
    method="radius",
    radius=30,
    similarityThreshold=0.5,
    verbose=False,
)
[c for c in adata.obs.columns if c.startswith("spatialSimilarityLookup")]

adata = sp.tl.spatialSimilarityLookup( adata, roiColumn="roi", method="radius", radius=30, similarityThreshold=0.5, verbose=False, ) [c for c in adata.obs.columns if c.startswith("spatialSimilarityLookup")]

Out[5]:

['spatialSimilarityLookup_tumour']

One obs column per ROI, holding whether each cell cleared the threshold, plus the raw scores in layers.

In [6]:

adata.obs["spatialSimilarityLookup_tumour"].value_counts()

adata.obs["spatialSimilarityLookup_tumour"].value_counts()

Out[6]:

spatialSimilarityLookup_tumour
other             8325
similar_to_ROI    2876
Name: count, dtype: int64

In [7]:

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

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

Step 3: tune the threshold¶

Higher is stricter. reuseSimilarityMatrix skips recomputing the spatial lag — the expensive part — so re-thresholding is nearly free.

In [8]:

for threshold in (0.3, 0.5, 0.7):
    result = sp.tl.spatialSimilarityLookup(
        adata.copy(),
        roiColumn="roi",
        similarityThreshold=threshold,
        reuseSimilarityMatrix="spatialSimilarityLookup",
        label=f"lookup_{threshold}",
        verbose=False,
    )
    counts = result.obs[f"lookup_{threshold}_tumour"].value_counts().to_dict()
    print(f"threshold={threshold}  {counts}")

for threshold in (0.3, 0.5, 0.7): result = sp.tl.spatialSimilarityLookup( adata.copy(), roiColumn="roi", similarityThreshold=threshold, reuseSimilarityMatrix="spatialSimilarityLookup", label=f"lookup_{threshold}", verbose=False, ) counts = result.obs[f"lookup_{threshold}_tumour"].value_counts().to_dict() print(f"threshold={threshold} {counts}")

threshold=0.3  {'similar_to_ROI': 7412, 'other': 3789}

threshold=0.5  {'other': 8325, 'similar_to_ROI': 2876}

threshold=0.7  {'other': 10808, 'similar_to_ROI': 393}

Where the scores live¶

The per-cell similarity scores are in layers[label] — one column per marker of the lag matrix — so you can threshold or plot them yourself.

In [9]:

adata.layers["spatialSimilarityLookup"].shape

adata.layers["spatialSimilarityLookup"].shape

Out[9]:

(11201, 9)

Next¶

• Add ROIs — getting hand-drawn regions into the object.
• Latent motifs — find recurring neighbourhoods without naming one first.