02. Cell re-segmentation with Proseg¶
Author: MD Shakhaowat Hossain, PhD Candidate, Tulane University School of Medicine · Project website · Code
Why: in the core workflow (Step 6), about 10% of cells were flagged mixed: their 10x cell outlines captured transcripts from neighbouring cells (segmentation spillover). The problem was not random. CD4 T cells never formed a high-confidence population and dropped out of the spatial statistics, and tumour-contact CD8 T cells carried epithelial transcripts (the Step 9b control row). Any finding that depends on single-cell expression at a tissue interface is therefore open to the question: biology, or segmentation?
What: re-segment the tissue from the transcripts themselves with Proseg (Jones et al., Nature Methods 2025), a probabilistic method that starts from the 10x nuclei and reassigns each transcript to the cell that most plausibly produced it.
How the comparison stays fair: every Proseg cell keeps its 10x cell id, so the same cells can be followed before and after re-segmentation with the same Step 6 labels, and the core analyses can be repeated with identical code. Differences then come from segmentation alone.
| Part | Content |
|---|---|
| A | Baseline: the metrics a new segmentation must improve |
| B | Run Proseg (or reuse the saved run) |
| C | Compare with the baseline: whole slide, paired cells, an example window |
| D | Export the Proseg counts in the core-workflow layout |
| E | Re-test the core findings on Proseg cells: what holds, what changes, what is newly visible |
Setup: imports, paths and output folder¶
Purpose: load the shared ist_analysis package and define where this notebook writes, so paths and constants match the core notebook and the Snakemake pipeline.
Approach: import annotate, cluster, proseg and utils (as xu) from ist_analysis (repository root; pip install -e .), check the kernel with xu.check_kernel(), and create the output folder with xu.output_dir("resegmentation").
Expected output: the output folder, printed as a ~/... path.
Note: this notebook reads the Step 6, 7 and 9 objects from the core notebook (annotated, clustered, spatial, time); run the core notebook first.
from pathlib import Path
from ist_analysis import annotate, cluster, proseg
from ist_analysis import utils as xu # paths and constants shared with the core workflow and the pipeline
xu.check_kernel()
import warnings
import numpy as np
import pandas as pd
import pyarrow.compute as pc
import pyarrow.parquet as pq
from IPython.display import display
warnings.filterwarnings("ignore", category=FutureWarning)
OUT_DIR = xu.output_dir("resegmentation")
print(f"Outputs go to {xu.display_path(OUT_DIR)}")
Outputs go to ~/data/xenium_lung/processed/resegmentation
Part A: Baseline quality of the 10x segmentation¶
Purpose: fix the numbers a re-segmentation must improve before running it, so success is judged against criteria set in advance rather than chosen after seeing the result.
Approach: four measures on the 10x segmentation:
- Transcripts assigned to cells: gene transcripts with
qv >= 20whosecell_idis notUNASSIGNED, fromtranscripts.parquet. More is better only while cells stay plausible. - Cells flagged
mixedin Step 6 (annotation_confidence): a direct measure of spillover. - Lineage-impossible co-expression: % of T / NK cells with any EPCAM or MALL transcript, and % of tumour cells with any CD3E or TRAC transcript. In a clean segmentation both approach background.
- Recovered populations: the number of high-confidence
CD4 Tcells.
Expected output: a one-column table (baseline_10x_segmentation.csv); inspect the spillover rows and the CD4 T count.
Note: reads three columns of the full transcripts.parquet (about 12 million rows) into memory; the table is freed after counting. Measures 2-4 use the Step 6 labels, which were themselves made on 10x counts.
# 1. Transcript assignment (gene transcripts, qv >= 20)
tx = pq.read_table(xu.OUTS_DIR / "transcripts.parquet", columns=["feature_name", "qv", "cell_id"])
is_gene = ~pc.match_substring_regex(tx["feature_name"], "^(NegControl|Unassigned|Deprecated)").to_numpy(zero_copy_only=False)
q20 = (tx["qv"].to_numpy() >= 20) & is_gene
assigned = pc.not_equal(tx["cell_id"], "UNASSIGNED").to_numpy(zero_copy_only=False)
frac_assigned = (q20 & assigned).sum() / q20.sum()
del tx
# 2-3. Spillover measures from the Step 6 annotation (all QC-passing cells)
ann = xu.load_step("annotated")
counts = ann.layers["counts"].tocsc()
def positive(genes):
idx = [ann.var_names.get_loc(g) for g in genes]
return np.asarray((counts[:, idx] > 0).sum(axis=1)).ravel() > 0
lineage = ann.obs["lineage"].astype(str).to_numpy()
is_t = lineage == "T / NK"
is_tumour = ann.obs["cell_type"].isin(xu.TUMOUR_TYPES).to_numpy()
baseline = pd.Series({
"gene transcripts (qv>=20) assigned to a cell": f"{frac_assigned:.1%}",
"cells after QC": f"{ann.n_obs:,}",
"cells flagged 'mixed' (Step 6)": f"{(ann.obs['annotation_confidence'] == 'mixed').mean():.1%}",
"T / NK cells with EPCAM or MALL": f"{positive(['EPCAM', 'MALL'])[is_t].mean():.1%}",
"tumour cells with CD3E or TRAC": f"{positive(['CD3E', 'TRAC'])[is_tumour].mean():.1%}",
"CD4 T cells, high confidence": f"{((ann.obs['cell_type'] == 'CD4 T') & (ann.obs['annotation_confidence'] == 'high')).sum():,}",
}, name="10x segmentation (baseline)")
display(baseline.to_frame())
baseline.to_csv(OUT_DIR / "baseline_10x_segmentation.csv")
| 10x segmentation (baseline) | |
|---|---|
| gene transcripts (qv>=20) assigned to a cell | 86.6% |
| cells after QC | 154,394 |
| cells flagged 'mixed' (Step 6) | 9.9% |
| T / NK cells with EPCAM or MALL | 18.8% |
| tumour cells with CD3E or TRAC | 8.2% |
| CD4 T cells, high confidence | 0 |
Interpretation. The 10x segmentation assigns 86.6% of high-quality gene transcripts to cells, so coverage is not the problem; purity is. Nearly one T / NK cell in five (18.8%) contains epithelial transcripts, 8.2% of tumour cells contain T-cell transcripts, and no CD4 T cell reaches high confidence. These are the targets for Parts C and E.
Part B: Run Proseg¶
Purpose: produce a transcript-driven segmentation of the whole slide, the input for every comparison below.
Approach: Proseg is a command-line tool written in Rust (cargo install proseg). Its --xenium preset reads transcripts.parquet directly, keeps gene transcripts with qv >= 20, initialises cells from the 10x nuclei, and iteratively reassigns transcripts while cell shapes grow or shrink to fit them. The cost was profiled on a 1 x 1 mm crop and projected linearly before the full run. Runtime and peak memory are read back from run.log.
Expected output: the command (paths shortened) and a run summary: transcripts used, starting cells, wall time, CPU time and peak memory. Check that the measured cost agrees with the projection.
Note: RUN_PROSEG = False reuses the saved run. Setting it to True (or deleting the output) re-runs Proseg: about 30 min on all CPU cores and up to ~13 GB RAM, and --overwrite replaces the previous output in resegmentation/proseg/.
import re
import subprocess
RUN_PROSEG = False # True = (re)run Proseg (~30 min, all CPU cores, up to ~13 GB RAM); False = reuse the saved run
PROSEG = Path("~/.cargo/bin/proseg").expanduser()
PROSEG_DIR = OUT_DIR / "proseg"
PROSEG_DIR.mkdir(exist_ok=True)
args = proseg.proseg_args(xu.OUTS_DIR / "transcripts.parquet")
print("Command (paths shortened):")
print(" proseg " + " ".join(xu.display_path(Path(a)) if a.startswith("/") else a for a in args), "\n")
if RUN_PROSEG or not (PROSEG_DIR / "proseg-output.zarr").exists():
# caffeinate keeps macOS awake; /usr/bin/time -l records runtime and peak memory in run.log
proseg.run_proseg(PROSEG, xu.OUTS_DIR / "transcripts.parquet", PROSEG_DIR)
log = (PROSEG_DIR / "run.log").read_text()
grab = lambda pattern: float(re.search(pattern, log).group(1))
run = pd.Series({
"transcripts used (qv >= 20, genes)": f"{int(grab(r'(\d+) transcripts')):,}",
"starting cells (10x nuclei)": f"{int(grab(r'(\d+) cells')):,}",
"wall time": f"{grab(r'([\d.]+) real') / 60:.1f} min",
"CPU time": f"{grab(r'([\d.]+) user') / 60:.0f} min (~{grab(r'([\d.]+) user') / grab(r'([\d.]+) real'):.1f} cores busy)",
"peak memory": f"{grab(r'(\d+)\s+peak memory footprint') / 1e9:.1f} GB",
}, name="Proseg 3.2.0, full slide")
display(run.to_frame())
print("Profiled first on a 1 x 1 mm crop (596k transcripts: 4.3 min, 0.8 GB); the linear projection for the\n"
"full slide (~220 CPU-min, <= 14 GB) matched the measured 221 CPU-min and 12.5 GB.")
Command (paths shortened): proseg --xenium ~/data/xenium_lung/Xenium_V1_humanLung_Cancer_FFPE/outs/transcripts.parquet --output-counts counts.csv.gz --output-cell-metadata cell-metadata.csv.gz --output-cell-polygons cell-polygons.geojson.gz --exclude-spatialdata-transcripts --output-spatialdata proseg-output.zarr --overwrite
| Proseg 3.2.0, full slide | |
|---|---|
| transcripts used (qv >= 20, genes) | 10,550,258 |
| starting cells (10x nuclei) | 155,468 |
| wall time | 30.2 min |
| CPU time | 221 min (~7.3 cores busy) |
| peak memory | 12.5 GB |
Profiled first on a 1 x 1 mm crop (596k transcripts: 4.3 min, 0.8 GB); the linear projection for the full slide (~220 CPU-min, <= 14 GB) matched the measured 221 CPU-min and 12.5 GB.
Interpretation. Proseg used 10.55 million transcripts and started from 155,468 nuclei. The full run took 30.2 min with about 7 cores busy and peaked at 12.5 GB, matching the projection from the crop (~220 CPU-min, at most 14 GB). The crop test is worth repeating for any new dataset or machine: memory, not time, is the limit that would stop a laptop run.
Part C: Compare with the baseline¶
C1. Whole-slide summary¶
Purpose: check that re-segmentation keeps the data intact (cells, transcripts, depth) before asking whether it is cleaner.
Approach: read the Proseg SpatialData output (proseg.read_proseg) and the 10x matrix (cell_feature_matrix.h5), then summarise each with proseg.summarise: cells, cells with at least 10 transcripts, transcripts in cells, median transcripts and genes per cell, and median 2D cell area.
Expected output: a side-by-side table with a Proseg / 10x ratio column. A large drop in transcripts or transcripts per cell would mean signal was lost; a drop in genes per cell at similar depth would mean foreign genes were removed.
Note: Proseg areas are 2D outline areas from its polygons; Proseg models cells in 3D, so area is only a rough comparison.
import anndata as ad
import scipy.sparse as sp
import scanpy as sc
import spatialdata as sd
# Proseg output (SpatialData: count table + cell polygons; 2D outline areas in um^2) and the 10x matrix
pro, pro_area, pro_sd = proseg.read_proseg(PROSEG_DIR / "proseg-output.zarr")
tenx = sc.read_10x_h5(xu.OUTS_DIR / "cell_feature_matrix.h5") # gene features only
cells10 = pd.read_parquet(xu.OUTS_DIR / "cells.parquet").set_index("cell_id")
summarise = proseg.summarise # cells, transcripts, medians per segmentation
compare = pd.DataFrame({
"10x segmentation": summarise(tenx, cells10.loc[tenx.obs_names, "cell_area"].to_numpy()),
"Proseg": summarise(pro, pro_area)})
compare["Proseg / 10x"] = (compare["Proseg"] / compare["10x segmentation"]).map("{:.2f}".format)
display(compare.style.format({"10x segmentation": "{:,.0f}", "Proseg": "{:,.0f}"}))
print("Proseg starts from the 10x nuclei, so cells without a nucleus (3.5%) are not re-created; it keeps 94% of the\n"
"transcripts in fewer, slightly larger cells, and every Proseg cell keeps its 10x id (original_cell_id).")
| 10x segmentation | Proseg | Proseg / 10x | |
|---|---|---|---|
| cells | 162,254 | 155,468 | 0.96 |
| cells with >= 10 transcripts | 154,472 | 144,445 | 0.94 |
| transcripts in cells | 9,182,510 | 8,643,227 | 0.94 |
| median transcripts / cell | 46 | 45 | 0.98 |
| median genes / cell | 29 | 24 | 0.83 |
| median cell area (um^2) | 61 | 72 | 1.19 |
Proseg starts from the 10x nuclei, so cells without a nucleus (3.5%) are not re-created; it keeps 94% of the transcripts in fewer, slightly larger cells, and every Proseg cell keeps its 10x id (original_cell_id).
Interpretation. Proseg keeps 94% of the transcripts held in 10x cells (8.64 vs 9.18 million) in 4% fewer cells (155,468 vs 162,254), because the 3.5% of 10x cells without a nucleus are not re-created. Depth is essentially unchanged (median 45 vs 46 transcripts per cell) while median genes per cell fall from 29 to 24, and cells are slightly larger (72 vs 61 µm²). Fewer genes at the same depth is the pattern expected if foreign transcripts were removed; C2 tests that directly.
C2. Paired before / after on the same cells¶
Purpose: measure spillover on the same cells before and after re-segmentation, so differences cannot come from a change in which cells are compared.
Approach: match cells by 10x id (original_cell_id) and group them by their Step 6 lineage. For each lineage, proseg.paired_spillover reports the % of cells containing any transcript of a lineage they should not express (for example EPCAM / MALL in T / NK cells). proseg.mixed_spillover repeats the test for the six largest Step 6 mixed groups, using each group's contaminant genes.
Expected output: two tables with 10x and Proseg % positive per row. Inspect which contaminations fall sharply, which barely change, and whether the mixed groups are cleaned or only reduced.
Note: "any transcript" is a strict test: a single misassigned transcript makes a cell positive, so residual rates include genuine low-level expression and background as well as spillover.
# ---- Paired before / after: the same cells, labelled by the core workflow (Step 6) -----------
ann = xu.load_step("annotated", backed="r")
labels = ann.obs[["cell_type", "lineage", "annotation_confidence"]].copy()
pro.obs_names = pro.obs["original_cell_id"].astype(str).to_numpy()
paired = labels.index.intersection(pro.obs_names).intersection(tenx.obs_names)
print(f"{len(paired):,} cells present in both segmentations and annotated in Step 6\n")
positive = proseg.positive # per cell: any transcript of the given genes
SPILL_TESTS = proseg.SPILL_TESTS # (Step 6 lineage, foreign genes it should NOT express)
spill = proseg.paired_spillover(tenx, pro, labels, paired)
display(spill.style.format({"n": "{:,}", "10x: % positive": "{:.1f}", "Proseg: % positive": "{:.1f}", "reduction": "{:.1f}x"}))
# Cells flagged 'mixed' in Step 6: does Proseg clean them up? (their contaminant genes, before / after)
mix = proseg.mixed_spillover(tenx, pro, labels, paired)
display(mix.style.format({"n": "{:,}", "10x: % positive": "{:.1f}", "Proseg: % positive": "{:.1f}"}))
print("Each row follows the SAME cells before and after re-segmentation. Falling contaminant rates mean Proseg\n"
"is reassigning neighbours' transcripts to the right cells, the problem behind the Step 6 'mixed' labels.")
148,992 cells present in both segmentations and annotated in Step 6
| cells (Step 6 lineage) | foreign markers | n | 10x: % positive | Proseg: % positive | reduction | |
|---|---|---|---|---|---|---|
| 0 | T / NK | EPCAM / MALL | 14,227 | 18.8 | 1.5 | 12.7x |
| 1 | T / NK | CD163 / CD68 | 14,227 | 16.4 | 15.2 | 1.1x |
| 2 | Myeloid | EPCAM / MALL | 17,809 | 46.5 | 23.2 | 2.0x |
| 3 | Myeloid | CD3E / TRAC | 17,809 | 26.7 | 20.0 | 1.3x |
| 4 | Epithelial | CD3E / TRAC | 81,008 | 8.5 | 4.1 | 2.1x |
| 5 | Stromal | EPCAM / MALL | 18,228 | 30.1 | 4.2 | 7.2x |
| Step 6 'mixed' group | n | contaminant genes | 10x: % positive | Proseg: % positive | |
|---|---|---|---|---|---|
| 0 | Myeloid / epithelial mixed | 3,086 | EPCAM / MALL | 95.0 | 84.8 |
| 1 | T / myeloid mixed | 2,922 | CD163 / MS4A6A | 65.0 | 51.8 |
| 2 | CD4 T | 2,570 | VWF / CAVIN1 | 40.2 | 20.9 |
| 3 | Myeloid / fibroblast mixed | 2,090 | FBN1 / PDGFRA | 53.2 | 42.8 |
| 4 | T / fibroblast mixed | 1,820 | FBN1 / PDGFRA | 51.4 | 42.9 |
| 5 | Myeloid / T mixed | 1,658 | CD3E / TRAC | 78.7 | 66.3 |
Each row follows the SAME cells before and after re-segmentation. Falling contaminant rates mean Proseg is reassigning neighbours' transcripts to the right cells, the problem behind the Step 6 'mixed' labels.
Interpretation.
- Tumour-to-neighbour spillover is largely removed. T / NK cells with epithelial transcripts drop from 18.8% to 1.5% (12.7x), and stromal cells with epithelial transcripts from 30.1% to 4.2% (7.2x).
- Other contaminations improve less. Epithelial transcripts in myeloid cells halve (46.5% to 23.2%), T-cell transcripts in epithelial cells halve (8.5% to 4.1%), and T-cell transcripts in myeloid cells fall only from 26.7% to 20.0%.
- T / NK cells with macrophage transcripts barely change (16.4% to 15.2%). This may reflect tight T cell-macrophage contacts in immune hubs or 3D overlap within the 5 µm section, which a 2D reassignment cannot separate; it is not resolved here.
- The Step 6
mixedgroups are reduced, not cleaned. Endothelial transcripts in the CD4 T group halve (40.2% to 20.9%); the other groups fall by 8-13 points and most cells keep their contaminant genes.
C3. A boundary window, before and after¶
Purpose: show what the paired numbers mean at the level of individual outlines, at the place where spillover is most likely.
Approach: choose the 120 µm tile with the highest minimum of epithelial (EPCAM, MALL) and T-cell (CD3E, TRAC) transcripts, i.e. the densest tumour-T cell contact. Draw 10x and Proseg outlines (yellow) and those transcripts on DAPI from morphology_focus_0000.ome.tif, reading only the tiles inside the window.
Expected output: the window coordinates, then two panels. Compare whether outlines at the tumour-T cell boundary contain transcripts of the other lineage.
Note: one window chosen for maximal contact is an illustration, not evidence; C2 is the quantitative test. Proseg outlines are drawn on a 1 µm voxel grid, so their edges are stepped and harder to read by eye than 10x outlines.
import matplotlib.pyplot as plt
import tifffile
import zarr
from matplotlib.collections import PolyCollection
# ---- Look at it: outlines and transcripts on DAPI, before vs after ---------------------------
# Choose the window from the data: the 120 um tile where tumour and T cells meet most densely
# (highest min(epithelial transcripts, T-cell transcripts)), i.e. where spillover is most likely
W = 120
border = pd.read_parquet(xu.OUTS_DIR / "transcripts.parquet", columns=["feature_name", "x_location", "y_location", "qv"],
filters=[("feature_name", "in", ["EPCAM", "MALL", "CD3E", "TRAC"]), ("qv", ">=", 20)])
border["kind"] = np.where(border.feature_name.isin(["EPCAM", "MALL"]), "epi", "T")
border["tile"] = list(zip(border.x_location // W, border.y_location // W))
tile_counts = border.groupby(["tile", "kind"]).size().unstack(fill_value=0)
best = tile_counts.min(axis=1).idxmax()
X0, Y0 = best[0] * W, best[1] * W
print(f"Window: x {X0:.0f}-{X0 + W:.0f} um, y {Y0:.0f}-{Y0 + W:.0f} um "
f"({tile_counts.loc[[best], 'epi'].iloc[0]} epithelial and {tile_counts.loc[[best], 'T'].iloc[0]} T-cell transcripts)")
win = lambda xy: (xy[:, 0] > X0) & (xy[:, 0] < X0 + W) & (xy[:, 1] > Y0) & (xy[:, 1] < Y0 + W)
px = xu.PIXEL_SIZE_UM
store = tifffile.imread(xu.OUTS_DIR / "morphology_focus" / "morphology_focus_0000.ome.tif", aszarr=True, level=0)
img = zarr.open(store, mode="r") # reads only the tiles inside the window (JPEG 2000 via imagecodecs)
r0, c0, n = int(Y0 / px), int(X0 / px), int(W / px)
dapi = np.asarray(img[0, r0:r0 + n, c0:c0 + n] if img.ndim == 3 else img[r0:r0 + n, c0:c0 + n]).astype(float)
store.close()
lo, hi = np.percentile(dapi, [1, 99.7])
tx = pd.read_parquet(xu.OUTS_DIR / "transcripts.parquet", columns=["feature_name", "x_location", "y_location", "qv"],
filters=[("x_location", ">", X0), ("x_location", "<", X0 + W), ("y_location", ">", Y0), ("y_location", "<", Y0 + W)])
tx = tx[tx.qv >= 20]
SHOW = {"EPCAM": "#e8743b", "MALL": "#e8743b", "CD3E": "#3fa7f5", "TRAC": "#3fa7f5"}
b10 = pd.read_parquet(xu.OUTS_DIR / "cell_boundaries.parquet")
in_win_ids = b10.loc[win(b10[["vertex_x", "vertex_y"]].to_numpy()), "cell_id"].unique()
polys10 = [g[["vertex_x", "vertex_y"]].to_numpy() for _, g in b10[b10.cell_id.isin(in_win_ids)].groupby("cell_id", sort=False)]
gdf = pro_sd.shapes["cell_boundaries"]
cent = np.column_stack([gdf.geometry.centroid.x, gdf.geometry.centroid.y])
polysP = [np.asarray(part.exterior.coords) for geom in gdf.geometry[win(cent)] for part in getattr(geom, "geoms", [geom])]
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
for ax, polys, title in [(axes[0], polys10, "10x segmentation"), (axes[1], polysP, "Proseg re-segmentation")]:
ax.imshow(np.clip((dapi - lo) / (hi - lo), 0, 1), cmap="gray", extent=[X0, X0 + W, Y0 + W, Y0])
ax.add_collection(PolyCollection(polys, facecolors="none", edgecolors="#f5d547", linewidths=0.8))
for gene, col in SHOW.items():
t = tx[tx.feature_name == gene]
ax.scatter(t.x_location, t.y_location, s=5, c=col, linewidths=0)
ax.set_xlim(X0, X0 + W); ax.set_ylim(Y0 + W, Y0); ax.set_title(title, loc="left", fontsize=11)
ax.set_xticks([]); ax.set_yticks([])
axes[0].plot([X0 + 5, X0 + 25], [Y0 + W - 5] * 2, color="white", lw=3)
axes[0].text(X0 + 15, Y0 + W - 8, "20 um", color="white", ha="center", fontsize=9)
fig.legend(handles=[plt.Line2D([], [], marker="o", ls="", color="#e8743b", label="epithelial (EPCAM, MALL)"),
plt.Line2D([], [], marker="o", ls="", color="#3fa7f5", label="T cell (CD3E, TRAC)"),
plt.Line2D([], [], color="#f5d547", label="cell outline")], loc="lower center", ncol=3, frameon=False)
plt.tight_layout(rect=(0, 0.05, 1, 1))
plt.show()
print("Same tissue, same transcripts. 10x outlines are fixed shapes grown from the stain images; Proseg outlines\n"
"follow where each cell's transcripts are, so a T cell's outline stops before the neighbouring tumour cell's\n"
"EPCAM / MALL transcripts. Proseg works on a 1 um voxel grid, hence its stepped edges.")
Window: x 4200-4320 um, y 1920-2040 um (621 epithelial and 841 T-cell transcripts)
Same tissue, same transcripts. 10x outlines are fixed shapes grown from the stain images; Proseg outlines follow where each cell's transcripts are, so a T cell's outline stops before the neighbouring tumour cell's EPCAM / MALL transcripts. Proseg works on a 1 um voxel grid, hence its stepped edges.
Interpretation. The window (x 4,200-4,320 µm, y 1,920-2,040 µm) is a dense lymphocyte field (blue) above large tumour nuclei (orange). In the 10x panel, outlines are smooth shapes drawn around nuclei, and those at the boundary row include EPCAM / MALL transcripts. Proseg outlines follow the transcripts instead, so the boundary between the two populations runs along the edge of the epithelial transcript cloud. The stepped outlines also show the cost of this approach: cell shapes become less regular and harder to interpret morphologically.
Part D: Export the Proseg counts¶
Purpose: save the Proseg count matrix in the layout the core workflow expects, so later analyses can use it without reading the Proseg output again.
Approach: proseg.export_counts builds an AnnData indexed by the original 10x cell ids, with raw counts in .X and layers["counts"], 2D cell area, Proseg volume and total counts in .obs, centroids in µm in obsm["spatial"], and the method and version in uns["segmentation"].
Expected output: the number of cells and genes saved and the path of *_proseg_counts.h5ad.
Note: overwrites the exported file on each run. Proseg stores its full command line, including local file paths, in uns["proseg_run"]; only the version and method are copied, so the exported file contains no local paths.
# ---- Export: Proseg counts as an AnnData the core workflow can read ------------------------
# indexed by the 10x cell ids; Proseg's own uns['proseg_run'] stores its command line with local paths and is not copied
out = proseg.export_counts(pro, pro_area)
EXPORT = OUT_DIR / f"{xu.SAMPLE}_proseg_counts.h5ad"
out.write_h5ad(EXPORT, compression="gzip")
print(f"Saved {out.n_obs:,} cells x {out.n_vars} genes to {xu.display_path(EXPORT)}")
print("Next: a SEGMENTATION = 'proseg' switch in the core workflow reads this file instead of the 10x matrix\n"
"and re-runs Steps 4-10 with identical code, so the comparison isolates the effect of segmentation.")
Saved 155,468 cells x 377 genes to ~/data/xenium_lung/processed/resegmentation/Xenium_V1_humanLung_Cancer_FFPE_proseg_counts.h5ad Next: a SEGMENTATION = 'proseg' switch in the core workflow reads this file instead of the 10x matrix and re-runs Steps 4-10 with identical code, so the comparison isolates the effect of segmentation.
Interpretation. All 155,468 Proseg cells and 377 genes are exported. Part E reads this file, applies the core QC and clustering, and re-tests the core findings.
Part E: Does re-segmentation change the biology?¶
Parts A-D show that Proseg cells are cleaner. Part E asks the question that matters for the conclusions: which core findings hold on cleaner cells, which were segmentation artefacts, and what becomes newly visible?
| Part | Approach |
|---|---|
| E1 | Re-discover cell types on Proseg cells with the core recipe; then re-cluster T / NK cells and look for clean CD4 T and Treg groups |
| E2 | Same cells, same Step 6 labels, 10x vs Proseg counts: re-measure the Step 9 CD8 T-cell gradient and checkpoint expression |
| E3 | Re-test CD8 T exclusion and the TLS aggregates on Proseg cell positions |
| E4 | One table: holds / changes / removed / newly visible |
E1. Cell types on Proseg cells¶
Purpose: confirm that the Proseg cells support the same cell-type structure as the 10x cells, before using them to re-test specific findings.
Approach: apply the core QC rules (proseg.qc_proseg: at least 10 counts, at least 5 genes, area within 3 MAD on the log scale), then the Step 5 recipe (cluster.normalise_and_pca, 30 PCs, 15 neighbours, Leiden at resolution 0.5, seed 0). Compare clusters with the 10x-based leiden_0.5 by adjusted Rand index (ARI) on shared cells, and describe each cluster by its top Wilcoxon markers and its most frequent Step 6 label.
Expected output: the QC summary, the number of clusters, the ARI, and a table of clusters with top markers, main Step 6 label and its share. Inspect whether each Proseg cluster has a coherent marker set and a recognisable Step 6 counterpart.
Note: the control-fraction QC rule cannot be applied because Proseg outputs no control probes. UMAP is skipped (viewing only, ~3 min). Low label shares are expected for whole-lineage clusters, because Step 6 split those lineages into several subtypes.
# ---- E1. Re-discover cell types on the Proseg cells (same recipe as core Steps 4-5) ----------
from sklearn.metrics import adjusted_rand_score
warnings.filterwarnings("ignore", message="zero-centering") # scanpy notice when scaling sparse data
pc, (area_lo, area_hi) = proseg.qc_proseg(sc.read_h5ad(EXPORT)) # core QC rules (counts, genes, 3-MAD area)
print(f"QC (counts >= 10, genes >= 5, area {area_lo:.0f}-{area_hi:.0f} um^2): {pc.n_obs:,} Proseg cells kept "
f"(10x: 154,394). Proseg output has no control probes, so the control-fraction rule does not apply.")
cluster.normalise_and_pca(pc) # same recipe as core Step 5
cluster.cluster(pc, n_pcs=30, n_neighbors=15, resolutions=[0.5], seed=0, umap=False)
pc.obs["leiden"] = pc.obs["leiden_0.5"]
print(f"Leiden (resolution 0.5): {pc.obs['leiden'].nunique()} clusters (10x-based: 14). UMAP skipped: it is for viewing only.\n")
# Labels from the core workflow, carried across by cell id
step6 = labels.reindex(pc.obs_names)
tenx_clusters = xu.load_step("clustered", backed="r").obs["leiden_0.5"]
shared = pc.obs_names.intersection(tenx_clusters.index)
print(f"Cluster agreement with the 10x-based clustering (ARI, {len(shared):,} shared cells): "
f"{adjusted_rand_score(tenx_clusters[shared], pc.obs.loc[shared, 'leiden']):.2f}\n")
sc.tl.rank_genes_groups(pc, "leiden", method="wilcoxon", layer="lognorm", key_added="markers")
top = sc.get.rank_genes_groups_df(pc, group=None, key="markers").query("pvals_adj < 0.01 and logfoldchanges > 1")
top6 = top.groupby("group", observed=True).head(6).groupby("group", observed=True)["names"].apply(", ".join)
purity = (pd.crosstab(pc.obs["leiden"], step6["cell_type"], normalize="index") * 100)
summary = pd.DataFrame({
"cells": pc.obs["leiden"].value_counts().reindex(top6.index),
"top markers": top6,
"main Step 6 label": purity.idxmax(axis=1).reindex(top6.index),
"% of cluster": purity.max(axis=1).reindex(top6.index).round(0),
})
display(summary)
QC (counts >= 10, genes >= 5, area 12-399 um^2): 143,975 Proseg cells kept (10x: 154,394). Proseg output has no control probes, so the control-fraction rule does not apply. Leiden (resolution 0.5): 14 clusters (10x-based: 14). UMAP skipped: it is for viewing only. Cluster agreement with the 10x-based clustering (ARI, 141,396 shared cells): 0.52
| cells | top markers | main Step 6 label | % of cluster | |
|---|---|---|---|---|
| group | ||||
| 0 | 19208 | MALL, EPCAM, SFTA2, TCIM, PCSK2, CFB | Tumour epithelial (MALL/TCIM) | 72.0 |
| 1 | 24336 | CYP2B6, TCIM, TFPI, MALL, CAPN8, CFTR | Tumour epithelial (CYP2B6/CFTR) | 84.0 |
| 2 | 12181 | CYP4B1, CAPN8, MYC, MLPH, EHF, CYP2B6 | Tumour epithelial (MYC/CAPN8) | 79.0 |
| 3 | 1324 | FCGR3A, MCEMP1, AIF1, MRC1, CD68, C1orf162 | Myeloid / epithelial mixed | 63.0 |
| 4 | 12968 | MS4A6A, AIF1, MRC1, CD163, MS4A4A, C1orf162 | cDC2 | 17.0 |
| 5 | 22491 | KRT7, EPCAM, GPRC5A, MYC, STEAP4, MET | Tumour epithelial (MYC/CAPN8) | 59.0 |
| 6 | 16506 | LTBP2, MYLK, FBN1, PDGFRA, PDGFRB, COL5A2 | Fibroblast | 78.0 |
| 7 | 17937 | TRAC, CXCR4, IL7R, PTPRC, CD3E, CD2 | CD8 T (GZMK+) | 19.0 |
| 8 | 7808 | PECAM1, VWF, CAVIN1, CD93, CLEC14A, CD34 | Endothelial | 92.0 |
| 9 | 2009 | CXCR4, MS4A1, BANK1, CD79A, PTPRC, CYTIP | B cell | 83.0 |
| 10 | 2789 | TENT5C, MZB1, SLAMF7, PRDM1, FKBP11, CD79A | Plasma cell | 89.0 |
| 11 | 1390 | CPA3, MS4A2, KIT, IL1RL1, SLC18A2, GATA2 | Mast cell | 91.0 |
| 12 | 2309 | MYH11, ACTA2, MYLK, CNN1, CAVIN1, ADAMTS1 | Smooth muscle | 88.0 |
| 13 | 719 | C20orf85, TSPAN19, SNTN, C1orf194, C6orf118, C... | Ciliated epithelial | 96.0 |
Interpretation. 143,975 Proseg cells pass QC (10x: 154,394), and Leiden again finds 14 clusters. ARI with the 10x-based clustering is 0.52 on 141,396 shared cells, about the same as the 10x-based clustering's own agreement with the 10x vendor clusters (0.51), so the partition is similar but not identical.
- Most clusters map to one Step 6 type: ciliated (96%), endothelial (92%), mast (91%), plasma (89%), smooth muscle (88%), B cells (83%), fibroblasts (78%) and the three tumour programmes (72-84%).
- The T / NK cluster (7) and the CD163+ myeloid cluster (4) have low label shares (19% and 17%) because they contain several Step 6 subtypes; E1b subclusters the T / NK cells.
- Cluster 3 (FCGR3A, MCEMP1, AIF1, CD68) consists mainly of cells Step 6 called
Myeloid / epithelial mixed(63%), and its top markers are myeloid rather than epithelial. On Proseg counts, many of these cells read as myeloid, though C2 shows most still carry some epithelial transcripts.
E1b. The key test: clean CD4 T and Treg populations¶
Purpose: test whether re-segmentation resolves the population that the 10x segmentation could not: conventional CD4 T cells, plus Tregs as a separate group.
Approach: subcluster the T / NK-majority clusters with the Step 6c recipe (annotate.subcluster, resolution 0.6). For each subcluster, proseg.tnk_probe_rates gives the % of cells with transcripts of CD4, CD8A, FOXP3, CD3E / TRAC and five foreign lineages (epithelial, macrophage, endothelial, fibroblast, mregDC). proseg.clean_cd4_treg then applies fixed thresholds: CD4 at least 30%, T at least 60%, CD8A below 20%, every foreign lineage below 25%; FOXP3 below 25% for conventional CD4 T, at least 25% for Treg.
Expected output: a table of subclusters with probe rates, top markers and size, then the subclusters that pass each rule (or "none").
Note: the rule must include fibroblast and mregDC probes and the FOXP3 exclusion. An earlier version without them reported four false-positive "clean CD4 T" subclusters, so validation rules need the same scrutiny as the data.
# ---- The key test: do CD4 T cells form a clean group on Proseg cells? ------------------------
lineage_of_cluster = pd.crosstab(pc.obs["leiden"], step6["lineage"]).idxmax(axis=1)
tnk_clusters = lineage_of_cluster.index[lineage_of_cluster == "T / NK"]
tnk = annotate.subcluster(pc, "leiden", list(tnk_clusters), resolution=0.6) # same recipe as core Step 6c
PROBES, FOREIGN = proseg.TNK_PROBES, proseg.FOREIGN
rates = proseg.tnk_probe_rates(tnk) # % of cells per subcluster with each lineage's transcripts
sub_top = sc.get.rank_genes_groups_df(tnk, group=None, key="sub_markers").groupby("group", observed=True).head(6) \
.groupby("group", observed=True)["names"].apply(", ".join)
sub_table = rates.join(sub_top.rename("top markers")).join(tnk.obs["sub"].value_counts().rename("cells"))
display(sub_table)
# Conventional CD4 T: CD4-leaning, T-cell positive, not Treg, and no other lineage's markers at high rates
# CD4 >= 30%, T >= 60%, CD8A < 20%, every foreign lineage < 25%; Treg = FOXP3 >= 25%
clean_cd4, clean_treg = proseg.clean_cd4_treg(sub_table)
print(f"T/NK cells re-clustered: {tnk.n_obs:,} in {tnk.obs['sub'].nunique()} subclusters.")
print("Clean conventional CD4 T (CD4 >= 30%, T >= 60%, CD8A < 20%, FOXP3 < 25%, every foreign lineage < 25%): "
+ (", ".join(f"{i} ({int(r.cells):,} cells)" for i, r in clean_cd4.iterrows()) if len(clean_cd4) else "none"))
print("Clean Treg (same, FOXP3 >= 25%): "
+ (", ".join(f"{i} ({int(r.cells):,} cells)" for i, r in clean_treg.iterrows()) if len(clean_treg) else "none"))
| CD4 | CD8A | FOXP3 | T (CD3E/TRAC) | epithelial | macrophage | endothelial | fibroblast | mregDC | top markers | cells | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 21.0 | 34.0 | 4.0 | 91.0 | 0.0 | 4.0 | 10.0 | 6.0 | 3.0 | CCL5, GZMA, GZMK, CD3E, CXCR4, CD8A | 4065 |
| 1 | 30.0 | 13.0 | 10.0 | 85.0 | 0.0 | 6.0 | 7.0 | 11.0 | 26.0 | CCL19, CCR7, CXCL9, LAMP3, IL7R, BASP1 | 2216 |
| 10 | 28.0 | 15.0 | 13.0 | 83.0 | 7.0 | 5.0 | 14.0 | 10.0 | 5.0 | STEAP4, TCIM, GPRC5A, CAPN8, CD2, MEF2C | 146 |
| 2 | 34.0 | 19.0 | 14.0 | 89.0 | 1.0 | 11.0 | 19.0 | 56.0 | 2.0 | SFRP4, COL5A2, VCAN, FBN1, ACTA2, ASPN | 1408 |
| 3 | 31.0 | 16.0 | 6.0 | 91.0 | 1.0 | 11.0 | 37.0 | 21.0 | 8.0 | MMRN1, ECSCR, CD34, ADGRL4, TFPI, CLEC14A | 563 |
| 4 | 13.0 | 33.0 | 4.0 | 73.0 | 2.0 | 8.0 | 18.0 | 8.0 | 3.0 | KLRD1, GNLY, PRF1, KLRC1, CCL5, NKG7 | 554 |
| 5 | 34.0 | 17.0 | 7.0 | 84.0 | 0.0 | 10.0 | 19.0 | 35.0 | 3.0 | PTGDS, FBLN1, INMT, MYLK, C7, PMP22 | 3535 |
| 6 | 42.0 | 17.0 | 8.0 | 88.0 | 0.0 | 52.0 | 13.0 | 14.0 | 2.0 | AIF1, MRC1, MS4A4A, CD14, IGSF6, CD163 | 3203 |
| 7 | 40.0 | 7.0 | 53.0 | 95.0 | 0.0 | 7.0 | 8.0 | 10.0 | 3.0 | CTLA4, FOXP3, TRAC, CD28, IL2RA, PTPRC | 940 |
| 8 | 32.0 | 22.0 | 12.0 | 92.0 | 1.0 | 11.0 | 10.0 | 11.0 | 3.0 | PCNA, MKI67, TRAC, TOP2A, CDK1, CD3D | 747 |
| 9 | 40.0 | 22.0 | 8.0 | 84.0 | 0.0 | 12.0 | 15.0 | 16.0 | 4.0 | CPA3, MS4A2, IL1RL1, KIT, CD4, SLC18A2 | 560 |
T/NK cells re-clustered: 17,937 in 11 subclusters. Clean conventional CD4 T (CD4 >= 30%, T >= 60%, CD8A < 20%, FOXP3 < 25%, every foreign lineage < 25%): none Clean Treg (same, FOXP3 >= 25%): 7 (940 cells)
Interpretation.
- A clean Treg population emerges (subcluster 7, 940 cells: FOXP3 53%, top markers CTLA4, FOXP3, IL2RA; every foreign lineage at or below 10%).
- No clean conventional CD4 T population, even after re-segmentation. The CD4-leaning candidates still carry mregDC transcripts (subcluster 1: 26%, with CCL19, CCR7, LAMP3), fibroblast transcripts (subclusters 2 and 5: 56% and 35%), endothelial transcripts (subcluster 3: 37%) or macrophage transcripts (subcluster 6: 52%).
- Proseg removes tumour-to-immune spillover but not the overlaps in dense TLS and stroma, where immune, stromal and dendritic cells are tightly packed in a 5 µm 2D section. CD4 is also weakly detected on this panel (30-42% of cells in the CD4-leaning subclusters), which limits any CD4 T call.
E2. Same cells, same labels: per-cell findings with 10x vs Proseg counts¶
Purpose: decide whether the Step 9b CD8 T-cell gradient is biology or spillover. With 10x counts, CD8 T cells lost CCR7 towards the tumour, but epithelial transcripts rose in the same direction, so the two could not be separated.
Approach: keep the Step 6 labels and Step 9 compartments (TLS, stroma > 50 µm, border 15-50 µm, tumour contact ≤ 15 µm) fixed by cell id and change only the counts. proseg.cd8_by_segmentation reports % positive per compartment for CCR7, IL7R, GZMA, CTLA4 and HAVCR2, with EPCAM and MALL as the spillover control. A second table repeats the Step 9a checkpoint positivity (CD274, PDCD1, CTLA4, HAVCR2) within cell types, against the tumour noise floor.
Expected output: two tables with 10x and Proseg rows per gene. The deciding pattern: if EPCAM / MALL flatten on Proseg cells while CCR7 and GZMA keep their gradient, the gradient is not spillover.
Note: cells were selected and labelled on 10x counts, which biases the comparison slightly in favour of 10x; a gradient that survives is therefore a conservative result.
# ---- E2. Same cells, same labels, different segmentation: re-measure per-cell findings --------
# Labels and locations (TLS / stroma / border / tumour contact) come from the core workflow by cell id;
# only the expression changes (10x counts vs Proseg counts).
time_obs = xu.load_step("time", backed="r").obs[["cell_type", "cd8_compartment"]]
COMPARTMENTS = list(time_obs["cd8_compartment"].cat.categories)
GENES = ["CCR7", "IL7R", "GZMA", "CTLA4", "HAVCR2", "EPCAM", "MALL"]
e2, n_cd8 = proseg.cd8_by_segmentation(tenx, pro, time_obs, GENES)
display(e2.loc[GENES])
print(f"{n_cd8:,} CD8 T cells (Step 6 label), % positive by location. EPCAM / MALL are the spillover control:\n"
"if they flatten on Proseg cells while CCR7 / GZMA keep their gradient, the gradient is real biology.")
# Checkpoint expression (Step 9a): % positive within cell type, 10x vs Proseg
types = ["CXCL9+ macrophage", "mregDC (LAMP3+)", "CD163+ macrophage", "CD8 T (GZMK+)", "Treg"]
ann_types = labels["cell_type"]
tumour = ann_types.index[ann_types.isin(xu.TUMOUR_TYPES)].intersection(pro.obs_names)
chk = []
for gene in ["CD274", "PDCD1", "CTLA4", "HAVCR2"]:
for seg, adata in [("10x", tenx), ("Proseg", pro)]:
row = {"gene": gene, "segmentation": seg,
"Tumour (noise floor)": 100 * (sp.csr_matrix(adata[tumour, [gene]].X) > 0).mean()}
for t in types:
cells = ann_types.index[ann_types == t].intersection(pro.obs_names)
row[t] = 100 * (sp.csr_matrix(adata[cells, [gene]].X) > 0).mean()
chk.append(row)
display(pd.DataFrame(chk).set_index(["gene", "segmentation"]).round(1))
| TLS | Stroma (> 50 um) | Border (15-50 um) | Tumour contact (<= 15 um) | ||
|---|---|---|---|---|---|
| gene | segmentation | ||||
| CCR7 | 10x | 35.9 | 29.5 | 22.0 | 11.0 |
| Proseg | 24.8 | 19.6 | 17.0 | 9.4 | |
| IL7R | 10x | 73.6 | 69.4 | 67.5 | 52.9 |
| Proseg | 69.6 | 65.9 | 64.4 | 56.3 | |
| GZMA | 10x | 29.9 | 40.3 | 38.1 | 44.8 |
| Proseg | 29.8 | 40.1 | 37.7 | 48.4 | |
| CTLA4 | 10x | 36.7 | 28.3 | 25.0 | 19.5 |
| Proseg | 33.4 | 27.0 | 25.1 | 20.9 | |
| HAVCR2 | 10x | 3.8 | 5.2 | 5.0 | 7.6 |
| Proseg | 3.7 | 6.9 | 5.0 | 7.4 | |
| EPCAM | 10x | 1.4 | 1.7 | 4.2 | 8.1 |
| Proseg | 0.3 | 0.0 | 0.4 | 0.7 | |
| MALL | 10x | 3.2 | 3.0 | 7.1 | 13.9 |
| Proseg | 0.0 | 0.2 | 0.1 | 0.0 |
3,781 CD8 T cells (Step 6 label), % positive by location. EPCAM / MALL are the spillover control: if they flatten on Proseg cells while CCR7 / GZMA keep their gradient, the gradient is real biology.
| Tumour (noise floor) | CXCL9+ macrophage | mregDC (LAMP3+) | CD163+ macrophage | CD8 T (GZMK+) | Treg | ||
|---|---|---|---|---|---|---|---|
| gene | segmentation | ||||||
| CD274 | 10x | 0.9 | 12.0 | 13.4 | 1.3 | 2.0 | 1.6 |
| Proseg | 0.6 | 9.6 | 10.6 | 1.4 | 2.0 | 1.6 | |
| PDCD1 | 10x | 0.1 | 0.4 | 0.6 | 0.2 | 0.9 | 0.4 |
| Proseg | 0.1 | 0.4 | 0.4 | 0.1 | 1.0 | 0.6 | |
| CTLA4 | 10x | 2.7 | 5.6 | 18.8 | 5.7 | 27.3 | 69.5 |
| Proseg | 2.0 | 6.0 | 17.7 | 5.2 | 26.6 | 65.5 | |
| HAVCR2 | 10x | 2.8 | 32.3 | 8.0 | 24.9 | 5.2 | 8.2 |
| Proseg | 1.0 | 32.7 | 7.1 | 30.0 | 5.6 | 7.8 |
Interpretation.
- The spillover control flattens. At tumour contact, EPCAM falls from 8.1% to 0.7% and MALL from 13.9% to 0.0%; elsewhere both are near zero on Proseg cells.
- The CCR7 gradient persists (Proseg: 24.8% in TLS to 9.4% at contact). Absolute rates are lower than with 10x and the gradient is shallower (about 2.6x vs 3.3x), so part of the 10x CCR7 signal in TLS came from neighbouring cells.
- GZMA still rises at contact (29.8% to 48.4%), and CTLA4 still falls (33.4% to 20.9%), both close to the 10x values.
- The effector-like shift of CD8 T cells at the tumour border is therefore not a spillover artefact.
- Checkpoint findings are unchanged. PD-L1 (CD274) stays on CXCL9+ macrophages and mregDCs (9.6% and 10.6% on Proseg) and at the noise floor on tumour (0.6%); PD-1 stays at or near 1% in every group; CTLA4 stays highest in Tregs (65.5%).
E3. Spatial findings on Proseg cell positions¶
Purpose: check whether the core spatial findings (CD8 T exclusion, the myeloid / stromal tumour border, the TLS-like aggregates) depend on where 10x placed the cells.
Approach: take Proseg centroids for high-confidence Step 6 cells (proseg.positions_high_confidence). Recompute the median distance to the nearest tumour cell and the % within 15 µm for five cell types against the all-non-tumour baseline (proseg.exclusion_on_positions), and re-detect B-cell aggregates with the Step 7 DBSCAN settings (eps 25 µm, at least 20 B cells). Match each aggregate to the nearest 10x-based TLS centre.
Expected output: the baseline, a distance table per cell type, and a table matching each Proseg aggregate to a 10x TLS with its distance.
Note: Proseg starts from the 10x nuclei, so cell positions move only a few µm and labels are unchanged. This part tests robustness to small position shifts, not an independent re-discovery of the spatial structure.
# ---- E3. Spatial findings with Proseg cell positions ------------------------------------------
from scipy.spatial import cKDTree
xy, types_hc = proseg.positions_high_confidence(out, labels) # Proseg centroids of high-confidence cells
excl, base = proseg.exclusion_on_positions(xy, types_hc, xu.TUMOUR_TYPES,
["CD8 T (GZMK+)", "B cell", "Alveolar macrophage", "cDC2", "Fibroblast"])
print(f"Proseg positions: baseline (all non-tumour cells) {100 * base:.0f}% within 15 um of a tumour cell")
display(excl.round(1))
centres = proseg.aggregates_on_positions(xy, types_hc) # DBSCAN eps 25 um, min 20 B cells
py_tls = xu.load_step("spatial", backed="r")
py_obs = py_tls.obs[["tls_id"]].assign(x=py_tls.obsm["spatial"][:, 0], y=py_tls.obsm["spatial"][:, 1])
py_centres = py_obs[py_obs.tls_id != "none"].groupby("tls_id", observed=True)[["x", "y"]].mean()
dist_to_py = cKDTree(py_centres.to_numpy()).query(centres, k=1)
tls_match = pd.DataFrame({"Proseg aggregate": [f"A{i + 1}" for i in range(len(centres))],
"nearest 10x TLS": py_centres.index[dist_to_py[1]], "distance (um)": dist_to_py[0].round(0)})
print(f"\nB-cell aggregates on Proseg positions (DBSCAN eps 25 um, min 20): {len(centres)}; "
f"{(dist_to_py[0] <= 100).sum()} lie within 100 um of a 10x-based TLS")
display(tls_match)
Proseg positions: baseline (all non-tumour cells) 34% within 15 um of a tumour cell
| median distance (um) | % within 15 um | |
|---|---|---|
| cell type | ||
| CD8 T (GZMK+) | 37.1 | 12.2 |
| B cell | 45.0 | 7.4 |
| Alveolar macrophage | 12.2 | 68.0 |
| cDC2 | 14.1 | 54.1 |
| Fibroblast | 15.2 | 49.0 |
B-cell aggregates on Proseg positions (DBSCAN eps 25 um, min 20): 6; 6 lie within 100 um of a 10x-based TLS
| Proseg aggregate | nearest 10x TLS | distance (um) | |
|---|---|---|---|
| 0 | A1 | TLS1 | 3.0 |
| 1 | A2 | TLS2 | 3.0 |
| 2 | A3 | TLS3 | 9.0 |
| 3 | A4 | TLS4 | 7.0 |
| 4 | A5 | TLS5 | 5.0 |
| 5 | A6 | TLS6 | 7.0 |
Interpretation. CD8 T exclusion holds: 12.2% of CD8 T cells lie within 15 µm of tumour (median 37 µm) against a 34% baseline, as in the core workflow. Alveolar macrophages (68%), cDC2 (54%) and fibroblasts (49%) still dominate the tumour border, and B cells remain distant (7.4%). DBSCAN finds six B-cell aggregates, each within 3-9 µm of a 10x-based TLS centre.
E4. Which conclusions survive re-segmentation?¶
Purpose: collect the re-tests into one verdict per finding, so the effect of segmentation on each conclusion is explicit.
Approach: a table of 10x vs Proseg values for each finding, with a verdict from fixed rules: exclusion holds if the CD8 T contact rate stays below 0.8x the baseline, the TLS holds if at least five aggregates match a 10x TLS within 100 µm, and the CCR7 gradient holds if the TLS rate stays more than twice the contact rate.
Expected output: one row per finding with columns 10x segmentation, Proseg and verdict.
Note: the 10x column restates values from the core workflow as fixed text; the Proseg column is computed in this notebook.
# ---- E4. Which conclusions survive re-segmentation? -------------------------------------------
cd8_contact = excl.loc["CD8 T (GZMK+)", "% within 15 um"]
ccr7 = e2.loc[("CCR7", "Proseg")]; epcam10, epcamP = e2.loc[("EPCAM", "10x")], e2.loc[("EPCAM", "Proseg")]
verdict = pd.DataFrame([
["Tumour-to-immune transcript spillover", "T/NK cells with epithelial transcripts 18.8%",
f"{spill.iloc[0]['Proseg: % positive']:.1f}%", "reduced 10x+ (segmentation artefact)"],
["CD8 T cells excluded from tumour", "12-16% within 15 um vs 35% baseline",
f"{cd8_contact:.0f}% vs {100 * base:.0f}% baseline", "holds" if cd8_contact < 0.8 * 100 * base else "changes"],
["TLS-like B-cell aggregates", "6 aggregates", f"{len(centres)} aggregates, {(dist_to_py[0] <= 100).sum()} at 10x TLS positions",
"holds" if (dist_to_py[0] <= 100).sum() >= 5 else "changes"],
["CD8 T cells lose CCR7 toward tumour", "36% (TLS) -> 11% (contact)",
f"{ccr7['TLS']:.0f}% -> {ccr7['Tumour contact (<= 15 um)']:.0f}%",
"holds" if ccr7["TLS"] > 2 * ccr7["Tumour contact (<= 15 um)"] else "weakens"],
["Epithelial spillover in tumour-contact CD8 T cells", f"EPCAM {epcam10['Tumour contact (<= 15 um)']:.0f}% at contact",
f"EPCAM {epcamP['Tumour contact (<= 15 um)']:.0f}% at contact", "removed by re-segmentation"],
["Clean conventional CD4 T population", "none (always 'mixed')",
f"{len(clean_cd4)} clean subcluster(s)", "newly visible" if len(clean_cd4) else "still not resolved"],
["Clean Treg population", "Treg subcluster (Python) / merged with CD4 T (R)",
f"{len(clean_treg)} clean subcluster(s)", "resolved" if len(clean_treg) else "not resolved"],
], columns=["finding", "10x segmentation", "Proseg", "verdict"]).set_index("finding")
display(verdict)
| 10x segmentation | Proseg | verdict | |
|---|---|---|---|
| finding | |||
| Tumour-to-immune transcript spillover | T/NK cells with epithelial transcripts 18.8% | 1.5% | reduced 10x+ (segmentation artefact) |
| CD8 T cells excluded from tumour | 12-16% within 15 um vs 35% baseline | 12% vs 34% baseline | holds |
| TLS-like B-cell aggregates | 6 aggregates | 6 aggregates, 6 at 10x TLS positions | holds |
| CD8 T cells lose CCR7 toward tumour | 36% (TLS) -> 11% (contact) | 25% -> 9% | holds |
| Epithelial spillover in tumour-contact CD8 T cells | EPCAM 8% at contact | EPCAM 1% at contact | removed by re-segmentation |
| Clean conventional CD4 T population | none (always 'mixed') | 0 clean subcluster(s) | still not resolved |
| Clean Treg population | Treg subcluster (Python) / merged with CD4 T (R) | 1 clean subcluster(s) | resolved |
Interpretation. Every headline finding tested here holds on Proseg cells: CD8 T exclusion, the six TLS-like aggregates, the CCR7 gradient and the checkpoint pattern. What re-segmentation changes is noise. Tumour-to-immune spillover is largely removed, which turns the Step 9b gradient from "possibly an artefact" into a supported result and resolves a clean Treg population. Conventional CD4 T cells remain unresolved, a limitation of 2D segmentation in dense immune-stromal zones and of weak CD4 detection on this panel.
Limits: one section and one re-segmentation method; Proseg initialises from the 10x nuclei, so it corrects transcript assignment but not missed or merged nuclei; labels and compartments come from the 10x-based core workflow.