The shared tumour — one multi-focal DCIS→IDC lesion behind every science notebook¶

Every science notebook in this folder (Integrating scDNA and scRNA with clonealign, Tree reconstruction from single-cell DNA with SCITE, Deconvolving Visium spots with RCTD, compartment_selection_confound, Integrating a patient cohort with scDEF) is built on this one spatially-structured tumour, grown deterministically by base_sim.grow_base_tumor(). The other tutorials (01_pipeline_walkthrough, the assay_* notebooks) cover the mechanics; these show the science — what each modality lets you recover, and how iscc's ground truth lets you score it honestly.

Rather than ship a giant CSV, we re-grow the tumour on demand (~15–30 s, byte-for-byte reproducible), so every notebook shares the identical substrate.

The model: a multi-focal breast DCIS → IDC on a ductal field¶

The substrate is a ductal field (an island model): a scatter of small

epithelial-ring glands in moderate-density stroma. A single cancer founder starts in one gland's

lumen and colonises the others through low-rate cross-gland (island) dispersal — so a section shows

several clonally related foci (multi-focal DCIS), exactly as a real duct system's out-of-plane

branches would produce. Two normal compartments then impose selection

the gland wall and the hostile stroma each add a local death hazard a clone must evolve a heritable trait to survive — breach (cross the wall) and stromal_survival (traverse the stroma). The lesion stays confined (in situ, DCIS) until a subclone evolves those escape traits and invades the stroma (IDC).

feature (notebooks/base_sim.py)	why it matters downstream
ductal field (n_glands=4, glands in stroma)	every notebook analyses a mixture — malignant + epithelial + stromal — the normals are the CNA-caller reference, the deconvolution background and the balanced-BAF anchor
one founder + cross-gland dispersal (cross_gland_kappa)	multi-focal, clonally-related foci → a clonal phylogeny to reconstruct (SCITE)
compartment selection (breach, stromal_survival; breach-gated invasion)	DCIS→IDC: confined foci that invade only once escape traits evolve — sequenceable, recoverable traits (compartment_selection_confound)
genetic-vs-niche emt confound (breach → emt + epithelial → emt)	the invasive program is driven by BOTH genotype and location — a confound no real dataset can resolve (compartment_selection_confound)
CINner selection (prop_driver)	a real fitness gradient of clones, not a saturated sheet
whole-genome duplication (wgd_rate=0.05)	a PCAWG-like fraction of malignant cells double their genome (the allele layer)
allele-specific expression	per-homolog dosage → cell_rna_baf (the allelic-imbalance signal)
gene programs	a known loading matrix to recover from scRNA (Integrating a patient cohort with scDEF)

It grows to ≥10 000 cancer cells. Everything below is ground truth iscc knows and no real dataset does — we only ever use it to score and colour.

In [1]:

%matplotlib inline
import os, sys, time
sys.path.insert(0, os.path.abspath("."))   # so `import base_sim` works when executed headless
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import iscc
import base_sim as B
from iscc.tumor.models import GenotypeTumor
from iscc.tumor import viz
print("iscc from:", iscc.__file__)

%matplotlib inline import os, sys, time sys.path.insert(0, os.path.abspath(".")) # so `import base_sim` works when executed headless import numpy as np import pandas as pd import matplotlib.pyplot as plt import iscc import base_sim as B from iscc.tumor.models import GenotypeTumor from iscc.tumor import viz print("iscc from:", iscc.__file__)

iscc from: /Users/pedroferreira/projects/iscc/repo/src/iscc/__init__.py

1. Grow the shared tumour¶

grow_base_tumor() grows a single seed to ≥10 000 cancer cells on the ductal field and materialises cell_data. We do not filter to cancer: the sample is malignant cells plus their microenvironment (epithelial + stromal), exactly as a real dissociation would be. Note the cell_gland label (which gland, or −1 for stroma) and the compartment traits (breach, stromal_survival) — ground truth the downstream notebooks lean on.

In [2]:

t0 = time.time()
t = B.grow_base_tumor()          # seed 3: grows a cm-scale IDC, materialises a max_cells biopsy
cd = t.cell_data
print(f"grown in {time.time()-t0:.0f}s")

ty = B.cell_types(t)
cancer = ty == "cancer"
n_cancer = int(cancer.sum())
is_wgd = cd["cell_wgd"]["is_wgd"].values
living_clones = B.cancer_clones(t)
colonised = B.glands_colonised(t)

print(f"total cells          : {len(ty):>7,d}")
for name in ("cancer", "epithelial", "stromal"):
    print(f"  {name:<11s}        : {int((ty==name).sum()):>7,d}")
print(f"tumour purity        : {cancer.mean():>7.1%}")
print(f"WGD fraction (malig) : {is_wgd[cancer].mean():>7.1%}")
print(f"living cancer clones : {living_clones:>7,d}")
print(f"glands colonised     : {len(colonised):>4d} / {t.n_glands}  (from ONE founder)")
print(f"cancer in stroma (IDC): {B.stroma_cancer_pct(t):>6.1f}%   ({100-B.stroma_cancer_pct(t):.0f}% still confined in glands, DCIS)")
print(f"breach genes {len(t.selection.get_breach())} · stromal_survival genes {len(t.selection.get_stromal_survival())} (both sequenceable)")
assert t.get_cancer_size() >= 50000, f"expected a large tumour, got {t.get_cancer_size()}"
print(f"\nOK — {t.get_cancer_size():,} cancer cells grown; {n_cancer:,} in the materialised biopsy")

t0 = time.time() t = B.grow_base_tumor() # seed 3: grows a cm-scale IDC, materialises a max_cells biopsy cd = t.cell_data print(f"grown in {time.time()-t0:.0f}s") ty = B.cell_types(t) cancer = ty == "cancer" n_cancer = int(cancer.sum()) is_wgd = cd["cell_wgd"]["is_wgd"].values living_clones = B.cancer_clones(t) colonised = B.glands_colonised(t) print(f"total cells : {len(ty):>7,d}") for name in ("cancer", "epithelial", "stromal"): print(f" {name:<11s} : {int((ty==name).sum()):>7,d}") print(f"tumour purity : {cancer.mean():>7.1%}") print(f"WGD fraction (malig) : {is_wgd[cancer].mean():>7.1%}") print(f"living cancer clones : {living_clones:>7,d}") print(f"glands colonised : {len(colonised):>4d} / {t.n_glands} (from ONE founder)") print(f"cancer in stroma (IDC): {B.stroma_cancer_pct(t):>6.1f}% ({100-B.stroma_cancer_pct(t):.0f}% still confined in glands, DCIS)") print(f"breach genes {len(t.selection.get_breach())} · stromal_survival genes {len(t.selection.get_stromal_survival())} (both sequenceable)") assert t.get_cancer_size() >= 50000, f"expected a large tumour, got {t.get_cancer_size()}" print(f"\nOK — {t.get_cancer_size():,} cancer cells grown; {n_cancer:,} in the materialised biopsy")

grown in 97s
total cells          :   7,838
  cancer             :   1,867
  epithelial         :     205
  stromal            :   5,766
tumour purity        :   23.8%
WGD fraction (malig) :  100.0%
living cancer clones :  20,124
glands colonised     :    8 / 8  (from ONE founder)
cancer in stroma (IDC):   74.9%   (25% still confined in glands, DCIS)
breach genes 117 · stromal_survival genes 121 (both sequenceable)

OK — 80,347 cancer cells grown; 1,867 in the materialised biopsy

2. Watch it grow — the multi-focal DCIS→IDC time-series¶

This is the picture the whole suite rests on. We grow a fresh copy of the same simulation (seed 3, identical parameters) and snapshot it at five milestones, from seeding to the final lesion.

• Top row — the tissue (cell resolution). Each deme is drawn as a block of its individual cells (a 2-D section, ~40 % of the deme's 3-D depth): green epithelial gland walls, pink stroma, red cancer. Cancer appears in one gland, seeds the others (multi-focal), fills the lumens (DCIS), then breaks through the walls into the stroma (IDC).
• Bottom row — cancer coloured by gland. Each colonised gland's cancer takes that gland's colour; cancer that has invaded the stroma is black. One founder (a single colour) lights up gland after gland, then the black stromal mass grows — the DCIS→IDC transition, made visual.

In [3]:

import matplotlib.patches as mpatches
MILESTONES = [300, 3000, 20000, 50000, 80000]   # cancer-cell counts to snapshot (seeding -> confluent IDC), cm-scale
tg = B.new_tumor(max_cells=8000)   # ungrown, cm-safe realistic ductal field; stepped below to build the movie
snaps, mi = [], 0
while mi < len(MILESTONES):
    ncan = B.n_cancer(tg)
    if ncan >= MILESTONES[mi]:
        tg.make_cell_data()
        img, _ = B.expanded_tissue_rgb(tg, section_frac=0.4)   # fast tissue section (no per-clone Muller colormap)
        snaps.append(dict(step=tg.step, ncan=ncan, ngl=len(B.glands_colonised(tg)),
                          stroma=B.stroma_cancer_pct(tg), img=img, gland_grid=B.cancer_gland_grid(tg)))
        mi += 1
        continue
    tg.grow(n_steps=(2 if ncan < 4000 else 5), seed=B.BASE_SEED)   # fine early so the DCIS phase is caught

ncol = len(snaps)
fig, axes = plt.subplots(2, ncol, figsize=(3.0 * ncol, 6.6), constrained_layout=True)
gland_cmap = plt.get_cmap("tab20")
tissue_legend = [mpatches.Patch(color=(0.17, 0.55, 0.24), label="epithelial"),
                 mpatches.Patch(color=(0.98, 0.80, 0.86), label="stromal"),
                 mpatches.Patch(color=(0.84, 0.15, 0.16), label="cancer")]
for k, s in enumerate(snaps):
    axes[0, k].imshow(s["img"], interpolation="nearest"); axes[0, k].set_xticks([]); axes[0, k].set_yticks([])
    axes[0, k].set_title(f"{s['ncan']:,} cancer · {s['ngl']}/{tg.n_glands} glands\nstroma {s['stroma']:.0f}%",
                         fontsize=9)
    grid = s["gland_grid"]
    rgb = np.ones((*grid.shape, 3))
    for r in range(grid.shape[0]):
        for c in range(grid.shape[1]):
            v = grid[r, c]
            if np.isnan(v):
                continue
            rgb[r, c] = (0.05, 0.05, 0.05) if v < 0 else gland_cmap(int(v) % 20)[:3]
    axes[1, k].imshow(rgb, interpolation="nearest"); axes[1, k].set_xticks([]); axes[1, k].set_yticks([])
axes[0, 0].legend(handles=tissue_legend, fontsize=7, loc="upper right", framealpha=0.85)
axes[0, 0].set_ylabel("tissue (cells)\ngreen wall · red cancer · pink stroma", fontsize=8)
axes[1, 0].set_ylabel("cancer by gland\ncolour = focus · black = stroma / IDC", fontsize=8)
fig.suptitle("Multi-focal DCIS → IDC growth on the ductal field (one founder)", fontsize=13)
print("snapshots (step, cancer, glands, stroma%):",
      [(s["step"], s["ncan"], s["ngl"], round(s["stroma"])) for s in snaps])

import matplotlib.patches as mpatches MILESTONES = [300, 3000, 20000, 50000, 80000] # cancer-cell counts to snapshot (seeding -> confluent IDC), cm-scale tg = B.new_tumor(max_cells=8000) # ungrown, cm-safe realistic ductal field; stepped below to build the movie snaps, mi = [], 0 while mi < len(MILESTONES): ncan = B.n_cancer(tg) if ncan >= MILESTONES[mi]: tg.make_cell_data() img, _ = B.expanded_tissue_rgb(tg, section_frac=0.4) # fast tissue section (no per-clone Muller colormap) snaps.append(dict(step=tg.step, ncan=ncan, ngl=len(B.glands_colonised(tg)), stroma=B.stroma_cancer_pct(tg), img=img, gland_grid=B.cancer_gland_grid(tg))) mi += 1 continue tg.grow(n_steps=(2 if ncan < 4000 else 5), seed=B.BASE_SEED) # fine early so the DCIS phase is caught ncol = len(snaps) fig, axes = plt.subplots(2, ncol, figsize=(3.0 * ncol, 6.6), constrained_layout=True) gland_cmap = plt.get_cmap("tab20") tissue_legend = [mpatches.Patch(color=(0.17, 0.55, 0.24), label="epithelial"), mpatches.Patch(color=(0.98, 0.80, 0.86), label="stromal"), mpatches.Patch(color=(0.84, 0.15, 0.16), label="cancer")] for k, s in enumerate(snaps): axes[0, k].imshow(s["img"], interpolation="nearest"); axes[0, k].set_xticks([]); axes[0, k].set_yticks([]) axes[0, k].set_title(f"{s['ncan']:,} cancer · {s['ngl']}/{tg.n_glands} glands\nstroma {s['stroma']:.0f}%", fontsize=9) grid = s["gland_grid"] rgb = np.ones((*grid.shape, 3)) for r in range(grid.shape[0]): for c in range(grid.shape[1]): v = grid[r, c] if np.isnan(v): continue rgb[r, c] = (0.05, 0.05, 0.05) if v < 0 else gland_cmap(int(v) % 20)[:3] axes[1, k].imshow(rgb, interpolation="nearest"); axes[1, k].set_xticks([]); axes[1, k].set_yticks([]) axes[0, 0].legend(handles=tissue_legend, fontsize=7, loc="upper right", framealpha=0.85) axes[0, 0].set_ylabel("tissue (cells)\ngreen wall · red cancer · pink stroma", fontsize=8) axes[1, 0].set_ylabel("cancer by gland\ncolour = focus · black = stroma / IDC", fontsize=8) fig.suptitle("Multi-focal DCIS → IDC growth on the ductal field (one founder)", fontsize=13) print("snapshots (step, cancer, glands, stroma%):", [(s["step"], s["ncan"], s["ngl"], round(s["stroma"])) for s in snaps])

snapshots (step, cancer, glands, stroma%): [(1248, 331, 7, 0), (1298, 3019, 8, 1), (1386, 20134, 8, 13), (1446, 55779, 8, 65), (1461, 83660, 8, 76)]

3. Clonal dynamics¶

The tumour is infinite-sites: most divisions spawn a new (passenger-differentiated) genotype, so there are tens of thousands of clones. A raw Muller plot would hang; min_freq merges every clone whose subtree never reaches that fraction of the population into its nearest ancestor (a Noble-style sensitivity threshold). by_drivers=True collapses passenger diversity, colouring by distinct driver-mutation combinations.

In [4]:

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6))
t.plot_muller(ax=axes[0], min_freq=0.03)
axes[0].set_title(f"Muller — clones ≥3% (of {living_clones:,} living clones)")
t.plot_muller(ax=axes[1], min_freq=0.03, by_drivers=True)
axes[1].set_title("Muller — coloured by driver-mutation combination")
plt.tight_layout()

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6)) t.plot_muller(ax=axes[0], min_freq=0.03) axes[0].set_title(f"Muller — clones ≥3% (of {living_clones:,} living clones)") t.plot_muller(ax=axes[1], min_freq=0.03, by_drivers=True) axes[1].set_title("Muller — coloured by driver-mutation combination") plt.tight_layout()

4. The ground-truth matrices¶

The per-cell matrices in cell_data are what every downstream notebook reads. Ordered normal-then- cancer, the structure is already visible: SNVs accumulate in malignant cells, copy number departs from the diploid (CN = 2) normal baseline, and expression carries both the copy-number dosage and the gene-program signal. Two ductal-field additions ride alongside: cell_gland (the focus each cell sits in) and the compartment traits in cell_evo (breach, stromal_survival).

In [5]:

rng = np.random.default_rng(0)
order = np.concatenate([rng.permutation(np.where(~cancer)[0])[:120],
                        rng.permutation(np.where(cancer)[0])[:280]])
split = 120

fig, axes = plt.subplots(1, 3, figsize=(15, 4.6))
axes[0].imshow(cd["cell_snv"].values[order], aspect="auto", cmap="binary", interpolation="none")
axes[0].set_title("cell_snv (mutation present)"); axes[0].set_xlabel("gene")
axes[0].set_ylabel("cell  (normal → cancer)")

im1 = axes[1].imshow(cd["cell_cnv"].values[order], aspect="auto", cmap="RdBu_r", vmin=0, vmax=4,
                     interpolation="none")
axes[1].set_title("cell_cnv (total copy number)"); axes[1].set_xlabel("gene")
plt.colorbar(im1, ax=axes[1], fraction=0.046, label="CN")

E = cd["cell_exp"].values[order]
im2 = axes[2].imshow(np.log1p(E / E.sum(1, keepdims=True) * 1e4), aspect="auto", cmap="viridis",
                     interpolation="none")
axes[2].set_title("cell_exp (log CP10k)"); axes[2].set_xlabel("gene")
plt.colorbar(im2, ax=axes[2], fraction=0.046, label="log1p")
for ax_ in axes:
    ax_.axhline(split, color="k", lw=1.2, ls="--")
plt.tight_layout()
print("dashed line = normal/cancer split;  cell_data keys:", list(cd.keys()))
br = B.cell_trait(t, "breach"); ss = B.cell_trait(t, "stromal_survival")
print(f"compartment traits (cancer): breach mean {br[cancer].mean():.2f} · "
      f"stromal_survival mean {ss[cancer].mean():.2f}  (both ≈0 in normal cells)")

rng = np.random.default_rng(0) order = np.concatenate([rng.permutation(np.where(~cancer)[0])[:120], rng.permutation(np.where(cancer)[0])[:280]]) split = 120 fig, axes = plt.subplots(1, 3, figsize=(15, 4.6)) axes[0].imshow(cd["cell_snv"].values[order], aspect="auto", cmap="binary", interpolation="none") axes[0].set_title("cell_snv (mutation present)"); axes[0].set_xlabel("gene") axes[0].set_ylabel("cell (normal → cancer)") im1 = axes[1].imshow(cd["cell_cnv"].values[order], aspect="auto", cmap="RdBu_r", vmin=0, vmax=4, interpolation="none") axes[1].set_title("cell_cnv (total copy number)"); axes[1].set_xlabel("gene") plt.colorbar(im1, ax=axes[1], fraction=0.046, label="CN") E = cd["cell_exp"].values[order] im2 = axes[2].imshow(np.log1p(E / E.sum(1, keepdims=True) * 1e4), aspect="auto", cmap="viridis", interpolation="none") axes[2].set_title("cell_exp (log CP10k)"); axes[2].set_xlabel("gene") plt.colorbar(im2, ax=axes[2], fraction=0.046, label="log1p") for ax_ in axes: ax_.axhline(split, color="k", lw=1.2, ls="--") plt.tight_layout() print("dashed line = normal/cancer split; cell_data keys:", list(cd.keys())) br = B.cell_trait(t, "breach"); ss = B.cell_trait(t, "stromal_survival") print(f"compartment traits (cancer): breach mean {br[cancer].mean():.2f} · " f"stromal_survival mean {ss[cancer].mean():.2f} (both ≈0 in normal cells)")

dashed line = normal/cancer split;  cell_data keys: ['cell_evo', 'cell_snv', 'cell_cnv', 'cell_exp', 'cell_crd', 'cell_type', 'cell_deme', 'cell_snv_p', 'cell_snv_m', 'cell_rna_vaf', 'cell_wgd', 'cell_gland', 'cell_exp_p', 'cell_exp_m', 'cell_rna_baf', 'cell_program']
compartment traits (cancer): breach mean 0.97 · stromal_survival mean 0.74  (both ≈0 in normal cells)

5. What each notebook uses¶

notebook	reads	recovers (scored vs iscc ground truth)
compartment_selection_confound	cell_gland, cell_evo, cell_program, microenv_truth	the DCIS→IDC transition + the genetic-vs-niche emt confound
Integrating scDNA and scRNA with clonealign	cell_cnv, cell_exp (scDNA + scRNA)	subclones from RNA via DNA-defined CN profiles (clonealign idea)
Tree reconstruction from single-cell DNA with SCITE	cell_snv via single-cell DNA	the mutation tree, from noisy genotype calls
Calling copy number from single-cell DNA with HMMcopy	cell_cnv via single-cell read depth	per-clone copy number, which clonealign then reads
Clone inference from bulk DNA with PyClone-VI	cell_snv pooled into bulk read counts	clonal clusters and their cancer-cell fractions
Deconvolving Visium spots with RCTD	Visium spots + scRNA reference	per-spot cell-type composition across foci by NNLS deconvolution
Integrating a patient cohort with scDEF	5 tumours (shared landscape)	programs shared across patients despite batch effects (scDEF / NMF)
Recovering progression constraints with TreeMHN	a cohort of mutation trees	the planted precedence constraints, from event order (TreeMHN)

Next: open any notebook above — each starts with the same import base_sim and grow_base_tumor(). The mechanics live in 01_pipeline_walkthrough, and the assay_* notebooks; the production-tool benchmarks and the full ductal-field / compartment validations live under validation/ (validate_ductal_field.py, validate_compartment_selection.py).