04. Spatial domains, and how a cell's state depends on where it is¶
Author: MD Shakhaowat Hossain, PhD Candidate, Tulane University School of Medicine · Project website · Code
Why: cell types say what is in the tissue; domains say how it is organised: tumour nests, tumour borders, stroma, vessels, lymphoid aggregates. The Step 7 niches did this with k-means on neighbourhood composition, a simple baseline. Here two methods built for spatial transcriptomics are applied, compared with that baseline, and then used to ask a central question of spatial biology: does the same cell type change its expression depending on where it sits?
Methods
- BANKSY (Singhal et al., Nature Genetics 2024): each cell's expression is augmented with the average expression of its spatial neighbours; λ sets the balance (λ = 0.8 is the recommended domain-finding mode). Then PCA and Leiden.
- CellCharter (Varrone et al., Nature Genetics 2024): each cell's embedding is concatenated with the mean embeddings of its 1st-, 2nd- and 3rd-order spatial neighbours, then clustered with a Gaussian mixture. The number of domains is chosen by stability across repeated runs.
- Pseudobulk differential expression with PyDESeq2, comparing one macrophage type between domains.
| Part | Content |
|---|---|
| A | Data |
| B | BANKSY, and how stable its domains are across parameters |
| C | CellCharter, with automatic choice of the number of domains |
| D | Compare the three domain sets (k-means niches, BANKSY, CellCharter) |
| E | Same cell type, different domain: pseudobulk differential expression |
| F | Export for Xenium Explorer |
Setup: imports and output folder¶
Purpose: load the shared ist_analysis package and define where results are written, so paths and constants match the core notebook and the Snakemake pipeline.
Approach: import domains (as dm) and utils (as xu) from ist_analysis, check the kernel, and create the output folder with xu.output_dir("spatial_domains").
Expected output: the output folder, printed as a ~/... path.
Note: this notebook reads the Step 9 object from the core notebook and the Proseg export from notebook 02; run both first. UserWarnings are silenced to keep BANKSY and CellCharter output readable.
from ist_analysis import domains as dm
from ist_analysis import utils as xu # paths and constants shared with the core workflow and the pipeline
xu.check_kernel()
import contextlib
import io
import logging
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import scipy.sparse as sp
from IPython.display import display
from sklearn.metrics import adjusted_rand_score as ari
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)
OUT_DIR = xu.output_dir("spatial_domains")
print(f"Outputs go to {xu.display_path(OUT_DIR)}")
Outputs go to ~/data/xenium_lung/processed/spatial_domains
Part A: Load the data¶
Purpose: start from the cells, labels and coordinates the core spatial analyses used, so the domains can be compared directly with the Step 7 niches.
Approach: load the Step 9 object (xu.load_step("time")): the high-confidence cells of the core workflow, with log-normalised expression set as .X, cell types and lineages, the Step 7 k-means niches, and centroids in µm. All domain labels are collected in one table, domains.
Expected output: cells x genes, the tissue extent in mm, and the number of Step 7 niches.
Note: Step 6 mixed cells are excluded, so domain composition describes high-confidence cells only.
adata = xu.load_step("time")
adata.X = adata.layers["lognorm"].copy()
xy = adata.obsm["spatial"]
adata.obs["x"], adata.obs["y"] = xy[:, 0], xy[:, 1]
domains = pd.DataFrame({"niche": adata.obs["niche"].astype(str)}, index=adata.obs_names) # all domain labels, collected
lineage = adata.obs["lineage"].astype(str)
print(f"{adata.n_obs:,} cells x {adata.n_vars} genes; tissue {np.ptp(xy[:, 0]) / 1000:.1f} x {np.ptp(xy[:, 1]) / 1000:.1f} mm")
print(f"Step 7 k-means niches: {domains['niche'].nunique()}")
139,129 cells x 377 genes; tissue 10.8 x 3.2 mm Step 7 k-means niches: 10
Interpretation. 139,129 cells x 377 genes over a 10.8 x 3.2 mm section, with the 10 Step 7 niches as the baseline.
Part B: BANKSY domains and their parameter stability¶
Purpose: find domains from each cell's own expression combined with its neighbours' expression, and check that the chosen map does not depend on a knife-edge parameter choice.
Approach: dm.banksy_grid builds the BANKSY matrix (15 nearest neighbours with Gaussian distance weights; no azimuthal Gabor filter, max_m = 0, as in the R workflow), then runs PCA (20 PCs) and Leiden for every combination of λ = 0.5 / 0.8 / 0.9 and resolution 0.1-0.3. Every setting is compared with the chosen one (λ = 0.8, resolution 0.2) by adjusted Rand index (ARI: 1 = identical partition, 0 = chance).
Expected output: the run time, a table of ARI with the chosen setting, and a table of the number of domains per setting. Look for high ARI around the chosen setting and how the partition changes with λ.
Note: about 1-2 min for the 12 settings. ARI near the chosen setting measures robustness, not correctness: a stable map can still be biologically uninformative.
LAMBDAS, RESOLUTIONS = [0.5, 0.8, 0.9], [0.1, 0.15, 0.2, 0.3]
CHOSEN = (0.8, 0.2)
t0 = time.time()
grid = dm.banksy_grid(adata, LAMBDAS, RESOLUTIONS) # (lambda, resolution) -> domain labels; pybanksy, 15 neighbours, max_m 0
print(f"BANKSY: {len(grid)} parameter settings in {time.time() - t0:.0f} s")
chosen = grid[CHOSEN]
stab, n_dom = dm.grid_stability(grid, CHOSEN, LAMBDAS, RESOLUTIONS)
display(stab.style.format("{:.2f}").background_gradient(cmap="Greens", vmin=0, vmax=1)
.set_caption(f"ARI with the chosen setting (lambda {CHOSEN[0]}, resolution {CHOSEN[1]})"))
display(n_dom.style.set_caption("Number of domains"))
domains["banksy"] = chosen
BANKSY: 12 parameter settings in 71 s
| res 0.1 | res 0.15 | res 0.2 | res 0.3 | |
|---|---|---|---|---|
| lambda 0.5 | 0.60 | 0.62 | 0.63 | 0.54 |
| lambda 0.8 | 0.65 | 0.81 | 1.00 | 0.82 |
| lambda 0.9 | 0.68 | 0.58 | 0.58 | 0.84 |
| res 0.1 | res 0.15 | res 0.2 | res 0.3 | |
|---|---|---|---|---|
| lambda 0.5 | 7 | 8 | 8 | 12 |
| lambda 0.8 | 8 | 8 | 10 | 12 |
| lambda 0.9 | 8 | 8 | 8 | 9 |
Interpretation. The chosen setting gives 10 domains. At λ = 0.8, the neighbouring resolutions agree with it at ARI 0.81-0.82, so the map is not a knife-edge result; resolution 0.1 merges domains (8 domains, ARI 0.65). λ = 0.9 is less consistent (ARI 0.58-0.84), and λ = 0.5, which gives each cell's own expression more weight, departs further (ARI 0.54-0.63), as expected from the parameter that controls the cell-versus-neighbourhood balance.
Part C: CellCharter domains with automatic choice of K¶
Purpose: find domains from multi-scale neighbourhood structure, and let the data choose the number of domains by stability rather than by eye.
Approach:
- Embedding (
dm.cellcharter_embedding): 30 PCs of the scaled expression, concatenated with their means over 3 layers of the Delaunay graph after removing long links (120 dimensions). The original method uses an scVI embedding; PCA is the lighter substitute for one sample without batch effects. - Number of domains (
dm.cellcharter_autok):ClusterAutoKfits Gaussian mixtures for K = 3-13, three runs each, and scores stability (mean Fowlkes-Mallows index between runs and neighbouring K) for K = 4-12. The most stable K is used.
Expected output: the embedding size, a stability curve with the chosen K marked, and the number of domains. Check whether the chosen K is a clear peak or one of several similar values.
Note: RECOMPUTE_AUTOK = False reuses the saved labels; recomputing takes ~9 min and overwrites cellcharter_autok.csv.gz and cellcharter_stability.csv. CellCharter's random_state does not seed PyTorch, so identical calls gave different domains until PyTorch was seeded globally (seed_everything). The fit runs on the CPU, which was faster than the Mac GPU for this model (~20 s vs ~45 s per fit).
RECOMPUTE_AUTOK = False
AUTOK_FILE = OUT_DIR / "cellcharter_autok.csv.gz"
dm.cellcharter_embedding(adata) # 30 PCs, aggregated over 3 Delaunay layers (long links removed) -> obsm['X_cellcharter']
print(f"CellCharter embedding: {adata.obsm['X_cellcharter'].shape[1]} dimensions (30 PCs x 4: the cell + 3 neighbour layers)")
if RECOMPUTE_AUTOK or not AUTOK_FILE.exists():
t0 = time.time()
labels, stability = dm.cellcharter_autok(adata, (4, 12), runs=3) # PyTorch seeded explicitly; CPU
pd.DataFrame({"cellcharter": labels}, index=adata.obs_names).to_csv(AUTOK_FILE)
stability.to_csv(OUT_DIR / "cellcharter_stability.csv", index_label="K")
print(f"AutoK fitted in {(time.time() - t0) / 60:.1f} min")
domains["cellcharter"] = pd.read_csv(AUTOK_FILE, index_col=0)["cellcharter"].astype(str).reindex(adata.obs_names)
stability = pd.read_csv(OUT_DIR / "cellcharter_stability.csv", index_col=0)["stability"]
best_k = int(stability.idxmax())
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(stability.index.astype(int), stability.to_numpy(), "o-", color="#2a78d6")
ax.axvline(best_k, color="black", ls="--", lw=1)
ax.set_xlabel("number of domains (K)")
ax.set_ylabel("stability (mean FMI)")
ax.set_title(f"CellCharter: most stable K = {best_k}")
plt.tight_layout()
plt.show()
print(f"CellCharter: {domains['cellcharter'].nunique()} domains")
INFO Creating graph using `None` transform and `1` libraries.
CellCharter embedding: 120 dimensions (30 PCs x 4: the cell + 3 neighbour layers)
CellCharter: 7 domains
Interpretation. Stability peaks at K = 7 (mean FMI 0.80), ahead of K = 6 (0.78) and K = 12 (0.77), and falls to its minimum at K = 9 (0.66). The peak is modest rather than sharp, so K = 7 is the most reproducible choice among close alternatives, not the only defensible one.
Part D: Compare the domain sets¶
D1. What each domain is made of¶
Purpose: give every domain an interpretable identity, so the three methods can be compared by content rather than by arbitrary cluster numbers.
Approach: dm.name_domains names each domain from its own lineage composition ("D: top lineage x% [+ second lineage y%]", the second shown if at least 15%) and numbers domains by decreasing epithelial share, the same rule as the Step 7 niches. A heatmap shows the lineage composition of every domain for each method.
Expected output: three heatmaps (k-means niches, BANKSY, CellCharter), one row per domain with its size. Identify the tumour, myeloid, lymphoid, stromal and endothelial domains in each method.
Note: names are generated from composition only; two domains with the same name can differ in expression or location, which D2 and D3 examine.
LINEAGES = list(xu.LINEAGE_COLORS)
def name_domains(labels):
"""'D<i>: <top lineage> x% [+ <second> y%]', numbered by decreasing epithelial share (dm.name_domains)."""
return dm.name_domains(labels, lineage, LINEAGES)
METHODS = {"niche": "k-means niches (Step 7)", "banksy": "BANKSY", "cellcharter": "CellCharter"}
for m in ["banksy", "cellcharter"]:
domains[m] = name_domains(domains[m])
fig, axes = plt.subplots(1, 3, figsize=(17, 5.5), gridspec_kw={"width_ratios": [1, 1, 1]})
for ax, (m, title) in zip(axes, METHODS.items()):
comp = pd.crosstab(domains[m], lineage, normalize="index").reindex(columns=LINEAGES).fillna(0) * 100
comp = comp.loc[sorted(comp.index, key=lambda s: int(s.split(":")[0][1:]))]
ax.imshow(comp.to_numpy(), cmap="Blues", vmin=0, vmax=100, aspect="auto")
ax.set_xticks(range(len(LINEAGES)), LINEAGES, rotation=45, ha="right")
ax.set_yticks(range(len(comp)), [f"{d} (n={n:,})" for d, n in domains[m].value_counts().reindex(comp.index).items()], fontsize=8)
for i in range(comp.shape[0]):
for j in range(comp.shape[1]):
if comp.iat[i, j] >= 5:
ax.text(j, i, f"{comp.iat[i, j]:.0f}", ha="center", va="center", fontsize=7,
color="white" if comp.iat[i, j] > 60 else "black")
ax.set_title(title)
plt.suptitle("Lineage composition of each domain (% of cells)")
plt.tight_layout()
plt.show()
Interpretation.
- Tumour: k-means splits the epithelium by how much stroma is mixed in (N1-N3: 94% to 61% epithelial); BANKSY has four almost pure epithelial domains (87-95%); CellCharter has two pure ones (100%) plus an epithelial-myeloid domain (D3) and a mixed epithelial border domain (D4, 64%).
- Myeloid: a myeloid-dominated domain appears in k-means (N4, 68% myeloid) and BANKSY (D5, 74%); in CellCharter the myeloid cells form the minority of the epithelial-myeloid domain D3 (31%).
- Lymphoid: k-means separates a T / NK-rich (N7) and a B / plasma-rich (N8) aggregate domain; BANKSY has one mixed lymphoid domain (D9, 43% T / NK, 42% B / plasma); CellCharter's lymphoid domain (D7) is larger and more diluted (29% T / NK, 28% B / plasma, 14,233 cells).
- Stroma and vessels: k-means and BANKSY each have an endothelial-majority domain (N10, D8); CellCharter places vessels in a stromal-endothelial domain (D6, 46% / 35%).
- One small artefact: BANKSY's D6 (141 cells, 70% stromal) lies almost entirely (96%) on the detached tissue sliver at the bottom of the section (see D2).
D2. Where the domains are, and how much the methods agree¶
Purpose: measure agreement between methods in a way that separates real disagreement about tissue structure from different ways of subdividing the tumour.
Approach: compute two measures for each pair of methods:
- ARI on domain labels: strict, and dominated by the tumour, which holds about 60% of cells.
- Compartment agreement: every domain is assigned its majority compartment (
dm.compartment: tumour = epithelial; immune = T / NK + B / plasma + myeloid + mast; stroma / vessels = stromal + endothelial), and the % of cells placed in the same compartment by both methods is reported.
The three domain maps are drawn on tissue coordinates.
Expected output: a pair-wise table of ARI and compartment agreement, and three tissue maps. Compare where the methods draw tumour subdivisions, immune structures and stroma.
Note: a domain is assigned to one compartment as a whole, so a mixed domain (for example epithelial 69% + myeloid 31%) counts all its cells as tumour.
COMPARTMENT = dm.COMPARTMENT # lineage -> tumour / immune / stroma + vessels
compartments = pd.DataFrame({m: dm.compartment(domains[m], lineage) for m in METHODS})
pairs = [("niche", "banksy"), ("niche", "cellcharter"), ("banksy", "cellcharter")]
agree = pd.DataFrame([{"pair": f"{METHODS[a]} vs {METHODS[b]}", "ARI (domains)": ari(domains[a], domains[b]),
"same compartment (% cells)": 100 * (compartments[a] == compartments[b]).mean()} for a, b in pairs]).set_index("pair")
display(agree.style.format({"ARI (domains)": "{:.2f}", "same compartment (% cells)": "{:.0f}"}))
fig, axes = plt.subplots(3, 1, figsize=(14, 13))
for ax, (m, title) in zip(axes, METHODS.items()):
labs = domains[m]
order = sorted(labs.unique(), key=lambda s: int(s.split(":")[0].lstrip("ND")))
cmap = plt.get_cmap("tab20" if len(order) > 10 else "tab10")
for i, d in enumerate(order):
sel = (labs == d).to_numpy()
ax.scatter(xy[sel, 0] / 1000, xy[sel, 1] / 1000, s=0.15, color=cmap(i), label=d, rasterized=True)
ax.set_aspect("equal")
ax.invert_yaxis()
ax.set_title(title)
ax.set_xlabel("x (mm)")
ax.set_ylabel("y (mm)")
ax.legend(markerscale=25, fontsize=7, loc="center left", bbox_to_anchor=(1.0, 0.5), frameon=False)
plt.tight_layout()
plt.show()
| ARI (domains) | same compartment (% cells) | |
|---|---|---|
| pair | ||
| k-means niches (Step 7) vs BANKSY | 0.17 | 86 |
| k-means niches (Step 7) vs CellCharter | 0.29 | 79 |
| BANKSY vs CellCharter | 0.26 | 85 |
Interpretation. ARI between methods is low (0.17-0.29), but compartment agreement is high (79-86% of cells), so the methods agree on the tissue layout and differ mainly in how they subdivide it. The maps show the difference: BANKSY divides the tumour into large contiguous regions (left, centre and right of the section), CellCharter into interleaved bands that follow nest structure, and the k-means niches into a fine mosaic. All three place stroma and vessels at the right-hand edge (x > 9.5 mm) and lymphoid structures near the TLS at ~(7.0, 2.0) mm. D3 tests what drives the tumour split.
D3. Why BANKSY and CellCharter disagree inside the tumour¶
Purpose: identify the axis along which each method divides the tumour, which explains the low ARI.
Approach: for each tumour (epithelial-named) domain, report the % of its tumour cells in each Step 6 tumour subtype, and the median depth in the nest: the distance from each epithelial cell to the nearest non-epithelial cell.
Expected output: one table per method. A domain dominated by one subtype points to an expression-programme split; domains with mixed subtypes but different depth point to an architecture split.
Note: depth is measured to the nearest non-epithelial high-confidence cell, so excluded mixed cells at nest edges slightly inflate depth.
from scipy.spatial import cKDTree
is_epi = (lineage == "Epithelial").to_numpy()
depth = pd.Series(cKDTree(xy[~is_epi]).query(xy[is_epi])[0], index=adata.obs_names[is_epi])
tumour = adata.obs_names[adata.obs["cell_type"].isin(xu.TUMOUR_TYPES).to_numpy()]
for m in ["banksy", "cellcharter"]:
epi_domains = [d for d in sorted(domains[m].unique(), key=lambda s: int(s.split(":")[0][1:])) if "Epithelial" in d.split(":")[1].split("%")[0]]
sub = pd.crosstab(domains.loc[tumour, m], adata.obs.loc[tumour, "cell_type"], normalize="index").reindex(epi_domains) * 100
sub.columns = [c.replace("Tumour epithelial ", "").replace("Proliferating tumour", "proliferating") for c in sub.columns]
sub["median depth in nest (um)"] = depth.groupby(domains.loc[depth.index, m]).median().reindex(epi_domains)
display(sub.style.format("{:.0f}").background_gradient(cmap="Blues", vmin=0, vmax=100, subset=sub.columns[:-1])
.set_caption(f"{METHODS[m]}: tumour domains, % of their tumour cells by subtype, and depth"))
| proliferating | (CYP2B6/CFTR) | (MALL/TCIM) | (MYC/CAPN8) | median depth in nest (um) | |
|---|---|---|---|---|---|
| banksy | |||||
| D1: Epithelial 95% | 2 | 78 | 16 | 4 | 22 |
| D2: Epithelial 93% | 8 | 7 | 12 | 73 | 23 |
| D3: Epithelial 90% | 3 | 21 | 67 | 9 | 18 |
| D4: Epithelial 87% | 2 | 66 | 22 | 10 | 18 |
| proliferating | (CYP2B6/CFTR) | (MALL/TCIM) | (MYC/CAPN8) | median depth in nest (um) | |
|---|---|---|---|---|---|
| cellcharter | |||||
| D1: Epithelial 100% | 3 | 33 | 23 | 41 | 30 |
| D2: Epithelial 100% | 7 | 30 | 38 | 25 | 16 |
| D3: Epithelial 69% + Myeloid 31% | 2 | 38 | 27 | 33 | 15 |
| D4: Epithelial 64% | 4 | 52 | 35 | 9 | 13 |
Interpretation.
- BANKSY follows tumour subtype. Each of its tumour domains is dominated (66-78%) by one tumour programme: D1 and D4 by CYP2B6/CFTR, D2 by MYC/CAPN8, D3 by MALL/TCIM. Their depths are similar (18-23 µm). The programmes occupy different regions of the section, which is why BANKSY's tumour domains are large and contiguous.
- CellCharter follows depth. Its two pure-tumour domains mix all subtypes but differ in depth: nest core (D1, median 30 µm) vs nest edge (D2, 16 µm); the border domains D3 and D4 are shallower still (13-15 µm).
- Why: BANKSY keeps a strong weight on each cell's own expression, so expression programmes drive the split; CellCharter aggregates three rings of neighbours, so architecture drives it. The k-means niches use lineage composition only, so they cannot see tumour subtype and split the tumour by stromal admixture.
- Practical rule: choose the method by the question: BANKSY to map regions with different tumour programmes, CellCharter to study architecture (core vs edge, borders), and report parameter stability either way.
Part E: Same cell type, different domain¶
E1. Design of the pseudobulk comparison¶
Purpose: test whether one annotated cell type carries different expression programmes in different domains, while guarding against the two main artefacts: pseudo-replication and spillover.
Approach:
- One cell type at a time. Macrophage types already sort by domain (most FCGR1A+ and alveolar macrophages sit in the CellCharter tumour-myeloid domain, most CD163+ macrophages in stroma), so pooling all macrophages would rediscover type differences. Two comparisons: FCGR1A+ macrophages in the tumour-myeloid domain vs stroma, and CD163+ macrophages in the lymphoid domain vs stroma. Domains are chosen by rule from the CellCharter compositions (
ROLE). - Spatial pseudo-replicates. Counts are summed per (1 x 1 mm block, domain) group with at least 10 cells; only blocks containing both domains are kept, and the model is
~ block + domain, so each block is compared with itself. - Spillover filter. A gene is flagged foreign if its mean expression in any other lineage is at least 4x its mean in myeloid cells (
dm.foreign_genes). The DE helper runs on 10x counts and on Proseg counts (notebook 02) for the same cells.
Expected output: the three domains used, and the number of genes flagged as foreign to myeloid cells.
Note: this is one section from one patient: blocks are pseudo-replicates, not independent samples, and neighbouring blocks are correlated, so p-values are optimistic. They are used to rank genes, with effect sizes read first.
pro = sc.read_h5ad(xu.output_dir("resegmentation") / f"{xu.SAMPLE}_proseg_counts.h5ad")
cc_labels = domains["cellcharter"]
comp_cc = pd.crosstab(cc_labels, lineage, normalize="index")
epi_major = comp_cc.index[comp_cc["Epithelial"] >= 0.5]
ROLE = {"tumour-myeloid": comp_cc.loc[epi_major, "Myeloid"].idxmax(), # epithelial-majority domain richest in myeloid cells
"stroma": cc_labels[compartments["cellcharter"] == "stroma / vessels"].value_counts().idxmax(), # largest stromal domain
"lymphoid": (comp_cc["T / NK"] + comp_cc["B / plasma"]).idxmax()}
print("Domains used:", *[f" {r}: {d}" for r, d in ROLE.items()], sep="\n")
block = dm.blocks(xy, 1000, adata.obs_names) # 1 x 1 mm spatial pseudo-replicates
# Foreign genes: mean expression in some non-myeloid lineage >= 4x the myeloid mean (10x counts, normalised)
foreign_of = dm.foreign_genes(adata, lineage)
print(f"\n{foreign_of.notna().sum()} of {adata.n_vars} genes flagged as foreign to myeloid cells")
def pseudobulk_de(cell_type, dom_a, dom_b, counts="10x", min_cells=10):
"""Paired pseudobulk DESeq2 (dm.pseudobulk_de): cell_type in domain dom_a vs dom_b, ~ block + domain."""
cells = adata.obs_names[(adata.obs["cell_type"] == cell_type).to_numpy() & cc_labels.isin([dom_a, dom_b]).to_numpy()]
if counts == "proseg":
cells = cells.intersection(pro.obs_names)
X, genes = pro[cells].X, pro.var_names
else:
X, genes = adata[cells].layers["counts"], adata.var_names
return dm.pseudobulk_de(X, genes, cells, block, cc_labels, dom_a, dom_b, min_cells)
Domains used: tumour-myeloid: D3: Epithelial 69% + Myeloid 31% stroma: D5: Stromal 45% + Endothelial 19% lymphoid: D7: T / NK 29% + B / plasma 28% 195 of 377 genes flagged as foreign to myeloid cells
Interpretation. The rules select D3 (epithelial 69% + myeloid 31%) as the tumour-myeloid domain, D5 (stromal 45% + endothelial 19%) as stroma, and D7 (T / NK 29% + B / plasma 28%) as lymphoid. 77% of FCGR1A+ macrophages and 90% of alveolar macrophages lie in D3, and 62% of CD163+ macrophages in D5. 195 of 377 genes are flagged as foreign to myeloid cells, so about half the panel cannot be interpreted as a macrophage-intrinsic change without further evidence.
E2. Differential expression results¶
Purpose: identify genes that change with domain within one macrophage type, and separate macrophage-intrinsic changes from spillover.
Approach: for each comparison, run the paired pseudobulk DESeq2 test (dm.pseudobulk_de) on 10x counts and on Proseg counts, then classify each gene (dm.call_genes, adjusted p < 0.05): n.s.; spillover (significant but foreign); not confirmed with Proseg (significant, not foreign, but not significant in the same direction with Proseg counts); or macrophage-intrinsic (passes all three).
Expected output: per comparison, the paired blocks and cells used, the correlation of log2 fold changes between 10x and Proseg counts, gene counts per call, and the macrophage-intrinsic genes with their 10x and Proseg statistics (positive = higher in the first domain).
Note: genes with too few counts are dropped before testing, and genes with undefined adjusted p-values count as n.s.; the foreign-gene filter only removes genes of other lineages, so spillover from neighbouring myeloid cells is not excluded.
COMPARISONS = [("FCGR1A+ macrophage", "tumour-myeloid", "stroma"), ("CD163+ macrophage", "lymphoid", "stroma")]
de = {}
for cell_type, a_role, b_role in COMPARISONS:
r10, info = pseudobulk_de(cell_type, ROLE[a_role], ROLE[b_role])
rP, _ = pseudobulk_de(cell_type, ROLE[a_role], ROLE[b_role], counts="proseg")
# macrophage-intrinsic = significant, not a foreign-lineage gene, and significant the same way with Proseg counts
r = dm.call_genes(r10, rP, foreign_of)
de[cell_type] = r
print(f"\n{cell_type}: {a_role} vs {b_role} | {info['blocks']} paired blocks, {info['cells']:,} cells | "
f"log2FC agreement 10x vs Proseg: r = {r[['log2FoldChange', 'log2FoldChange_proseg']].corr().iloc[0, 1]:.2f}")
display(r["call"].value_counts().to_frame("genes").T)
top = r[r["call"] == "macrophage-intrinsic"].sort_values("log2FoldChange")
show = pd.concat([top.head(8), top.tail(8)]).drop_duplicates()
display(show[["baseMean", "log2FoldChange", "padj", "log2FoldChange_proseg", "padj_proseg"]]
.style.format({"baseMean": "{:.0f}", "log2FoldChange": "{:+.2f}", "log2FoldChange_proseg": "{:+.2f}",
"padj": "{:.1e}", "padj_proseg": "{:.1e}"})
.set_caption(f"Macrophage-intrinsic genes (positive = higher in {a_role})"))
FCGR1A+ macrophage: tumour-myeloid vs stroma | 13 paired blocks, 1,501 cells | log2FC agreement 10x vs Proseg: r = 0.88
| call | n.s. | spillover (foreign gene) | macrophage-intrinsic | not confirmed with Proseg |
|---|---|---|---|---|
| genes | 238 | 42 | 15 | 3 |
| baseMean | log2FoldChange | padj | log2FoldChange_proseg | padj_proseg | |
|---|---|---|---|---|---|
| ALDH1A3 | 1 | -2.05 | 2.8e-03 | -3.13 | 3.6e-02 |
| CTSK | 18 | -1.38 | 2.0e-04 | -1.54 | 5.2e-04 |
| TREM2 | 38 | -0.90 | 2.0e-05 | -0.94 | 1.1e-07 |
| PLA2G7 | 76 | -0.89 | 5.3e-04 | -0.73 | 5.2e-04 |
| LY86 | 11 | -0.85 | 8.5e-04 | -0.74 | 1.3e-02 |
| LPL | 14 | -0.74 | 3.0e-03 | -0.87 | 3.9e-03 |
| IGSF6 | 79 | -0.53 | 1.4e-02 | -0.53 | 8.6e-04 |
| AIF1 | 75 | +0.31 | 3.3e-02 | +0.43 | 6.6e-04 |
| CD163 | 63 | +0.36 | 3.1e-02 | +0.50 | 5.1e-04 |
| BCL2L11 | 12 | +0.70 | 1.6e-02 | +0.74 | 3.1e-02 |
| AQP9 | 5 | +1.61 | 5.9e-03 | +1.95 | 5.2e-04 |
| VCAN | 29 | +1.65 | 1.1e-04 | +1.68 | 2.0e-06 |
| MCEMP1 | 7 | +1.88 | 8.8e-05 | +1.58 | 5.2e-04 |
| MARCO | 11 | +2.20 | 1.2e-06 | +2.49 | 1.3e-08 |
| CXCL2 | 8 | +2.89 | 3.4e-04 | +2.73 | 4.4e-04 |
CD163+ macrophage: lymphoid vs stroma | 14 paired blocks, 1,705 cells | log2FC agreement 10x vs Proseg: r = 0.75
| call | n.s. | spillover (foreign gene) | macrophage-intrinsic | not confirmed with Proseg |
|---|---|---|---|---|
| genes | 255 | 10 | 9 | 6 |
| baseMean | log2FoldChange | padj | log2FoldChange_proseg | padj_proseg | |
|---|---|---|---|---|---|
| VCAN | 6 | -1.42 | 9.3e-04 | -1.02 | 2.7e-03 |
| MNDA | 9 | -1.29 | 3.0e-04 | -0.74 | 4.0e-02 |
| FCGR3A | 20 | -1.02 | 9.2e-05 | -1.14 | 9.3e-05 |
| MS4A4A | 94 | +0.29 | 4.8e-02 | +0.33 | 1.5e-03 |
| MPEG1 | 50 | +0.42 | 2.3e-02 | +0.46 | 1.5e-03 |
| KCNMA1 | 26 | +0.43 | 5.0e-02 | +0.41 | 4.9e-02 |
| CXCR4 | 28 | +0.44 | 3.9e-02 | +0.44 | 2.7e-02 |
| PLA2G7 | 37 | +0.91 | 6.5e-04 | +0.90 | 4.1e-03 |
| LYVE1 | 3 | +1.76 | 1.0e-03 | +2.21 | 9.3e-05 |
Interpretation. Both comparisons rest on 13-14 paired blocks and about 1,500-1,700 cells. Fold changes agree well between segmentations (r = 0.88 and 0.75). In the FCGR1A+ comparison, 42 of the 60 significant genes are spillover, 15 are macrophage-intrinsic and 3 are not confirmed with Proseg; in the CD163+ comparison, 10 of 25 are spillover and 9 are intrinsic. The volcano plots below show where these groups fall.
E3. Volcano plots by call¶
Purpose: show how the filters change the picture a naive differential expression analysis would give.
Approach: plot log2 fold change against −log10 adjusted p for every tested gene, coloured by call, and label the six most significant macrophage-intrinsic genes and the four most significant spillover genes.
Expected output: one volcano plot per comparison; the legend gives the number of genes per call. Compare which genes would rank highest without the spillover filter.
Note: adjusted p-values are clipped at 1e-30 for display.
CALL_COLORS = {"macrophage-intrinsic": "#c0392b", "spillover (foreign gene)": "#95a5a6",
"not confirmed with Proseg": "#f39c12", "n.s.": "#d5d8dc"}
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
for ax, (cell_type, a_role, b_role) in zip(axes, COMPARISONS):
r = de[cell_type].dropna(subset=["padj"])
for call, col in CALL_COLORS.items():
s = r[r["call"] == call]
ax.scatter(s["log2FoldChange"], -np.log10(s["padj"].clip(lower=1e-30)), s=14, color=col, label=f"{call} ({len(s)})")
lab = pd.concat([r[r["call"] == "macrophage-intrinsic"].nsmallest(6, "padj"),
r[r["call"] == "spillover (foreign gene)"].nsmallest(4, "padj")])
for g, row in lab.iterrows():
ax.annotate(g, (row["log2FoldChange"], -np.log10(max(row["padj"], 1e-30))), fontsize=7, xytext=(3, 2), textcoords="offset points")
ax.axhline(-np.log10(0.05), color="black", lw=0.6, ls=":")
ax.axvline(0, color="black", lw=0.6)
ax.set_xlabel(f"log2 fold change ({a_role} / {b_role})")
ax.set_ylabel("-log10 adjusted p")
ax.set_title(f"{cell_type}: {a_role} vs {b_role}")
ax.legend(fontsize=7, frameon=False, loc="upper center", bbox_to_anchor=(0.5, -0.13), ncol=2)
plt.tight_layout()
plt.show()
Interpretation.
FCGR1A+ macrophages, tumour-myeloid domain vs stroma.
- Higher in the tumour-myeloid domain: MARCO and MCEMP1 (alveolar-macrophage genes), CXCL2 (a neutrophil-attracting chemokine), VCAN and AQP9 (monocyte-associated, inflammatory genes). Alveolar macrophages are the most common myeloid cell in this domain (42%).
- Higher in stroma: TREM2, CTSK, LPL and PLA2G7, genes of the lipid-handling, matrix-remodelling programme described for TREM2+ tumour-associated macrophages.
- Spillover dominates the naive result: 8 of the 10 most significant genes are foreign, mostly fibroblast and smooth-muscle genes in stromal macrophages (THBS2, FBN1, ACTA2, LTBP2, MYLK) and epithelial genes in tumour-zone macrophages (KRT7, CYP2B6, TCIM).
- Caveat: MARCO and MCEMP1 are myeloid genes, so spillover from the many neighbouring alveolar macrophages cannot be excluded by the foreign-gene filter. Proseg confirming them (log2 fold change +2.49 and +1.58) makes this less likely, but not impossible.
CD163+ macrophages, lymphoid domain vs stroma. A smaller shift: LYVE1, PLA2G7, CXCR4, MPEG1 and MS4A4A higher in the lymphoid domain; FCGR3A, MNDA and VCAN (monocyte-associated genes) lower. This pattern is consistent with more tissue-resident (LYVE1+), less monocyte-like macrophages in lymphoid aggregates. The spillover genes here include T- and B-cell genes (CD3D, CD79A, SLAMF7) higher in the lymphoid domain, as expected from the neighbourhood.
What this shows: the same annotated macrophage type carries different programmes depending on its domain, and much of the naive differential expression in imaging data is spillover from the neighbourhood. Filtering foreign genes and confirming with a second segmentation separates the two.
Part F: Export for Xenium Explorer¶
Purpose: make the cell types, niches and domains browsable on the tissue, over the morphology and H&E images, by anyone with the Xenium bundle.
Approach: write each labelling as a cell groups CSV (cell_id, group, color; dm.explorer_groups), the format Xenium Explorer imports under Cell > Add cell categorization. Lineages use the project palette; other labellings use tab20. All domain labels and both differential expression tables are also saved.
Expected output: the number of cell-group files and the paths of the saved tables.
Note: overwrites the files in xenium_explorer/, domains.csv.gz and the pseudobulk_de_*.csv tables. Only high-confidence cells are included, so mixed cells appear unlabelled in Xenium Explorer.
EXPLORER_DIR = OUT_DIR / "xenium_explorer"
EXPLORER_DIR.mkdir(exist_ok=True)
groupings = {"cell_types": adata.obs["cell_type"].astype(str), "lineages": lineage, "niches_kmeans": domains["niche"],
"domains_banksy": domains["banksy"], "domains_cellcharter": domains["cellcharter"]}
for name, labels in groupings.items(): # Xenium Explorer cell groups: cell_id, group, color
dm.explorer_groups(labels, adata.obs_names, xu.LINEAGE_COLORS if name == "lineages" else None) \
.to_csv(EXPLORER_DIR / f"{name}.csv", index=False)
domains.to_csv(OUT_DIR / "domains.csv.gz", index_label="cell_id")
print(f"Saved {len(groupings)} cell-group files to {xu.display_path(EXPLORER_DIR)}")
print(f"Saved all domain labels to {xu.display_path(OUT_DIR / 'domains.csv.gz')}")
for name, r in de.items():
r.to_csv(OUT_DIR / f"pseudobulk_de_{name.split('+')[0].lower()}.csv")
Saved 5 cell-group files to ~/data/xenium_lung/processed/spatial_domains/xenium_explorer Saved all domain labels to ~/data/xenium_lung/processed/spatial_domains/domains.csv.gz
Interpretation. Five labellings are exported (cell types, lineages, k-means niches, BANKSY and CellCharter domains), together with the domain table used by the project's interactive viewer.
Conclusions¶
- Domains are reproducible at the compartment level. k-means niches, BANKSY and CellCharter agree on the tumour / immune / stroma compartment of 79-86% of cells, and all three recover lymphoid aggregates, a myeloid-rich tumour zone and stroma with vessels.
- The methods split the tumour along different axes. BANKSY by tumour programme (regions dominated by one subtype), CellCharter by architecture (nest core vs edge). Low ARI between methods (0.17-0.29) reflects that difference, not failure; the choice of method should follow the question.
- Stability was measured, not assumed. BANKSY's chosen map agrees with neighbouring resolutions at ARI ~0.8 (λ = 0.8); CellCharter chose K = 7 as its most stable number of domains, a modest peak, and needed explicit PyTorch seeding to be reproducible.
- Macrophage state depends on location. FCGR1A+ macrophages express an alveolar / inflammatory programme (MARCO, MCEMP1, CXCL2) in the tumour-myeloid domain and a TREM2 / CTSK / LPL programme in stroma; CD163+ macrophages in lymphoid aggregates are LYVE1-high and monocyte-gene-low.
- Spillover is the dominant artefact in spatial differential expression. 42 of 60 significant genes in the FCGR1A+ comparison were neighbouring-lineage genes. Results are reported only when they survive a foreign-gene filter and the Proseg counts.
- Limits: one section, so blocks are pseudo-replicates and p-values are optimistic; within-lineage spillover (for example from alveolar macrophages) is not excluded by the filter; domains are computed on high-confidence cells only.
The cell-group files in xenium_explorer/ load into Xenium Explorer for interactive browsing of all labellings.