The analysis dataset and its ground truth¶
Every notebook in this section analyses simulated data with a real, published tool and then asks
did it get the right answer. That question only has meaning because iscc knows the answer. This
notebook is where that answer lives: it grows a tumour through the public API, and shows what the
simulator records about it that no real experiment could tell you.
The other notebooks in this section reference this one rather than re-explaining it:
- Calling copy number from single-cell DNA with HMMcopy
- Tree reconstruction from single-cell DNA with SCITE
- Clone inference from bulk DNA with PyClone-VI
- Integrating scDNA and scRNA with clonealign
- Calling copy number from scRNA with Numbat
- Deconvolving Visium spots with RCTD
- Recovering progression constraints with TreeMHN
- Progression constraints from cross-sectional data with MHN
- Integrating a patient cohort with scDEF
Why the analysis notebooks load files instead of simulating¶
A benchmark should not be able to peek. The datasets those notebooks read contain the tool's input and nothing else, with the truth kept in a separate file that is opened only at scoring. Splitting generation from analysis also means an analysis notebook has no Python in it at all — which is why the four above can be R notebooks, running the R you would write yourself rather than a Python wrapper around it.
Generate them with:
python validation/make_analysis_data.py
In [1]:
Copied!
%matplotlib inline
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from iscc.tumor.models import GenotypeTumor
CONFIG = os.path.join(os.getcwd(), "example_config.yaml")
# A breast lesion on a ductal field: many small epithelial-ring glands in stroma, grown from ONE
# founder. It spends a long opening phase confined inside its founding acinus — that is where the
# clonal backbone forms — and only then spreads and invades.
tumor = GenotypeTumor(config=CONFIG, seed=2)
# Grow to a cancer-cell target rather than a fixed step count: the opening confined phase takes a
# variable number of generations, and stopping too early leaves a lesion of a handful of cells.
while tumor.get_cancer_size() < 20_000:
tumor.grow(n_steps=4, seed=2)
print(f"tumour: {tumor.get_cancer_size():,} cancer cells")
# `cell_data` is a BOUNDED, representative subsample of the tumour, not the whole thing — a cm-scale
# lesion has millions of cells and materialising every one of them is neither necessary nor
# affordable. So the counts below are the sample's, while get_cancer_size() above is the tumour's.
cd = tumor.cell_data
print(f"materialised sample: {len(cd['cell_type']):,} cells "
f"(of a {tumor.get_cancer_size():,}-cancer-cell tumour)\n")
print("ground-truth layers iscc records for every sampled cell:")
for k, v in cd.items():
shape = getattr(v, "shape", None)
print(f" {k:<16} {str(shape) if shape is not None else type(v).__name__}")
%matplotlib inline
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from iscc.tumor.models import GenotypeTumor
CONFIG = os.path.join(os.getcwd(), "example_config.yaml")
# A breast lesion on a ductal field: many small epithelial-ring glands in stroma, grown from ONE
# founder. It spends a long opening phase confined inside its founding acinus — that is where the
# clonal backbone forms — and only then spreads and invades.
tumor = GenotypeTumor(config=CONFIG, seed=2)
# Grow to a cancer-cell target rather than a fixed step count: the opening confined phase takes a
# variable number of generations, and stopping too early leaves a lesion of a handful of cells.
while tumor.get_cancer_size() < 20_000:
tumor.grow(n_steps=4, seed=2)
print(f"tumour: {tumor.get_cancer_size():,} cancer cells")
# `cell_data` is a BOUNDED, representative subsample of the tumour, not the whole thing — a cm-scale
# lesion has millions of cells and materialising every one of them is neither necessary nor
# affordable. So the counts below are the sample's, while get_cancer_size() above is the tumour's.
cd = tumor.cell_data
print(f"materialised sample: {len(cd['cell_type']):,} cells "
f"(of a {tumor.get_cancer_size():,}-cancer-cell tumour)\n")
print("ground-truth layers iscc records for every sampled cell:")
for k, v in cd.items():
shape = getattr(v, "shape", None)
print(f" {k:<16} {str(shape) if shape is not None else type(v).__name__}")
tumour: 21,514 cancer cells materialised sample: 8,178 cells (of a 21,514-cancer-cell tumour) ground-truth layers iscc records for every sampled cell: cell_evo (8178, 19) cell_snv (8178, 6000) cell_cnv (8178, 6000) cell_exp (8178, 6000) cell_crd (8178, 2) cell_type (8178, 1) cell_deme (8178, 1) cell_snv_p (8178, 6000) cell_snv_m (8178, 6000) cell_rna_vaf (8178, 6000) cell_wgd (8178, 1) cell_gland (8178, 1)
What "ground truth" means here¶
A real study infers all of the below. iscc assigned it, so each is exact.
The four columns a tool is scored against are the ones a real experiment cannot observe directly: which cells are malignant, which clone each belongs to, its copy number gene by gene, and — where there is a spatial assay — which cells ended up in which spot.
In [2]:
Copied!
gid = cd["cell_type"]["cell_id"].astype(str).values
types = np.array([tumor.genotypes[g].type if g in tumor.genotypes else "?" for g in gid])
is_cancer = types == "cancer"
cnv = cd["cell_cnv"].values
snv = cd["cell_snv"].values
print(f"cells: {len(types):,} malignant {is_cancer.sum():,} normal {(~is_cancer).sum():,} "
f"(purity {is_cancer.mean():.1%})")
print(f"genome: {cnv.shape[1]:,} genes over {tumor.n_segments} segments")
print("\ncell types present:", ", ".join(f"{t} {int((types==t).sum()):,}"
for t in sorted(set(types))))
# Segment boundaries: copy number is piecewise-constant along a segment, so the segments are the
# unit a CNA caller is really trying to find. Marking them makes the block structure legible.
bounds = np.cumsum(tumor.selection.segment_sizes)[:-1]
# Order the cells BY CLONE, so a clone's cells are contiguous and its copy-number profile reads as a
# block. Any real CNV heatmap clusters its rows; here the true clone label does it exactly.
mal = np.where(is_cancer)[0]
mal = mal[np.argsort(gid[mal], kind="stable")][:400] # the same cells, same order, in both panels
from matplotlib.colors import TwoSlopeNorm
fig, axes = plt.subplots(1, 3, figsize=(15, 4), width_ratios=[1, 1.3, 1.3])
crd = cd["cell_crd"]
axes[0].scatter(crd["col"], crd["row"], c=np.where(is_cancer, "#c1344e", "#c9c9c9"),
s=3, linewidths=0)
axes[0].set_title("where the malignant cells are"); axes[0].set_aspect("equal")
axes[0].set_xlabel("x (demes)"); axes[0].set_ylabel("y (demes)")
# Copy number is a DIVERGING quantity around the diploid baseline, so it gets a diverging map
# centred at 2 — and the scale never goes below zero, because negative copy number is not a thing.
im = axes[1].imshow(cnv[mal], aspect="auto", cmap="bwr", interpolation="nearest",
norm=TwoSlopeNorm(vcenter=2.0, vmin=0.0, vmax=max(4.0, cnv.max())))
axes[1].set_title("TRUE copy number (cells x genes)")
fig.colorbar(im, ax=axes[1], label="copy number")
# SNVs are sparse: most of the 6,000 genes are mutated in nobody, and plotting those columns hides
# the ones that matter. Restrict to genes mutated in at least one shown cell — what any mutation
# heatmap does — and keep them in genome order so the segments still line up.
mut = (snv[mal] > 0)
seen = np.where(mut.any(axis=0))[0]
im = axes[2].imshow(mut[:, seen].astype(float), aspect="auto", cmap="Greys",
interpolation="nearest", vmin=0, vmax=1)
axes[2].set_title("TRUE somatic mutations (cells x mutated genes)")
axes[2].set_xlabel(f"mutated gene ({len(seen):,} of {snv.shape[1]:,} carry any mutation)")
fig.colorbar(im, ax=axes[2], label="mutated", ticks=[0, 1])
for b in bounds:
axes[1].axvline(b - 0.5, color="0.25", lw=0.6, alpha=0.7)
axes[1].set_xlabel(f"gene ({tumor.n_segments} segments, boundaries marked)")
for ax in axes[1:]:
ax.set_ylabel(f"malignant cell (n={len(mal)}, ordered by clone)")
fig.tight_layout()
gid = cd["cell_type"]["cell_id"].astype(str).values
types = np.array([tumor.genotypes[g].type if g in tumor.genotypes else "?" for g in gid])
is_cancer = types == "cancer"
cnv = cd["cell_cnv"].values
snv = cd["cell_snv"].values
print(f"cells: {len(types):,} malignant {is_cancer.sum():,} normal {(~is_cancer).sum():,} "
f"(purity {is_cancer.mean():.1%})")
print(f"genome: {cnv.shape[1]:,} genes over {tumor.n_segments} segments")
print("\ncell types present:", ", ".join(f"{t} {int((types==t).sum()):,}"
for t in sorted(set(types))))
# Segment boundaries: copy number is piecewise-constant along a segment, so the segments are the
# unit a CNA caller is really trying to find. Marking them makes the block structure legible.
bounds = np.cumsum(tumor.selection.segment_sizes)[:-1]
# Order the cells BY CLONE, so a clone's cells are contiguous and its copy-number profile reads as a
# block. Any real CNV heatmap clusters its rows; here the true clone label does it exactly.
mal = np.where(is_cancer)[0]
mal = mal[np.argsort(gid[mal], kind="stable")][:400] # the same cells, same order, in both panels
from matplotlib.colors import TwoSlopeNorm
fig, axes = plt.subplots(1, 3, figsize=(15, 4), width_ratios=[1, 1.3, 1.3])
crd = cd["cell_crd"]
axes[0].scatter(crd["col"], crd["row"], c=np.where(is_cancer, "#c1344e", "#c9c9c9"),
s=3, linewidths=0)
axes[0].set_title("where the malignant cells are"); axes[0].set_aspect("equal")
axes[0].set_xlabel("x (demes)"); axes[0].set_ylabel("y (demes)")
# Copy number is a DIVERGING quantity around the diploid baseline, so it gets a diverging map
# centred at 2 — and the scale never goes below zero, because negative copy number is not a thing.
im = axes[1].imshow(cnv[mal], aspect="auto", cmap="bwr", interpolation="nearest",
norm=TwoSlopeNorm(vcenter=2.0, vmin=0.0, vmax=max(4.0, cnv.max())))
axes[1].set_title("TRUE copy number (cells x genes)")
fig.colorbar(im, ax=axes[1], label="copy number")
# SNVs are sparse: most of the 6,000 genes are mutated in nobody, and plotting those columns hides
# the ones that matter. Restrict to genes mutated in at least one shown cell — what any mutation
# heatmap does — and keep them in genome order so the segments still line up.
mut = (snv[mal] > 0)
seen = np.where(mut.any(axis=0))[0]
im = axes[2].imshow(mut[:, seen].astype(float), aspect="auto", cmap="Greys",
interpolation="nearest", vmin=0, vmax=1)
axes[2].set_title("TRUE somatic mutations (cells x mutated genes)")
axes[2].set_xlabel(f"mutated gene ({len(seen):,} of {snv.shape[1]:,} carry any mutation)")
fig.colorbar(im, ax=axes[2], label="mutated", ticks=[0, 1])
for b in bounds:
axes[1].axvline(b - 0.5, color="0.25", lw=0.6, alpha=0.7)
axes[1].set_xlabel(f"gene ({tumor.n_segments} segments, boundaries marked)")
for ax in axes[1:]:
ax.set_ylabel(f"malignant cell (n={len(mal)}, ordered by clone)")
fig.tight_layout()
cells: 8,178 malignant 605 normal 7,573 (purity 7.4%) genome: 6,000 genes over 12 segments cell types present: cancer 605, epithelial 269, stromal 7,304
The datasets, and what each tool is allowed to see¶
make_analysis_data.py turns a tumour like the one above into plain tables. The split is the point:
the tool reads the left column, and the right column is opened only once the tool has committed to an
answer.
| dataset | what the tool reads | ground truth, held back |
|---|---|---|
clonealign |
scRNA counts (cells × genes), clone copy number (genes × clones) | each cell's true clone |
numbat |
expression counts, a normal reference, phased allele counts, gene annotation | malignant vs normal, and each cell's clone |
rctd |
a Visium section (counts + coordinates) and a paired scRNA reference | the true cell-type proportion of every spot |
treemhn |
SCITE-reconstructed mutation trees, and a noisy single-cell genotype matrix per patient | the planted precedence network, and iscc's own true trees |
mhn |
the same kind of cohort read cross-sectionally, plus a no-constraint control arm | the planted precedence network |
dna |
single-cell read depth in bins, single-cell genotype calls, bulk read counts | the true copy number, the uncorrupted genotypes, and every mutation's true CCF and clonal identity |
cohort |
pooled scRNA counts from five patients, with a batch label | the shared program dictionary and every cell's true program activities |
Each dataset also ships a meta.json recording the shape of the problem and the baselines that
make a score readable — chance, and the majority-class rate. A method that beats chance but not
the majority baseline has not done anything useful, and the number alone will not tell you that.
In [3]:
Copied!
import json
DATA = os.path.join("..", "analysis_data")
if not os.path.isdir(DATA):
print("datasets not generated yet — run: python validation/make_analysis_data.py")
else:
man = json.load(open(os.path.join(DATA, "manifest.json")))
for name, m in man["datasets"].items():
keep = {k: v for k, v in m.items()
if k in ("n_cells", "n_genes", "n_clones", "n_spots", "n_patients",
"chance_baseline", "majority_baseline", "spot_purity_mean",
"frac_spots_multitype", "cells_per_spot_mean")}
pretty = ", ".join(f"{k}={v:,}" if isinstance(v, int) else f"{k}={v:.2f}"
for k, v in keep.items())
print(f" {name:<12} {pretty}")
import json
DATA = os.path.join("..", "analysis_data")
if not os.path.isdir(DATA):
print("datasets not generated yet — run: python validation/make_analysis_data.py")
else:
man = json.load(open(os.path.join(DATA, "manifest.json")))
for name, m in man["datasets"].items():
keep = {k: v for k, v in m.items()
if k in ("n_cells", "n_genes", "n_clones", "n_spots", "n_patients",
"chance_baseline", "majority_baseline", "spot_purity_mean",
"frac_spots_multitype", "cells_per_spot_mean")}
pretty = ", ".join(f"{k}={v:,}" if isinstance(v, int) else f"{k}={v:.2f}"
for k, v in keep.items())
print(f" {name:<12} {pretty}")
rctd n_spots=2,542, cells_per_spot_mean=4.89, spot_purity_mean=0.72, frac_spots_multitype=0.75 clonealign n_cells=1,501, n_genes=6,000, n_clones=4, majority_baseline=0.49, chance_baseline=0.25 numbat n_clones=4, n_cells=1,677, n_genes=6,000 treemhn n_patients=36 cohort n_patients=5, n_cells=1,500, n_genes=6,000 mhn n_patients=295 dna n_clones=4
The true clone phylogeny¶
The one piece of ground truth that no experiment recovers directly: which clone descends from which. Every tree-inference method estimates this from mutations, and every published tree is a hypothesis. Here it is recorded.
iscc.inference.clone_tree contracts the genotype genealogy onto the clone phylogeny — each node is
a block of genotypes sharing one driver combination, sized by how many cells it holds. This is the
answer the tree notebooks are scored against.
In [4]:
Copied!
from io import StringIO
from Bio import Phylo
from iscc.inference import clone_tree
from iscc.integrations.lineage import LineageTree
parents, sizes = clone_tree(tumor)
tot = max(sum(sizes.values()), 1)
print(f"true clone phylogeny: {len(sizes):,} clones over {tot:,} cells")
# Prune to the clones a real study could resolve. The full genealogy has one node per genotype and
# most are singletons, which no method recovers and no figure can show; a phylogeny is normally
# drawn over the clones that actually carry cells. `LineageTree` takes the node set and keeps the
# ancestors needed to connect them.
MIN_FRAC = 0.01
major = [g for g, n in sizes.items() if n / tot >= MIN_FRAC]
major_set = set(major)
# internal clones are only labelled when they are DOMINANT, not merely resolvable
DOMINANT_FRAC = 0.05
dominant = {g for g, n in sizes.items() if n / tot >= DOMINANT_FRAC}
print(f" clones above {MIN_FRAC:.0%} of the tumour: {len(major)}")
print(f" dominant clones (>= {DOMINANT_FRAC:.0%} of the tumour), largest first:")
for g in sorted(dominant, key=lambda g: -sizes[g]):
kind = "tip" if g not in set(parents.values()) else "ancestor of other clones"
print(f" {g:>6} {sizes[g]/tot:>5.0%} {kind}")
def clade_id(c):
"""The clone id of a clade, wherever Bio.Phylo put it.
Newick has no separate slot for an internal node NAME, so the parser files numeric internal
labels under `confidence` and leaves `name` empty -- which is why matching on `.name` alone
labels the tips and silently skips every internal clone.
"""
if c.name:
return str(c.name)
return str(int(c.confidence)) if c.confidence is not None else None
tree = LineageTree(dict(tumor.genotypes_parents), major)
phylo = Phylo.read(StringIO(tree.newick()), "newick")
phylo.ladderize()
# Newick + Bio.Phylo is the standard way to render a phylogeny, and it is what any tree-inference
# tool downstream will emit too — so the true tree and a reconstructed one are drawn the same way.
fig, ax = plt.subplots(figsize=(9.5, max(3.2, 0.30 * len(major))))
# The two biggest clones here (36% and 16%) are INTERNAL nodes -- ancestors that still carry most of
# the cells -- so tip labels alone would show a tree whose largest label is 3% while the text reports
# 36%. Tips get `label_func` (drawn to the right); internal clones get `branch_labels` (drawn above
# the branch), because putting both through `label_func` renders them on top of each other.
# show_confidence=False: Newick has no separate slot for an internal node NAME, so Bio.Phylo reads
# numeric internal labels as bootstrap confidences and prints them along every branch. Our clone ids
# are numeric, so without this the tree is covered in meaningless numbers.
Phylo.draw(phylo, axes=ax, do_show=False, show_confidence=False,
label_func=lambda c: (f" {clade_id(c)} ({sizes[clade_id(c)]/tot:.0%})"
if c.is_terminal() and clade_id(c) in major_set else ""),
)
ax.set_title(f"true clone phylogeny (clones >= {MIN_FRAC:.0%} of the tumour)")
ax.set_xlabel("branch length (divisions from the founder)")
fig.tight_layout()
from io import StringIO
from Bio import Phylo
from iscc.inference import clone_tree
from iscc.integrations.lineage import LineageTree
parents, sizes = clone_tree(tumor)
tot = max(sum(sizes.values()), 1)
print(f"true clone phylogeny: {len(sizes):,} clones over {tot:,} cells")
# Prune to the clones a real study could resolve. The full genealogy has one node per genotype and
# most are singletons, which no method recovers and no figure can show; a phylogeny is normally
# drawn over the clones that actually carry cells. `LineageTree` takes the node set and keeps the
# ancestors needed to connect them.
MIN_FRAC = 0.01
major = [g for g, n in sizes.items() if n / tot >= MIN_FRAC]
major_set = set(major)
# internal clones are only labelled when they are DOMINANT, not merely resolvable
DOMINANT_FRAC = 0.05
dominant = {g for g, n in sizes.items() if n / tot >= DOMINANT_FRAC}
print(f" clones above {MIN_FRAC:.0%} of the tumour: {len(major)}")
print(f" dominant clones (>= {DOMINANT_FRAC:.0%} of the tumour), largest first:")
for g in sorted(dominant, key=lambda g: -sizes[g]):
kind = "tip" if g not in set(parents.values()) else "ancestor of other clones"
print(f" {g:>6} {sizes[g]/tot:>5.0%} {kind}")
def clade_id(c):
"""The clone id of a clade, wherever Bio.Phylo put it.
Newick has no separate slot for an internal node NAME, so the parser files numeric internal
labels under `confidence` and leaves `name` empty -- which is why matching on `.name` alone
labels the tips and silently skips every internal clone.
"""
if c.name:
return str(c.name)
return str(int(c.confidence)) if c.confidence is not None else None
tree = LineageTree(dict(tumor.genotypes_parents), major)
phylo = Phylo.read(StringIO(tree.newick()), "newick")
phylo.ladderize()
# Newick + Bio.Phylo is the standard way to render a phylogeny, and it is what any tree-inference
# tool downstream will emit too — so the true tree and a reconstructed one are drawn the same way.
fig, ax = plt.subplots(figsize=(9.5, max(3.2, 0.30 * len(major))))
# The two biggest clones here (36% and 16%) are INTERNAL nodes -- ancestors that still carry most of
# the cells -- so tip labels alone would show a tree whose largest label is 3% while the text reports
# 36%. Tips get `label_func` (drawn to the right); internal clones get `branch_labels` (drawn above
# the branch), because putting both through `label_func` renders them on top of each other.
# show_confidence=False: Newick has no separate slot for an internal node NAME, so Bio.Phylo reads
# numeric internal labels as bootstrap confidences and prints them along every branch. Our clone ids
# are numeric, so without this the tree is covered in meaningless numbers.
Phylo.draw(phylo, axes=ax, do_show=False, show_confidence=False,
label_func=lambda c: (f" {clade_id(c)} ({sizes[clade_id(c)]/tot:.0%})"
if c.is_terminal() and clade_id(c) in major_set else ""),
)
ax.set_title(f"true clone phylogeny (clones >= {MIN_FRAC:.0%} of the tumour)")
ax.set_xlabel("branch length (divisions from the founder)")
fig.tight_layout()
true clone phylogeny: 987 clones over 21,514 cells
clones above 1% of the tumour: 12
dominant clones (>= 5% of the tumour), largest first:
741 36% ancestor of other clones
932 16% ancestor of other clones
821 5% ancestor of other clones
The truth each tool is scored against¶
Everything below is loaded from the truth_* files of the actual datasets the analysis notebooks
read. This is where the answer key is drawn. Those notebooks show what their tool produced, using
the tool's own plotting; the panels to hold that output against are here, so the comparison is a
comparison of two figures rather than a figure the benchmark drew for you.
The copy number itself is not redrawn per tool — it is the cells x genes panel above, which is the ground truth all of them are scored against.
In [5]:
Copied!
import numpy as np
DATA = os.path.join("..", "analysis_data")
# --- clonealign: which clone each cell really belongs to, and what copy number that clone carries
ca_truth = pd.read_csv(os.path.join(DATA, "clonealign", "truth.csv"))
ca_meta = json.load(open(os.path.join(DATA, "clonealign", "meta.json")))
cn = np.array(ca_meta["cn_consensus"])
# --- numbat: which cells are malignant, and the consensus copy number of each malignant clone
nb_truth = pd.read_csv(os.path.join(DATA, "numbat", "truth.csv"))
nb_cn = np.loadtxt(os.path.join(DATA, "numbat", "truth_consensus_cn.csv"), delimiter=",", ndmin=2)
print("clonealign — true clone sizes:", ca_truth["true_clone"].value_counts().sort_index().to_dict())
print(f" majority baseline {ca_meta['majority_baseline']:.2f}, chance {ca_meta['chance_baseline']:.2f}")
print(f"numbat — {int(nb_truth['is_malignant'].sum()):,} malignant of {len(nb_truth):,} cells "
f"(purity {nb_truth['is_malignant'].mean():.1%}), {nb_cn.shape[0]} malignant clones")
import numpy as np
DATA = os.path.join("..", "analysis_data")
# --- clonealign: which clone each cell really belongs to, and what copy number that clone carries
ca_truth = pd.read_csv(os.path.join(DATA, "clonealign", "truth.csv"))
ca_meta = json.load(open(os.path.join(DATA, "clonealign", "meta.json")))
cn = np.array(ca_meta["cn_consensus"])
# --- numbat: which cells are malignant, and the consensus copy number of each malignant clone
nb_truth = pd.read_csv(os.path.join(DATA, "numbat", "truth.csv"))
nb_cn = np.loadtxt(os.path.join(DATA, "numbat", "truth_consensus_cn.csv"), delimiter=",", ndmin=2)
print("clonealign — true clone sizes:", ca_truth["true_clone"].value_counts().sort_index().to_dict())
print(f" majority baseline {ca_meta['majority_baseline']:.2f}, chance {ca_meta['chance_baseline']:.2f}")
print(f"numbat — {int(nb_truth['is_malignant'].sum()):,} malignant of {len(nb_truth):,} cells "
f"(purity {nb_truth['is_malignant'].mean():.1%}), {nb_cn.shape[0]} malignant clones")
clonealign — true clone sizes: {0: 575, 1: 730, 2: 87, 3: 109}
majority baseline 0.49, chance 0.25
numbat — 1,527 malignant of 1,677 cells (purity 91.1%), 3 malignant clones
Spatial: what is actually in each Visium spot¶
iscc placed the cells, so it knows each spot's composition exactly — no reference-based estimate
standing in for truth. This is the panel to compare RCTD's weight maps against.
In [6]:
Copied!
sp_truth = pd.read_csv(os.path.join(DATA, "rctd", "truth_type_composition.csv"), index_col=0)
sp_coords = pd.read_csv(os.path.join(DATA, "rctd", "sp_coords.csv"), index_col=0)
types = [c for c in sp_truth.columns if sp_truth[c].sum() > 0]
print(f"{len(sp_truth):,} spots; mean composition:",
{c: round(float(sp_truth[c].mean()), 3) for c in types})
fig, axes = plt.subplots(1, len(types), figsize=(3.4 * len(types), 3.4))
for ax, k in zip(np.atleast_1d(axes), types):
s = ax.scatter(sp_coords["x"], sp_coords["y"], c=sp_truth.loc[sp_coords.index, k],
cmap="Purples", vmin=0, vmax=1, s=4, linewidths=0) # fractions -> Purples
ax.set_title(f"true {k} fraction"); ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
fig.colorbar(s, ax=ax, fraction=0.046)
fig.tight_layout()
sp_truth = pd.read_csv(os.path.join(DATA, "rctd", "truth_type_composition.csv"), index_col=0)
sp_coords = pd.read_csv(os.path.join(DATA, "rctd", "sp_coords.csv"), index_col=0)
types = [c for c in sp_truth.columns if sp_truth[c].sum() > 0]
print(f"{len(sp_truth):,} spots; mean composition:",
{c: round(float(sp_truth[c].mean()), 3) for c in types})
fig, axes = plt.subplots(1, len(types), figsize=(3.4 * len(types), 3.4))
for ax, k in zip(np.atleast_1d(axes), types):
s = ax.scatter(sp_coords["x"], sp_coords["y"], c=sp_truth.loc[sp_coords.index, k],
cmap="Purples", vmin=0, vmax=1, s=4, linewidths=0) # fractions -> Purples
ax.set_title(f"true {k} fraction"); ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
fig.colorbar(s, ax=ax, fraction=0.046)
fig.tight_layout()
2,542 spots; mean composition: {'cancer': 0.21, 'epithelial': 0.004, 'stromal': 0.542, 'immune': 0.243}
The gene programs, and which genes are in them¶
The cohort dataset plants one shared program dictionary that every patient draws from, while each patient's tumour evolves privately. Two things are recorded: the loading matrix (which genes each program spans, and how strongly) and the per-cell activities (how much each cell used each program). A method that claims to have recovered a program can be held against both.
In [7]:
Copied!
loading = pd.read_csv(os.path.join(DATA, "cohort", "truth_loading.csv.gz"), index_col=0)
activity = pd.read_csv(os.path.join(DATA, "cohort", "truth_activity.csv.gz"), index_col=0)
obs = pd.read_csv(os.path.join(DATA, "cohort", "obs.csv"), index_col=0)
print(f"{loading.shape[0]} programs over {loading.shape[1]} genes; "
f"{(loading > 0).sum(1).to_dict()} genes carried per program\n")
for prog in loading.index: # the genes that DEFINE each program
top = loading.loc[prog].sort_values(ascending=False).head(6)
print(f" {prog}: " + ", ".join(f"{g} {w:.2f}" for g, w in top.items()))
order = np.argsort(obs.loc[activity.index, "patient"].values, kind="stable")
fig, axes = plt.subplots(1, 2, figsize=(12, 3.6), width_ratios=[2, 1])
im = axes[0].imshow(loading.values, aspect="auto", cmap="cividis") # programs -> cividis
axes[0].set_yticks(range(loading.shape[0])); axes[0].set_yticklabels(loading.index)
axes[0].set_xlabel("gene"); axes[0].set_title("planted loading matrix (programs x genes)")
fig.colorbar(im, ax=axes[0], label="loading")
im = axes[1].imshow(activity.values[order].T, aspect="auto", cmap="cividis")
axes[1].set_yticks(range(activity.shape[1])); axes[1].set_yticklabels(activity.columns)
axes[1].set_xlabel("cell (grouped by patient)"); axes[1].set_title("true per-cell activities")
fig.colorbar(im, ax=axes[1], label="activity")
fig.tight_layout()
loading = pd.read_csv(os.path.join(DATA, "cohort", "truth_loading.csv.gz"), index_col=0)
activity = pd.read_csv(os.path.join(DATA, "cohort", "truth_activity.csv.gz"), index_col=0)
obs = pd.read_csv(os.path.join(DATA, "cohort", "obs.csv"), index_col=0)
print(f"{loading.shape[0]} programs over {loading.shape[1]} genes; "
f"{(loading > 0).sum(1).to_dict()} genes carried per program\n")
for prog in loading.index: # the genes that DEFINE each program
top = loading.loc[prog].sort_values(ascending=False).head(6)
print(f" {prog}: " + ", ".join(f"{g} {w:.2f}" for g, w in top.items()))
order = np.argsort(obs.loc[activity.index, "patient"].values, kind="stable")
fig, axes = plt.subplots(1, 2, figsize=(12, 3.6), width_ratios=[2, 1])
im = axes[0].imshow(loading.values, aspect="auto", cmap="cividis") # programs -> cividis
axes[0].set_yticks(range(loading.shape[0])); axes[0].set_yticklabels(loading.index)
axes[0].set_xlabel("gene"); axes[0].set_title("planted loading matrix (programs x genes)")
fig.colorbar(im, ax=axes[0], label="loading")
im = axes[1].imshow(activity.values[order].T, aspect="auto", cmap="cividis")
axes[1].set_yticks(range(activity.shape[1])); axes[1].set_yticklabels(activity.columns)
axes[1].set_xlabel("cell (grouped by patient)"); axes[1].set_title("true per-cell activities")
fig.colorbar(im, ax=axes[1], label="activity")
fig.tight_layout()
6 programs over 6000 genes; {'program0': 25, 'program1': 25, 'program2': 25, 'program3': 25, 'program4': 25, 'program5': 25} genes carried per program
program0: G_2_341 2.23, G_8_161 2.08, G_11_273 1.48, G_11_214 1.25, G_2_407 1.12, G_4_217 0.98
program1: G_9_350 2.40, G_4_402 2.19, G_7_188 1.90, G_7_195 1.80, G_10_468 1.62, G_10_138 1.52
program2: G_0_370 2.09, G_11_137 1.70, G_1_322 1.54, G_9_243 1.36, G_9_161 1.20, G_4_68 1.17
program3: G_9_87 3.39, G_4_68 2.07, G_10_424 2.03, G_9_468 1.61, G_2_430 1.01, G_4_408 0.85
program4: G_5_143 3.21, G_7_368 2.01, G_8_280 1.66, G_0_401 1.66, G_7_217 1.33, G_1_427 1.00
program5: G_6_213 6.08, G_1_483 4.74, G_0_343 4.48, G_2_337 3.70, G_11_21 2.64, G_0_470 2.16
Progression: the planted network, and the trees it produced¶
The progression cohort is grown under a DAG stating which events must precede which, gated on accessibility so the constraint acts on the mutation process itself and therefore shows up in tree topology. Two ground truths come out of it: the DAG, and each patient's true mutation tree.
The tree matters because the trees TreeMHN reads are not these — they are reconstructed by SCITE from noisy single-cell genotypes, exactly as a real study would. These are what that reconstruction is scored against.
In [8]:
Copied!
tm = os.path.join(DATA, "treemhn")
truth_net = json.load(open(os.path.join(tm, "truth_network.json")))
dag = truth_net["true_dag_edges"]
X = pd.read_csv(os.path.join(tm, "X_presence.csv"), index_col=0)
F = pd.read_csv(os.path.join(tm, "X_cellfraction.csv"), index_col=0)
ev = list(X.columns)
print(f"planted precedence constraints ({truth_net['dependency_params']['gating_mode']}-gated):")
for p, c in dag:
print(f" {ev[p]} must precede {ev[c]}")
# The binary column is at min_freq=0 -- ANY cell carrying the event counts. The continuous
# cancer-cell-fraction matrix is the same cohort without that floor collapsed in.
print("\nevent presence at min_freq=0:", X.mean(0).round(2).to_dict())
print("median cancer-cell fraction: ", F.median(0).round(3).to_dict())
truth_trees = pd.read_csv(os.path.join(tm, "truth_trees.csv"))
print(f"\ntrue mutation trees: {truth_trees['Patient_ID'].nunique()} patients, "
f"{len(truth_trees)} nodes total")
fig, axes = plt.subplots(1, 2, figsize=(11, 3.4))
axes[0].bar(ev, X.mean(0).values, color="#7a5195")
axes[0].plot(ev, F.median(0).values, "o-", color="#c1344e", label="median cell fraction")
axes[0].set_ylim(0, 1.05); axes[0].set_ylabel("fraction of patients / of cells")
axes[0].set_title("event presence vs cancer-cell fraction"); axes[0].legend(fontsize=8)
sizes = truth_trees.groupby("Patient_ID").size() - 1 # nodes below the root = events acquired
axes[1].hist(sizes, bins=range(0, int(sizes.max()) + 2), color="#7a5195", align="left")
axes[1].set_xlabel("events on the patient's true tree"); axes[1].set_ylabel("patients")
axes[1].set_title("size of the true mutation trees")
fig.tight_layout()
# The SAME kind of constraint, in the cohort built for the cross-sectional view. It is a separate
# cohort because the two observables need different densities: the trees above need events common
# enough to leave an ordering trace, which fixes the parents in every patient and leaves a
# presence/absence method nothing to estimate.
mhn_dir = os.path.join(DATA, "mhn")
if os.path.isdir(mhn_dir):
Xm = pd.read_csv(os.path.join(mhn_dir, "X_presence.csv"), index_col=0)
Xc = pd.read_csv(os.path.join(mhn_dir, "X_presence_control.csv"), index_col=0)
mm = json.load(open(os.path.join(mhn_dir, "meta.json")))
print(f"\nMHN cohort: {mm['n_patients']} patients ({mm['n_control_patients']} control), "
f"event_size {mm['event_size']}, {mm['grow_steps']} growth steps")
print(" presence:", Xm.mean(axis=0).round(2).to_dict(),
" control:", Xc.mean(axis=0).round(2).to_dict())
tm = os.path.join(DATA, "treemhn")
truth_net = json.load(open(os.path.join(tm, "truth_network.json")))
dag = truth_net["true_dag_edges"]
X = pd.read_csv(os.path.join(tm, "X_presence.csv"), index_col=0)
F = pd.read_csv(os.path.join(tm, "X_cellfraction.csv"), index_col=0)
ev = list(X.columns)
print(f"planted precedence constraints ({truth_net['dependency_params']['gating_mode']}-gated):")
for p, c in dag:
print(f" {ev[p]} must precede {ev[c]}")
# The binary column is at min_freq=0 -- ANY cell carrying the event counts. The continuous
# cancer-cell-fraction matrix is the same cohort without that floor collapsed in.
print("\nevent presence at min_freq=0:", X.mean(0).round(2).to_dict())
print("median cancer-cell fraction: ", F.median(0).round(3).to_dict())
truth_trees = pd.read_csv(os.path.join(tm, "truth_trees.csv"))
print(f"\ntrue mutation trees: {truth_trees['Patient_ID'].nunique()} patients, "
f"{len(truth_trees)} nodes total")
fig, axes = plt.subplots(1, 2, figsize=(11, 3.4))
axes[0].bar(ev, X.mean(0).values, color="#7a5195")
axes[0].plot(ev, F.median(0).values, "o-", color="#c1344e", label="median cell fraction")
axes[0].set_ylim(0, 1.05); axes[0].set_ylabel("fraction of patients / of cells")
axes[0].set_title("event presence vs cancer-cell fraction"); axes[0].legend(fontsize=8)
sizes = truth_trees.groupby("Patient_ID").size() - 1 # nodes below the root = events acquired
axes[1].hist(sizes, bins=range(0, int(sizes.max()) + 2), color="#7a5195", align="left")
axes[1].set_xlabel("events on the patient's true tree"); axes[1].set_ylabel("patients")
axes[1].set_title("size of the true mutation trees")
fig.tight_layout()
# The SAME kind of constraint, in the cohort built for the cross-sectional view. It is a separate
# cohort because the two observables need different densities: the trees above need events common
# enough to leave an ordering trace, which fixes the parents in every patient and leaves a
# presence/absence method nothing to estimate.
mhn_dir = os.path.join(DATA, "mhn")
if os.path.isdir(mhn_dir):
Xm = pd.read_csv(os.path.join(mhn_dir, "X_presence.csv"), index_col=0)
Xc = pd.read_csv(os.path.join(mhn_dir, "X_presence_control.csv"), index_col=0)
mm = json.load(open(os.path.join(mhn_dir, "meta.json")))
print(f"\nMHN cohort: {mm['n_patients']} patients ({mm['n_control_patients']} control), "
f"event_size {mm['event_size']}, {mm['grow_steps']} growth steps")
print(" presence:", Xm.mean(axis=0).round(2).to_dict(),
" control:", Xc.mean(axis=0).round(2).to_dict())
planted precedence constraints (accessibility-gated):
E3 must precede E1
E0 must precede E2
event presence at min_freq=0: {'E0': 1.0, 'E1': 0.44, 'E2': 0.33, 'E3': 1.0}
median cancer-cell fraction: {'E0': 0.07, 'E1': 0.0, 'E2': 0.0, 'E3': 0.134}
true mutation trees: 36 patients, 150 nodes total
MHN cohort: 295 patients (295 control), event_size 7, 150 growth steps
presence: {'E0': 0.78, 'E1': 0.06, 'E2': 0.25, 'E3': 0.83} control: {'E0': 0.8, 'E1': 0.63, 'E2': 0.94, 'E3': 0.64}
The DNA chain: copy number, genotypes and clonal identity¶
Three notebooks read one tumour — HMMcopy calls its copy number, SCITE reconstructs a tree from its single-cell genotypes, and PyClone-VI clusters its bulk read counts. Here is what each of them is held against.
In [9]:
Copied!
dna_dir = os.path.join(DATA, "dna")
dmeta = json.load(open(os.path.join(dna_dir, "meta.json")))
print(f"one tumour, {dmeta['n_cancer_cells']:,} cancer cells in {dmeta['n_clones']} clones; "
f"{dmeta['hmmcopy']['n_cells']} cells sequenced single-cell over "
f"{dmeta['hmmcopy']['n_bins']} bins, {dmeta['pyclonevi']['n_mutations']} mutations in bulk "
f"at purity {dmeta['pyclonevi']['purity']:.2f}\n")
# what SCITE is scored on: the genotypes BEFORE the assay corrupted them
gtruth = pd.read_csv(os.path.join(dna_dir, "scite", "truth_genotypes.csv"), index_col=0)
gcall = pd.read_csv(os.path.join(dna_dir, "scite", "sc_mutations.csv"), index_col=0)
fn = ((gtruth.values == 1) & (gcall.values == 0)).sum() / max((gtruth.values == 1).sum(), 1)
fp = ((gtruth.values == 0) & (gcall.values == 1)).sum() / max((gtruth.values == 0).sum(), 1)
print(f"single-cell genotypes: {gtruth.shape[0]} mutations x {gtruth.shape[1]} cells; "
f"the assay lost {fn:.1%} of real calls and invented {fp:.2%}")
# what PyClone-VI is scored on: each mutation's true CCF and which clones carry it
ptruth = pd.read_csv(os.path.join(dna_dir, "pyclonevi", "truth.csv"))
print(f"\nbulk mutations by true clonal identity (which clones carry it):")
print(ptruth.groupby("carrier_clones").agg(n=("mutation_id", "size"),
mean_ccf=("true_ccf", "mean")).round(3).to_string())
fig, axes = plt.subplots(1, 2, figsize=(11, 3.4))
order = np.argsort(-gtruth.values.sum(1))
axes[0].imshow(gtruth.values[order][:, np.argsort(gtruth.values.sum(0), kind="stable")],
aspect="auto", cmap="Greys", interpolation="nearest")
axes[0].set_title("true single-cell genotypes (what SCITE must recover)")
axes[0].set_xlabel("cell"); axes[0].set_ylabel("mutation")
axes[1].hist(ptruth["true_ccf"], bins=40, color="#7a5195")
axes[1].set_xlabel("true cancer-cell fraction"); axes[1].set_ylabel("mutations")
axes[1].set_title("true CCF (what PyClone-VI must recover)")
fig.tight_layout()
dna_dir = os.path.join(DATA, "dna")
dmeta = json.load(open(os.path.join(dna_dir, "meta.json")))
print(f"one tumour, {dmeta['n_cancer_cells']:,} cancer cells in {dmeta['n_clones']} clones; "
f"{dmeta['hmmcopy']['n_cells']} cells sequenced single-cell over "
f"{dmeta['hmmcopy']['n_bins']} bins, {dmeta['pyclonevi']['n_mutations']} mutations in bulk "
f"at purity {dmeta['pyclonevi']['purity']:.2f}\n")
# what SCITE is scored on: the genotypes BEFORE the assay corrupted them
gtruth = pd.read_csv(os.path.join(dna_dir, "scite", "truth_genotypes.csv"), index_col=0)
gcall = pd.read_csv(os.path.join(dna_dir, "scite", "sc_mutations.csv"), index_col=0)
fn = ((gtruth.values == 1) & (gcall.values == 0)).sum() / max((gtruth.values == 1).sum(), 1)
fp = ((gtruth.values == 0) & (gcall.values == 1)).sum() / max((gtruth.values == 0).sum(), 1)
print(f"single-cell genotypes: {gtruth.shape[0]} mutations x {gtruth.shape[1]} cells; "
f"the assay lost {fn:.1%} of real calls and invented {fp:.2%}")
# what PyClone-VI is scored on: each mutation's true CCF and which clones carry it
ptruth = pd.read_csv(os.path.join(dna_dir, "pyclonevi", "truth.csv"))
print(f"\nbulk mutations by true clonal identity (which clones carry it):")
print(ptruth.groupby("carrier_clones").agg(n=("mutation_id", "size"),
mean_ccf=("true_ccf", "mean")).round(3).to_string())
fig, axes = plt.subplots(1, 2, figsize=(11, 3.4))
order = np.argsort(-gtruth.values.sum(1))
axes[0].imshow(gtruth.values[order][:, np.argsort(gtruth.values.sum(0), kind="stable")],
aspect="auto", cmap="Greys", interpolation="nearest")
axes[0].set_title("true single-cell genotypes (what SCITE must recover)")
axes[0].set_xlabel("cell"); axes[0].set_ylabel("mutation")
axes[1].hist(ptruth["true_ccf"], bins=40, color="#7a5195")
axes[1].set_xlabel("true cancer-cell fraction"); axes[1].set_ylabel("mutations")
axes[1].set_title("true CCF (what PyClone-VI must recover)")
fig.tight_layout()
one tumour, 1,501 cancer cells in 4 clones; 260 cells sequenced single-cell over 240 bins, 300 mutations in bulk at purity 0.50
single-cell genotypes: 20 mutations x 200 cells; the assay lost 10.8% of real calls and invented 0.38%
bulk mutations by true clonal identity (which clones carry it):
n mean_ccf
carrier_clones
c0000 242 0.011
c0001 2 0.071
c0010 16 0.084
c1110 8 0.921
c1111 32 0.999
Reading a score honestly¶
Two habits carried through every notebook in this section, both learned from results that looked fine and were not:
Report the baseline, not just the score. A clone-assignment accuracy of 0.67 sounds reasonable until you notice the largest clone holds 70% of the cells. Every scoring cell prints chance and majority beside the result.
Report a threshold-free measure too. Accuracy depends on where a cut-off falls and on how imbalanced the classes are; AUC and ARI do not. A method can rank cells correctly — real, usable signal — and still lose on hard assignment. Reporting only the first number would call that a failure.