Neighbourhood lag¶

[sp.tl.neighCount][scimappro.tl.neighCount] asks what kinds of cell surround this one. [sp.tl.neighExp][scimappro.tl.neighExp] asks what those neighbours express — the spatial lag of each marker.

The difference matters when the interesting signal is not captured by your phenotype labels: a tumour cell surrounded by PD-L1-high neighbours is in a different environment from one surrounded by PD-L1-low neighbours, even though both neighbourhoods are "tumour".

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

Computing the lag¶

Every cell gets the average expression of its neighbours, giving a cell x marker matrix in which each value describes the cell's surroundings rather than the cell itself.

In [3]:

adata = sp.tl.neighExp(adata, method="radius", radius=30, layer="raw", log=True,
                       verbose=False)
adata.uns["neighExp"].head().round(3)

adata = sp.tl.neighExp(adata, method="radius", radius=30, layer="raw", log=True, verbose=False) adata.uns["neighExp"].head().round(3)

Out[3]:

	ELANE	CD57	CD45	CD11B	SMA	CD16	ECAD	FOXP3	NCAM
quant_1	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
quant_2	7.225	5.161	6.873	6.941	6.371	6.293	7.384	5.733	6.920
quant_3	7.234	5.198	6.662	6.876	6.253	6.247	7.265	5.821	6.887
quant_4	7.217	5.506	6.789	6.951	7.544	5.973	7.002	5.817	6.802
quant_5	7.422	5.097	5.845	6.862	5.972	5.980	7.216	5.623	6.802

Clustering the environment¶

Same shape as the motif tutorial, but the input is expression rather than composition.

In [4]:

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

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

Out[4]:

expressionNeighbourhood
0    6131
5    3663
3    1287
2      48
1      38
4      34
Name: count, dtype: int64

In [5]:

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

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

Cell versus environment¶

The point of the lag is that the two differ. Cross-tabulating a cell's own phenotype against its expression neighbourhood shows how much:

In [6]:

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

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

Out[6]:

expressionNeighbourhood	0	1	2	3	4	5
phenotype
Dendritic cells	0.97	0.00	0.00	0.00	0.00	0.03
ECAD+	0.37	0.00	0.00	0.15	0.00	0.47
Immune	0.00	1.00	0.00	0.00	0.00	0.00
NK cells	0.00	0.03	0.57	0.14	0.23	0.03
Other myeloid cells	0.86	0.00	0.00	0.05	0.00	0.09
SMA+	0.77	0.01	0.00	0.12	0.00	0.10
Treg	0.88	0.00	0.00	0.05	0.00	0.06
Unknown	0.00	0.00	0.59	0.00	0.33	0.09

Rows that spread across several columns are cell types whose environment varies — which is exactly the signal that a phenotype label alone throws away.

Counts alongside expression¶

The two are complementary; nothing stops you computing both and comparing.

In [7]:

adata = sp.tl.neighCount(adata, phenotype="phenotype", method="radius",
                         radius=30, verbose=False)
adata = sp.tl.cluster(adata, mode="spatial", layer="neighCount", method="kmeans",
                      k=6, label="compositionNeighbourhood", verbose=False)

pd.crosstab(adata.obs["compositionNeighbourhood"],
            adata.obs["expressionNeighbourhood"])

adata = sp.tl.neighCount(adata, phenotype="phenotype", method="radius", radius=30, verbose=False) adata = sp.tl.cluster(adata, mode="spatial", layer="neighCount", method="kmeans", k=6, label="compositionNeighbourhood", verbose=False) pd.crosstab(adata.obs["compositionNeighbourhood"], adata.obs["expressionNeighbourhood"])

Out[7]:

expressionNeighbourhood	0	1	2	3	4	5
compositionNeighbourhood
0	1298	0	1	241	4	839
1	953	0	1	0	4	11
2	680	0	0	818	8	2592
3	1508	0	0	98	0	123
4	486	38	46	91	18	96
5	1206	0	0	39	0	2

Nearest neighbours instead of a radius¶

With variable cell density, knn gives every cell the same number of neighbours, which keeps the lag comparable across sparse and dense regions.

In [8]:

knnLag = sp.tl.neighExp(adata.copy(), method="knn", knn=10, layer="raw",
                        label="neighExp_knn", verbose=False)
knnLag.uns["neighExp_knn"].head().round(3)

knnLag = sp.tl.neighExp(adata.copy(), method="knn", knn=10, layer="raw", label="neighExp_knn", verbose=False) knnLag.uns["neighExp_knn"].head().round(3)

Out[8]:

	ELANE	CD57	CD45	CD11B	SMA	CD16	ECAD	FOXP3	NCAM
quant_1	7.269	5.823	6.066	6.617	7.023	6.281	6.890	5.550	6.763
quant_2	7.272	5.215	6.710	6.924	6.601	6.277	7.338	5.757	6.893
quant_3	7.271	5.219	6.619	6.884	6.422	6.248	7.290	5.797	6.882
quant_4	7.213	5.446	6.689	6.964	7.332	5.996	7.006	5.826	6.809
quant_5	7.393	5.199	6.251	6.886	6.117	6.123	7.307	5.703	6.842

Next¶

• Latent motifs — the composition-based counterpart.
• Search patterns — the lag matrix is also what powers neighbourhood similarity search.