03. Reference mapping: validate the manual cell types against a lung cancer atlas

Author: MD Shakhaowat Hossain, PhD Candidate, Tulane University School of Medicine · Project website · Code

Why: the cell types in the core workflow (Step 6) were assigned by hand from cluster marker genes. That is standard practice, but it rests on one analyst's judgement. Here every cell is also labelled by a classifier trained on an independent, published reference, and the two labelings are compared cell type by cell type. Agreement turns "labelled by hand" into "supported by a reference"; disagreement shows where the manual labels, the panel or the segmentation are weak.

Reference: the core Lung Cancer Atlas (LuCA; Salcher et al., Cancer Cell 2022), an expert-annotated single-cell RNA-seq atlas of non-small cell lung cancer and normal lung (892,296 cells, 298 donors), downloaded from CELLxGENE (CC BY 4.0).

Method: the CellTypist recipe (Domínguez Conde et al., Science 2022): logistic regression on log-normalised, scaled expression, trained only on the genes shared by the atlas and the Xenium panel.

Part	Content
A	Prepare the reference: download, restrict to panel genes, check sequencing depth
B	Train the classifier and measure its ceiling on held-out atlas donors
C	Annotate both segmentations (10x and Proseg, from notebook 02)
D	Agreement with the manual labels: lineage, cell type, segmentation, causes of disagreement
E	Export and conclusions

Setup: imports, paths and output folders

Purpose: load the shared ist_analysis package and define where the reference and the results are stored, so paths match the core notebook and the Snakemake pipeline.

Approach: import proseg, reference and utils (as xu) from ist_analysis (repository root; pip install -e .), check the kernel, and create the reference folder (xu.DATA_ROOT / "reference") and the output folder (reference_mapping).

Expected output: both folders, printed as ~/... paths.

Note: this notebook reads the Step 6 object from the core notebook and the Proseg export from notebook 02; run both first.

from ist_analysis import proseg, reference
from ist_analysis import utils as xu   # paths and constants shared with the core workflow and the pipeline

xu.check_kernel()

import pickle
import time
import warnings

import anndata as ad
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


warnings.filterwarnings("ignore", category=FutureWarning)
REF_DIR = xu.DATA_ROOT / "reference"
REF_DIR.mkdir(exist_ok=True)
OUT_DIR = xu.output_dir("reference_mapping")
print(f"Reference files: {xu.display_path(REF_DIR)}\nOutputs: {xu.display_path(OUT_DIR)}")

Reference files: ~/data/xenium_lung/reference
Outputs: ~/data/xenium_lung/processed/reference_mapping

Part A: Prepare the reference

A1. Download the atlas and restrict it to the panel genes

Purpose: obtain the reference and reduce it to the genes the Xenium panel measures, since a classifier for Xenium cells can only use genes both data sets share.

Approach: download the LuCA core atlas .h5ad from CELLxGENE if requested, then stream it in blocks of 50,000 cells (reference.build_panel_reference), keep raw counts (raw.X) for the panel genes and cache the result. Genes are matched by Ensembl ID, not symbol, so renamed genes are not lost.

Expected output: atlas size (cells, donors, datasets), the number of panel genes found, and the list of absent genes. Check that the absent genes are not markers needed for lung cell types.

Note: the atlas is a 12.9 GB download; DOWNLOAD_ATLAS = False (default) stops with an error if neither the atlas nor the cache exists. Building the panel-gene cache (~80 MB) takes ~17 min and runs only if the cache is missing. anndata's backed mode was far slower than block streaming for this row selection.

import requests
from tqdm.auto import tqdm

ATLAS_URL = "https://datasets.cellxgene.cziscience.com/46e0287b-9a33-4e83-99f3-8c044131bfdc.h5ad"   # LuCA core atlas
ATLAS = REF_DIR / "luca_core_atlas.h5ad"
PANEL_REF = REF_DIR / "luca_core_panel_genes.h5ad"
DOWNLOAD_ATLAS = False   # True = download (12.9 GB) when the file is missing

if not ATLAS.exists() and not PANEL_REF.exists():
    if not DOWNLOAD_ATLAS:
        raise FileNotFoundError(f"{xu.display_path(ATLAS)} not found: set DOWNLOAD_ATLAS = True to download it (12.9 GB).")
    with requests.get(ATLAS_URL, stream=True, timeout=60) as r, open(ATLAS, "wb") as fh:
        r.raise_for_status()
        with tqdm(total=int(r.headers["Content-Length"]), unit="B", unit_scale=True, desc="LuCA core atlas") as bar:
            for chunk in r.iter_content(chunk_size=8 << 20):
                bar.update(fh.write(chunk))

ann = xu.load_step("annotated")   # Step 6: manual labels + raw counts for all QC-passing cells

if not PANEL_REF.exists():
    t0 = time.time()
    # stream the atlas in 50,000-cell blocks, keep the panel genes (matched by Ensembl id), cache the result
    reference.build_panel_reference(ATLAS, ann.var["gene_ids"], PANEL_REF)
    print(f"Built the panel-gene cache in {(time.time() - t0) / 60:.0f} min")

panel_ref = ad.read_h5ad(PANEL_REF)
missing = ann.var_names[~ann.var_names.isin(panel_ref.var_names)]
print(f"Atlas: {panel_ref.n_obs:,} cells, {panel_ref.obs['donor_id'].nunique()} donors, {panel_ref.obs['dataset'].nunique()} datasets")
print(f"Panel genes found in the atlas: {panel_ref.n_vars} of {ann.n_vars}")
print(f"Absent ({len(missing)}, non-lung markers of the multi-tissue panel): {', '.join(missing)}")

Atlas: 892,296 cells, 298 donors, 21 datasets
Panel genes found in the atlas: 354 of 377
Absent (23, non-lung markers of the multi-tissue panel): ADIPOQ, AMY2A, APOA5, ASCL3, CCL27, CFHR1, CLCA1, FCN2, GCG, GLYATL1, GYPA, INS, KRT20, MLANA, MTRNR2L11, NAT8, OPRPN, SCGN, SLC22A8, SLC26A3, SST, TMEM174, UMOD

Interpretation. The atlas holds 892,296 cells from 298 donors in 21 datasets. 354 of the 377 panel genes are present. The 23 absent genes are markers of non-lung tissues on this multi-tissue panel (for example INS, GCG, UMOD, MLANA, GYPA), so no lung cell type loses a defining marker.

A2. Training cells and sequencing depth

Purpose: choose atlas cells comparable to the Xenium section and check that both data sets have similar depth on the shared genes. A classifier trained on much deeper cells would expect markers that Xenium rarely detects, and the reference would then need down-sampling first.

Approach: keep atlas cells from lung tissue (no lymph node or metastasis sites), from UMI-based platforms (Smart-seq2 counts reads, not molecules) and with a label (cell_type_major, excluding other), via reference.eligible. Compare the distribution of counts per cell on the 354 shared genes between eligible atlas cells and Xenium cells (10x segmentation).

Expected output: the number of eligible training cells and cell types, and a percentile table of counts per cell. Compare the medians and the spread.

Note: depth on the shared genes is the relevant comparison, not total depth: scRNA-seq spreads its counts over ~18,000 genes, while Xenium spends all of its counts on the panel.

o = panel_ref.obs
eligible = reference.eligible(o)   # lung tissue, UMI platforms, labelled cells
print(f"Eligible training cells: {eligible.sum():,} of {len(o):,} "
      f"(lung tissue, UMI platforms, labelled), {o.loc[eligible, 'cell_type_major'].nunique()} cell types\n")

shared_genes = panel_ref.var_names
ref_depth = np.asarray(panel_ref[eligible.to_numpy()].X.sum(axis=1)).ravel()
xen_depth = np.asarray(ann[:, shared_genes].layers["counts"].sum(axis=1)).ravel()
depth = pd.DataFrame({f"{name}": np.percentile(v, [10, 25, 50, 75, 90]) for name, v in
                      [("LuCA atlas (scRNA-seq)", ref_depth), ("Xenium (10x segmentation)", xen_depth)]},
                     index=["10th", "25th", "median", "75th", "90th"]).T
display(depth.style.format("{:.0f}").set_caption("Counts per cell on the shared panel genes (percentiles)"))
print("The medians match: the atlas has about as many counts per cell on these genes as Xenium does,\n"
      "so no down-sampling is needed. (scRNA-seq spreads its counts over ~18,000 genes; Xenium spends\n"
      "all of its counts on the panel, so the two end up at similar depth on the panel genes.)")

Eligible training cells: 756,622 of 892,296 (lung tissue, UMI platforms, labelled), 23 cell types

Counts per cell on the shared panel genes (percentiles)

	10th	25th	median	75th	90th
LuCA atlas (scRNA-seq)	14	25	47	89	153
Xenium (10x segmentation)	19	29	47	76	113

The medians match: the atlas has about as many counts per cell on these genes as Xenium does,
so no down-sampling is needed. (scRNA-seq spreads its counts over ~18,000 genes; Xenium spends
all of its counts on the panel, so the two end up at similar depth on the panel genes.)

Interpretation. 756,622 atlas cells in 23 types are eligible. Median depth on the shared genes is identical (47 counts per cell in both), so no down-sampling is needed. The atlas has the wider spread (10th-90th percentile 14-153 vs 19-113), so the training cells cover the full Xenium depth range.

Part B: Train the classifier

B1. Training, with held-out donors as a ceiling

Purpose: train the reference classifier and measure how well 354 panel genes can separate the 23 atlas types at all, before any Xenium-specific problem (segmentation, platform) enters. This ceiling sets the best agreement a manual type can reach in Part D.

Approach:

• Donor split: 20% of donors (not cells) are held out (reference.donor_split), because cells from one donor resemble each other and a random-cell split would overstate accuracy.
• Balanced training set: at most 5,000 cells per type, so abundant types (alveolar macrophages: 114,163 eligible cells) do not drown rare ones (mature DCs: 1,928).
• Normalisation: each cell scaled to 10,000 counts over the shared genes only, then log1p, identically for atlas and Xenium cells (reference.lognorm).
• Model: the CellTypist recipe (scale each gene, L2-regularised multinomial logistic regression), fit directly with scikit-learn (reference.train_model). The model is saved as luca_panel_genes_model.pkl.

Expected output: training size and time, overall accuracy on unseen donors, and a bar chart of per-type accuracy (dashed line = overall). Inspect which types fall well below the line.

Note: the celltypist package itself is not called. Its default trainer (v1.7.1) passes an option (multi_class) that scikit-learn 1.9 removed, and its SGD fallback fits one-vs-rest models whose probabilities were not calibrated on Xenium cells. A single multinomial model gives probabilities that sum to 1 across the 23 types. Overwrites the saved model.

CAP = 5000   # training cells per cell type, at most
train_idx, test_idx, test_donors, donors = reference.donor_split(o, eligible, test_fraction=0.2, max_per_type=CAP)
lognorm = reference.lognorm   # raw counts on the shared genes -> 10,000 per cell -> log1p (CellTypist's input)

t0 = time.time()
model = reference.train_model(lognorm(panel_ref[train_idx], shared_genes), o.loc[train_idx, "cell_type_major"].astype(str).to_numpy())
print(f"Trained on {len(train_idx):,} cells ({len(donors) - len(test_donors)} donors) in {time.time() - t0:.0f} s")
with open(OUT_DIR / "luca_panel_genes_model.pkl", "wb") as fh:
    pickle.dump({"model": model, "genes": list(shared_genes)}, fh)

held = reference.predict(model, lognorm(panel_ref[test_idx], shared_genes), test_idx)["atlas_label"].to_numpy()
truth = o.loc[test_idx, "cell_type_major"].astype(str).to_numpy()
ceiling = pd.Series(held == truth).groupby(truth).mean().sort_values()
print(f"Accuracy on {len(test_idx):,} cells from {len(test_donors)} unseen donors: {(held == truth).mean():.0%}\n")

fig, ax = plt.subplots(figsize=(6, 5.5))
ax.barh(ceiling.index, 100 * ceiling.to_numpy(), color="#7f8c8d")
ax.axvline(100 * (held == truth).mean(), color="black", lw=1, ls="--")
ax.set_xlabel("% correct on unseen donors")
ax.set_title("Ceiling: how well 354 panel genes separate atlas cell types")
ax.set_xlim(0, 100)
plt.tight_layout()
plt.show()

Trained on 101,745 cells (180 donors) in 7 s
Accuracy on 19,110 cells from 45 unseen donors: 87%

Interpretation. Trained on 101,745 cells from 180 donors in 7 s, the classifier labels 87% of 19,110 cells from 45 unseen donors correctly. Types with clean markers on the panel are recognised at 91-100% (stromal, endothelial, ciliated, B, mast, NK, plasma, neutrophils, AT1 and AT2). The hard ones are types that differ in degree rather than kind: transitional club/AT2 (64%), CD4 T (68%), Treg (75%), tumour cells (78%) and macrophages (81%). Disagreement between those types in Part D is expected from the panel alone and is not, by itself, evidence that the manual labels are wrong.

Part C: Annotate the Xenium cells

Purpose: label every Xenium cell with the reference classifier, in both segmentations, so the manual labels and the effect of segmentation can be tested against the same model.

Approach: apply the model (reference.predict) after the same normalisation as in training to:

• 10x segmentation: the 154,394 QC-passing cells of the core workflow (Step 6);
• Proseg: the re-segmented cells from notebook 02, after the same QC as notebook 02 Part E1 (proseg.qc_proseg).

Each cell receives the most probable atlas type and its probability.

Expected output: a summary (cells annotated, median probability, % of cells with probability ≥ 0.5) and the number of cells per atlas type in each segmentation. Compare confidence between segmentations.

Note: probabilities express confidence among the 23 atlas types only; a cell type absent from the atlas (for example lymphatic endothelium) is forced onto the nearest available type.

# ---- 10x segmentation: Step 6 cells
atlas10 = reference.predict(model, lognorm(ann, shared_genes), ann.obs_names)

# ---- Proseg segmentation (notebook 02 export), same QC as notebook 02 Part E1
pro, _ = proseg.qc_proseg(sc.read_h5ad(xu.output_dir("resegmentation") / f"{xu.SAMPLE}_proseg_counts.h5ad"))
atlasP = reference.predict(model, lognorm(pro, shared_genes), pro.obs_names)

summary = pd.DataFrame({
    "cells annotated": [len(atlas10), len(atlasP)],
    "median probability": [atlas10["atlas_prob"].median(), atlasP["atlas_prob"].median()],
    "cells with probability >= 0.5": [(atlas10["atlas_prob"] >= 0.5).mean(), (atlasP["atlas_prob"] >= 0.5).mean()],
}, index=["10x segmentation", "Proseg"])
display(summary.style.format({"cells annotated": "{:,}", "median probability": "{:.2f}",
                              "cells with probability >= 0.5": "{:.0%}"}))
counts = pd.DataFrame({"10x": atlas10["atlas_label"].value_counts(), "Proseg": atlasP["atlas_label"].value_counts()}).fillna(0).astype(int)
display(counts.sort_values("10x", ascending=False).T)

	cells annotated	median probability	cells with probability >= 0.5
10x segmentation	154,394	0.85	91%
Proseg	143,975	0.86	92%

atlas_label	transitional club/AT2	Tumor cells	Stromal	Endothelial cell	Plasma cell	Macrophage alveolar	Alveolar cell type 2	Macrophage	T cell CD4	T cell regulatory	...	B cell	cDC2	T cell CD8	Mast cell	DC mature	cDC1	Alveolar cell type 1	Neutrophils	pDC	NK cell
10x	47593	31572	22128	6796	6049	6029	5701	5650	5072	2616	...	2018	1629	1094	967	950	455	395	195	160	118
Proseg	44134	27402	20429	6948	5152	3933	4859	6779	6346	3680	...	1673	2162	1596	1010	585	374	328	199	107	148

2 rows × 23 columns

Interpretation. Both segmentations are labelled with similar confidence: median probability 0.85 (10x) and 0.86 (Proseg), with 91% and 92% of cells at ≥ 0.5. Confidence alone therefore does not distinguish the segmentations; Part D tests whether the labels agree with the manual annotation.

Part D: Agreement with the manual labels

D1. Lineage level

Purpose: test the broadest claim first: do the manual lineages match an independent reference?

Approach: map atlas types to the 7 Step 6 lineages (reference.ATLAS_LINEAGE; for example Tumor cells, Alveolar cell type 2 and Ciliated are Epithelial) and cross-tabulate against the manual lineage for high-confidence cells only, since mixed cells are by definition not one cell type.

Expected output: a heatmap of % of each manual lineage (rows) by atlas lineage (columns), and the overall agreement. Read across each row to see where the atlas places those cells.

Note: the overall figure is weighted by abundance (epithelium is about 60% of cells), so per-row agreement matters more than the total.

ATLAS_LINEAGE = reference.ATLAS_LINEAGE   # atlas type -> the 7 Step 6 lineages
assert set(ATLAS_LINEAGE) == set(model.classes_), "every atlas type needs a lineage"
LINEAGES = list(xu.LINEAGE_COLORS)

manual = ann.obs[["cell_type", "lineage", "annotation_confidence"]].astype(str)
high = manual.index[manual["annotation_confidence"] == "high"]
both = pd.DataFrame({"manual": manual.loc[high, "lineage"],
                     "atlas": atlas10.loc[high, "atlas_label"].map(ATLAS_LINEAGE)})
lin_tab = pd.crosstab(both["manual"], both["atlas"], normalize="index").reindex(index=LINEAGES, columns=LINEAGES).fillna(0) * 100

fig, ax = plt.subplots(figsize=(6.5, 5))
im = ax.imshow(lin_tab.to_numpy(), cmap="Blues", vmin=0, vmax=100)
ax.set_xticks(range(len(LINEAGES)), LINEAGES, rotation=45, ha="right")
ax.set_yticks(range(len(LINEAGES)), LINEAGES)
for i in range(len(LINEAGES)):
    for j in range(len(LINEAGES)):
        v = lin_tab.iat[i, j]
        if v >= 1:
            ax.text(j, i, f"{v:.0f}", ha="center", va="center", fontsize=8, color="white" if v > 60 else "black")
ax.set_xlabel("Atlas (LuCA) lineage")
ax.set_ylabel("Manual lineage (Step 6)")
ax.set_title("% of each manual lineage, by atlas lineage (10x cells)")
fig.colorbar(im, ax=ax, shrink=0.8, label="% of row")
plt.tight_layout()
plt.show()
print(f"Lineage agreement on {len(high):,} high-confidence cells: {(both['manual'] == both['atlas']).mean():.1%}")

Lineage agreement on 139,129 high-confidence cells: 92.0%

Interpretation. 92.0% of 139,129 high-confidence cells receive the same lineage from the atlas.

• Strong agreement: epithelial (99%), myeloid (91%), stromal (90%) and B / plasma (83%).
• Weaker: endothelial (67%, with 18% called stromal), T / NK (63%, with 15% called B / plasma and 15% myeloid) and mast (36%, with 24% called B / plasma and 15% T / NK).
• The weaker rows involve cells that sit among other immune or stromal cells; D4 tests whether spillover or the platform explains them.

D2. Cell type by cell type

Purpose: test each manual cell type individually, because lineage agreement can hide failures at the cell-type level.

Approach: for every manual type, the atlas types that would confirm it were written down before comparing (reference.CONFIRMED_BY; for example CD163+ macrophage is confirmed by Macrophage, Alveolar macrophage by Macrophage alveolar, the S100A12+ monocyte group by Monocyte or Neutrophils). reference.agreement reports, per type, the % of cells the atlas confirms, the most common atlas call, the median probability, and the atlas ceiling from Part B for the matching type.

Expected output: a table with one row per manual type, coloured by % confirmed. Compare each type's agreement with its ceiling: a manual type cannot be confirmed more often than the atlas recognises its own version of that type.

Note: the assertion stops the cell if any manual type lacks a written expectation, so no type is judged by a rule chosen after seeing the result.

CONFIRMED_BY = reference.CONFIRMED_BY     # manual type (Step 6) -> atlas types that agree with it
ann_types = set(manual.loc[high, "cell_type"])
assert ann_types <= set(CONFIRMED_BY), f"no expectation written for: {ann_types - set(CONFIRMED_BY)}"


def agreement(atlas, cells):
    """Per manual type: % confirmed by the atlas, most common call, median probability, atlas ceiling."""
    return reference.agreement(atlas, manual, cells, ceiling)


agree10 = agreement(atlas10, high)
display(agree10.style.format({"cells": "{:,}", "% confirmed by atlas": "{:.0f}", "median probability": "{:.2f}",
                              "atlas ceiling (%)": "{:.0f}"})
        .background_gradient(subset=["% confirmed by atlas"], cmap="RdYlGn", vmin=0, vmax=100))

	lineage	cells	% confirmed by atlas	most common atlas call	median probability	atlas ceiling (%)
manual type
Tumour epithelial (MALL/TCIM)	Epithelial	25,301	30	transitional club/AT2 (61%)	0.81	78
Tumour epithelial (CYP2B6/CFTR)	Epithelial	27,886	20	transitional club/AT2 (65%)	0.80	78
Tumour epithelial (MYC/CAPN8)	Epithelial	27,699	49	Tumor cells (49%)	0.81	78
Proliferating tumour	Epithelial	3,826	87	Tumor cells (87%)	0.98	78
Ciliated epithelial	Epithelial	1,028	87	Ciliated (87%)	1.00	98
CD8 T (GZMK+)	T / NK	3,811	18	T cell CD4 (29%)	0.67	87
Treg	T / NK	681	52	T cell regulatory (52%)	0.87	75
NK	T / NK	633	16	T cell CD8 (18%)	0.73	95
Proliferating T	T / NK	872	56	T cell regulatory (28%)	0.72	77
B cell	B / plasma	2,313	80	B cell (80%)	0.98	97
Plasma cell	B / plasma	3,607	81	Plasma cell (81%)	1.00	94
CD163+ macrophage	Myeloid	2,467	61	Macrophage (61%)	0.70	81
FCGR1A+ macrophage	Myeloid	1,918	66	Macrophage (66%)	0.79	81
LYVE1+ macrophage	Myeloid	368	54	Macrophage (54%)	0.79	81
CXCL9+ macrophage	Myeloid	261	49	Macrophage (49%)	0.75	81
Alveolar macrophage	Myeloid	1,925	84	Macrophage alveolar (84%)	0.92	86
Proliferating myeloid	Myeloid	414	57	Macrophage alveolar (31%)	0.76	84
Monocyte	Myeloid	202	72	Neutrophils (46%)	0.81	88
cDC1 (IRF8+)	Myeloid	618	39	cDC1 (39%)	0.72	86
cDC2	Myeloid	2,653	40	cDC2 (40%)	0.74	83
mregDC (LAMP3+)	Myeloid	470	77	DC mature (77%)	0.89	85
Mast cell	Mast	2,443	36	Mast cell (36%)	0.86	97
Fibroblast	Stromal	15,944	90	Stromal (90%)	1.00	100
Smooth muscle	Stromal	2,379	92	Stromal (92%)	1.00	100
Endothelial	Endothelial	9,094	69	Endothelial cell (69%)	1.00	98
Lymphatic endothelial	Endothelial	316	14	T cell CD4 (31%)	0.73	98

Interpretation.

• Confirmed (≥ 75%, close to the ceiling): smooth muscle (92%), fibroblast (90%), ciliated (87%), proliferating tumour (87%), alveolar macrophage (84%), plasma (81%), B (80%) and mregDC (77%). The S100A12+ monocyte group is confirmed at 72%, most often as Neutrophils (46%), consistent with the R workflow's SingleR result that this group is neutrophil-like.
• Right lineage, finer call uncertain (39-69%): the macrophage subtypes (49-66%), Treg (52%), proliferating T and myeloid cells (56-57%), cDC1 and cDC2 (39-40%). Their ceilings (75-86%) show these types are hard to separate on 354 genes. Endothelial cells (69%) are the exception: their ceiling is 98%, and most of the shortfall is cells called Stromal (18%), so the endothelial-stromal boundary rather than the panel limits agreement there.
• Tumour subtypes: the atlas calls most MALL/TCIM and CYP2B6/CFTR cells transitional club/AT2 (61-65%), and confirms only the MYC/CAPN8 programme as Tumor cells (49%). This does not contradict the manual call. LuCA identified malignant cells by copy-number inference on whole transcriptomes, which a 354-gene panel cannot reproduce, and 77% of its transitional club/AT2 cells come from primary tumour samples. Expression confirms these cells are lung epithelium with an alveolar / club-like programme; whether they are malignant needs copy-number evidence.
• Poorly confirmed (< 40%): CD8 T (18%), NK (16%), mast (36%) and lymphatic endothelial (14%). D4 examines why.

D3. Does cleaner segmentation raise agreement?

Purpose: test whether spillover was lowering agreement, by holding everything fixed except the counts.

Approach: pair cells by their 10x id and recompute D2's % confirmed with 10x counts and with Proseg counts, using the same manual labels and the same classifier. A second table follows the Step 6 mixed cells and reports the most common atlas lineage under each segmentation.

Expected output: a per-type table sorted by change in points (negative = lower agreement with Proseg counts), and a table of the mixed groups with their top atlas lineage before and after.

Note: this comparison is biased in favour of 10x. The manual labels were made from the 10x counts, so a cell was called mregDC partly because its 10x counts contained LAMP3; measured any other way, the same cells regress toward the average. The mixed cells, which were not selected for a clean marker, are the fairer test.

paired = manual.index.intersection(atlasP.index)
print(f"{len(paired):,} cells present in both segmentations\n")
seg = pd.DataFrame({
    "% confirmed, 10x counts": agreement(atlas10, paired.intersection(high))["% confirmed by atlas"],
    "% confirmed, Proseg counts": agreement(atlasP, paired.intersection(high))["% confirmed by atlas"],
})
seg["change (points)"] = seg.iloc[:, 1] - seg.iloc[:, 0]
seg.insert(0, "cells", agreement(atlas10, paired.intersection(high))["cells"])
display(seg.sort_values("change (points)").style.format({"cells": "{:,}", "% confirmed, 10x counts": "{:.0f}",
                                                         "% confirmed, Proseg counts": "{:.0f}", "change (points)": "{:+.0f}"})
        .background_gradient(subset=["change (points)"], cmap="RdBu", vmin=-30, vmax=30))

# Step 6 'mixed' cells: atlas lineage with 10x vs Proseg counts
mixed = paired[manual.loc[paired, "annotation_confidence"].to_numpy() == "mixed"]
mix = pd.DataFrame({"manual": manual.loc[mixed, "cell_type"],
                    "10x": atlas10.loc[mixed, "atlas_label"].map(ATLAS_LINEAGE),
                    "Proseg": atlasP.loc[mixed, "atlas_label"].map(ATLAS_LINEAGE)})
mix_tab = (mix.groupby("manual")[["10x", "Proseg"]]
           .agg(lambda s: f"{s.value_counts(normalize=True).index[0]} ({100 * s.value_counts(normalize=True).iloc[0]:.0f}%)"))
mix_tab.insert(0, "cells", mix["manual"].value_counts())
display(mix_tab.sort_values("cells", ascending=False))

141,396 cells present in both segmentations

	cells	% confirmed, 10x counts	% confirmed, Proseg counts	change (points)
manual type
mregDC (LAMP3+)	431	78	51	-27
Ciliated epithelial	930	87	73	-14
B cell	2,110	79	66	-14
Alveolar macrophage	1,807	84	74	-10
Monocyte	183	72	62	-10
Plasma cell	3,398	81	72	-9
Lymphatic endothelial	290	15	8	-7
Proliferating tumour	3,347	88	82	-6
Fibroblast	14,943	91	86	-6
Tumour epithelial (CYP2B6/CFTR)	25,482	21	17	-4
cDC1 (IRF8+)	590	39	37	-3
Proliferating myeloid	389	57	55	-2
Tumour epithelial (MALL/TCIM)	21,940	30	29	-1
Smooth muscle	2,128	93	92	-1
NK	609	16	16	+0
Tumour epithelial (MYC/CAPN8)	25,373	49	50	+0
CXCL9+ macrophage	236	50	50	+0
CD8 T (GZMK+)	3,565	18	19	+1
Mast cell	2,320	36	38	+2
Endothelial	8,655	70	72	+2
Treg	648	53	55	+2
CD163+ macrophage	2,332	62	71	+8
LYVE1+ macrophage	348	55	64	+8
FCGR1A+ macrophage	1,735	67	76	+10
cDC2	2,514	41	52	+11
Proliferating T	830	57	71	+14

	cells	10x	Proseg
manual
Myeloid / epithelial mixed	2905	Epithelial (74%)	Epithelial (79%)
T / myeloid mixed	2795	Myeloid (43%)	T / NK (52%)
CD4 T	2409	T / NK (72%)	T / NK (75%)
Myeloid / fibroblast mixed	1974	Myeloid (52%)	Myeloid (64%)
T / fibroblast mixed	1739	Stromal (42%)	T / NK (52%)
Myeloid / T mixed	1561	Myeloid (81%)	Myeloid (68%)
T / epithelial mixed	880	Epithelial (41%)	T / NK (65%)

Interpretation. Two opposite effects appear in the first table.

• Selection bias against Proseg. Types defined by one strong marker drop most: mregDC (−27 points), ciliated and B cells (−14). The defining markers themselves are detected less often in the same cells with Proseg counts (LAMP3 in mregDCs 92% → 63%, MS4A1 in B cells 93% → 64%), which is the regression expected from selecting cells on their 10x counts. These drops are therefore not evidence that Proseg is worse.
• Spillover removal in favour of Proseg. Types embedded among other cells gain despite that bias: macrophage subtypes (+8 to +10), cDC2 (+11) and proliferating T cells (+14).
• The mixed cells support Proseg. With Proseg counts, the T / epithelial, T / fibroblast and T / myeloid mixed groups are called T / NK most often (52-65%), where 10x counts called them epithelial, stromal or myeloid. Removing spillover lets the reference recognise the T cell inside the mixed outline. Myeloid / epithelial mixed cells remain mostly epithelial (74% → 79%).

D4. Why do some types disagree? Spillover vs platform

Purpose: separate two causes of poor agreement that need different responses:

• Spillover (segmentation): Xenium cells contain transcripts of a neighbouring lineage that the atlas version of the type does not express;
• Platform difference (chemistry): a defining marker is detected far less often by Xenium probes than by scRNA-seq, so the cell looks like a different type to a model trained on scRNA-seq.

Approach: for four poorly confirmed types, compare the % of cells with at least one count of three own markers and three foreign markers, in high-confidence Xenium cells and in the matching eligible atlas cells, and report the ratio.

Expected output: a table per type and gene with Xenium %, atlas % and the Xenium / atlas ratio. Own markers well below 1x point to a platform gap; foreign markers well above 1x point to spillover.

Note: the atlas has no lymphatic endothelial type, so lymphatic cells are compared with all atlas endothelial cells. The foreign markers tested are T-cell genes (or CD4-lineage genes for CD8 T); other contaminants are not tested here.

CHECKS = {   # manual type -> (matching atlas type, own markers, foreign markers)
    "Mast cell": ("Mast cell", ["CPA3", "MS4A2", "GATA2"], ["CD3E", "TRAC", "IL7R"]),
    "Lymphatic endothelial": ("Endothelial cell", ["PROX1", "MMRN1", "PDPN"], ["CD3E", "TRAC", "IL7R"]),
    "NK": ("NK cell", ["NKG7", "GNLY", "KLRD1"], ["CD3E", "TRAC", "CD8A"]),
    "CD8 T (GZMK+)": ("T cell CD8", ["CD8A", "NKG7", "GZMK"], ["CD4", "CCR7", "IL7R"]),
}
rows = []
for t, (atlas_type, own, foreign) in CHECKS.items():
    xc = ann[high[manual.loc[high, "cell_type"].to_numpy() == t]].layers["counts"]
    rc = panel_ref[(o["cell_type_major"] == atlas_type).to_numpy() & eligible.to_numpy()].X
    for kind, genes in [("own marker", own), ("foreign marker", foreign)]:
        for g in genes:
            rows.append({"manual type": t, "atlas type": atlas_type, "gene": g, "kind": kind,
                         "Xenium % detected": 100 * (xc[:, ann.var_names.get_loc(g)] > 0).mean(),
                         "atlas % detected": 100 * (rc[:, panel_ref.var_names.get_loc(g)] > 0).mean()})
detect = pd.DataFrame(rows).set_index(["manual type", "kind", "gene"])
detect["Xenium / atlas"] = detect["Xenium % detected"] / detect["atlas % detected"].clip(lower=1)
display(detect.drop(columns="atlas type").style.format({"Xenium % detected": "{:.0f}", "atlas % detected": "{:.0f}",
                                                        "Xenium / atlas": "{:.1f}x"}))

			Xenium % detected	atlas % detected	Xenium / atlas
manual type	kind	gene
Mast cell	own marker	CPA3	72	92	0.8x
		MS4A2	57	76	0.7x
		GATA2	33	70	0.5x
	foreign marker	CD3E	27	2	15.5x
		TRAC	43	3	13.0x
		IL7R	39	2	16.6x
Lymphatic endothelial	own marker	PROX1	45	9	5.1x
		MMRN1	53	21	2.5x
		PDPN	39	9	4.5x
	foreign marker	CD3E	57	1	56.6x
		TRAC	74	1	57.0x
		IL7R	78	20	4.0x
NK	own marker	NKG7	40	98	0.4x
		GNLY	56	89	0.6x
		KLRD1	54	85	0.6x
	foreign marker	CD3E	41	19	2.1x
		TRAC	44	9	4.8x
		CD8A	36	9	3.8x
CD8 T (GZMK+)	own marker	CD8A	37	51	0.7x
		NKG7	13	76	0.2x
		GZMK	41	38	1.1x
	foreign marker	CD4	25	4	6.6x
		CCR7	25	5	4.8x
		IL7R	67	34	2.0x

Interpretation.

• Lymphatic endothelial: spillover. Own markers are detected more often than in atlas endothelial cells (PROX1, MMRN1, PDPN at 2.5-5.1x, as expected for a lymphatic subset), but CD3E and TRAC appear 57x more often. The most common atlas calls are T-cell types (CD4 T 31%, Treg 17%), so the manual label is likely right and the outlines also contain parts of T cells.
• Mast cells: spillover plus weaker own markers. CD3E, TRAC and IL7R appear 13-17x more often than in atlas mast cells, while own markers are detected somewhat less often (CPA3, MS4A2, GATA2 at 0.5-0.8x). The most common wrong call is Plasma cell (24%), a contaminant this table does not test.
• CD8 T and NK: platform gap plus spillover. NKG7, a strong NK / CD8 marker in scRNA-seq (98% of atlas NK and 76% of atlas CD8 T cells), is detected in only 40% and 13% of the Xenium cells. CD8 T cells are most often called CD4 T (29%) and NK cells CD8 T (18%). CD8 T cells also carry CD4, CCR7 and IL7R 2-7x more often than atlas CD8 T cells, consistent with transcripts from neighbouring CD4 T cells and dendritic cells.

Lesson: a reference trained on scRNA-seq transfers well for types with robustly detected markers and poorly where a key marker is weak on the panel. Per-type agreement, read together with marker detection rates, shows which manual labels the reference supports and why the others cannot be confirmed.

Part E: Export

Purpose: make the atlas labels available to later analyses, for example to keep only cells where the manual and atlas labels agree.

Approach: join the manual labels with the atlas label, probability and lineage from both segmentations, indexed by 10x cell id, and save the per-type agreement table.

Expected output: the number of cells saved and the output folder, which also holds the trained model.

Note: overwrites atlas_labels_luca.csv.gz and agreement_by_cell_type.csv. Cells removed by Proseg QC have empty proseg_* columns.

export = (manual.join(atlas10.add_prefix("tenx_"), how="left")
          .join(atlasP.add_prefix("proseg_"), how="left"))
export["tenx_atlas_lineage"] = export["tenx_atlas_label"].map(ATLAS_LINEAGE)
export["proseg_atlas_lineage"] = export["proseg_atlas_label"].map(ATLAS_LINEAGE)
export.index.name = "cell_id"
path = OUT_DIR / "atlas_labels_luca.csv.gz"
export.to_csv(path)
agree10.to_csv(OUT_DIR / "agreement_by_cell_type.csv")
print(f"Saved {len(export):,} cells to {xu.display_path(path)}")
print(f"Saved the per-type agreement table and the trained model (luca_panel_genes_model.pkl) to {xu.display_path(OUT_DIR)}")

Saved 154,394 cells to ~/data/xenium_lung/processed/reference_mapping/atlas_labels_luca.csv.gz
Saved the per-type agreement table and the trained model (luca_panel_genes_model.pkl) to ~/data/xenium_lung/processed/reference_mapping

Conclusions

1. Lineages are largely supported. 92% of high-confidence cells receive the same lineage from an independent reference as from the manual annotation. That figure is carried by epithelium (99%), myeloid (91%) and stroma (90%); T / NK cells agree least (63%), with 15% called B / plasma and 15% myeloid. The ceiling for this panel on held-out atlas donors is 87% at the 23-type level.
2. Confirmed at cell-type level: smooth muscle, fibroblast, ciliated, proliferating tumour, alveolar macrophage, plasma, B and mregDC cells, and the S100A12+ group as neutrophil / monocyte.
3. Not confirmable from this panel: malignancy. The atlas places most tumour-labelled cells in an AT2 / club-like epithelial state. The tumour call remains an inference from expression programmes and spatial context, as noted in Step 6, and needs copy-number evidence to confirm.
4. Poor agreement has identifiable causes: spillover (lymphatic endothelial, mast) and a detection gap for NKG7 combined with spillover (CD8 T, NK). These are properties of the data and do not by themselves show that the labels are wrong.
5. Segmentation: on cells not selected for a clean marker (Step 6 mixed), Proseg counts let the reference recognise the T cell inside the outline, independent support for notebook 02's conclusion that Proseg removes spillover.

Limits: one reference and one classifier at the cell_type_major level; manual labels were made on 10x counts, which biases the segmentation comparison toward 10x; agreement is measured on high-confidence cells only.

Possible next steps: scANVI / scArches mapping (needs a GPU; better suited to the 51-type cell_type_tumor level), and copy-number inference on the Xenium cells to test the tumour call directly.

← Project home