Co-occurrence analysis¶

Distance tells you how far apart two cell types are. It does not tell you whether that is surprising: in a tissue that is 60% tumour, everything is near tumour.

[sp.tl.spatialCooccurrence][scimappro.tl.spatialCooccurrence] answers the harder question by comparing observed adjacency to a permutation null — shuffle the phenotype labels, keep the geometry, and see how often each pair still ends up adjacent.

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

Running the test¶

permutation sets how many shuffles sit behind each p-value. Runtime is linear in it — use 100 while exploring, 1000 for a result you will publish.

In [3]:

adata = sp.tl.spatialCooccurrence(
    adata,
    phenotype="phenotype",
    method="radius",
    radius=30,
    permutation=1000,
    verbose=False,
)
adata.uns["spatialCooccurrence"].head(10)

adata = sp.tl.spatialCooccurrence( adata, phenotype="phenotype", method="radius", radius=30, permutation=1000, verbose=False, ) adata.uns["spatialCooccurrence"].head(10)

Out[3]:

	phenotype	neighbour_phenotype	quant	pvalue_quant
0	Dendritic cells	Dendritic cells	1.000000	0.000000e+00
1	Dendritic cells	ECAD+	-0.240656	0.000000e+00
2	Dendritic cells	NK cells	-0.000325	1.986503e-05
3	Dendritic cells	Other myeloid cells	-0.364641	9.926223e-28
4	Dendritic cells	SMA+	-0.030874	9.894845e-37
5	Dendritic cells	Treg	0.047774	1.485168e-04
6	Dendritic cells	Unknown	-0.000487	6.995860e-06
7	ECAD+	Dendritic cells	-0.047771	0.000000e+00
8	ECAD+	ECAD+	1.000000	0.000000e+00
9	ECAD+	NK cells	-0.001322	7.358747e-11

The result is a long table: one row per ordered pair of phenotypes, a score column named after the image, and a matching pvalue_<image> column.

The score is not a raw count. Counts are divided by the number of cells of the source type, rescaled so each row's largest value is 1, and signed by whether the pair occurred more or less often than expected. So +1 is strong attraction, -1 strong avoidance.

The heatmap¶

Pairs that fail the significance threshold are drawn grey rather than on the colour scale, so only real structure carries colour.

In [4]:

sp.pl.spatialCooccurrence(adata, pVal=0.05)

sp.pl.spatialCooccurrence(adata, pVal=0.05)

Attraction and avoidance only, ignoring magnitude:

In [5]:

sp.pl.spatialCooccurrence(adata, binaryView=True, pVal=0.05)

sp.pl.spatialCooccurrence(adata, binaryView=True, pVal=0.05)

The network view¶

Faster to read when the question is simply which cell types interact at all. An edge is drawn wherever the pair is significant.

In [6]:

sp.pl.spatialInteractionNetwork(adata, pVal=0.01)

sp.pl.spatialInteractionNetwork(adata, pVal=0.01)

Nearest neighbours instead of a radius¶

method="radius" asks "who is within 30 pixels". method="knn" asks "who are the 10 closest cells", regardless of how far away they are. The second is more robust when cell density varies a lot across the tissue.

In [7]:

knn = sp.tl.spatialCooccurrence(
    adata.copy(), phenotype="phenotype", method="knn", knn=10,
    permutation=1000, label="cooccurrence_knn", verbose=False,
)
sp.pl.spatialCooccurrence(knn, spatialCooccurrence="cooccurrence_knn")

knn = sp.tl.spatialCooccurrence( adata.copy(), phenotype="phenotype", method="knn", knn=10, permutation=1000, label="cooccurrence_knn", verbose=False, ) sp.pl.spatialCooccurrence(knn, spatialCooccurrence="cooccurrence_knn")

Focusing on a few cell types¶

In [8]:

sp.pl.spatialCooccurrence(
    adata,
    subsetPhenotype=["Treg", "NK cells", "Dendritic cells"],
    subsetNeighbourPhenotype=["ECAD+", "SMA+"],
)

sp.pl.spatialCooccurrence( adata, subsetPhenotype=["Treg", "NK cells", "Dendritic cells"], subsetNeighbourPhenotype=["ECAD+", "SMA+"], )

Getting the numbers out¶

In [9]:

scores = sp.pl.spatialCooccurrence(adata, returnData=True, show=False)
scores.round(2)

scores = sp.pl.spatialCooccurrence(adata, returnData=True, show=False) scores.round(2)

Out[9]:

neighbour_phenotype	Dendritic cells	ECAD+	NK cells	Other myeloid cells	SMA+	Treg	Unknown
phenotype
Dendritic cells	1.00	-0.24	-0.00	-0.36	-0.03	0.05	-0.00
ECAD+	-0.05	1.00	-0.00	-0.24	-0.05	-0.02	-0.00
NK cells	-0.03	-0.53	1.00	-0.03	NaN	-0.03	0.35
Other myeloid cells	-0.25	-0.80	-0.00	1.00	0.10	0.08	NaN
SMA+	-0.12	-1.00	NaN	0.61	0.33	-0.05	-0.00
Treg	0.41	-0.82	-0.00	1.00	-0.10	0.16	NaN
Unknown	-0.03	-0.32	0.25	NaN	-0.01	NaN	1.00

Next¶

• Proximity scores — quantify how much of a sample is in an interaction zone.
• Latent motifs — beyond pairs, to recurring multi-cell neighbourhoods.