iscc overview¶

iscc simulates a tumour from its first transformed cell to the molecular data a study would actually collect — and records the truth behind every step, which no real experiment can.

This notebook is the whole arc in one pass: grow a lesion, watch it evolve, look at the ground truth, cut a specimen, and run the assays over it — bulk DNA, single-cell DNA, single-cell RNA and spatial transcriptomics — with the answer sitting beside each result.

Everything below runs inline through the public API on one tumour, so what you read is what you would write. The other tutorials go deeper on each stage: Exploring the tumor on the evolution, the Data generation notebooks on each assay, and the Data analysis examples on running real published tools over the output.

In [1]:

%matplotlib inline
import os
from collections import Counter

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scanpy as sc
import anndata as ad
import squidpy as sq

import iscc
from iscc.constants import normal_names
from iscc.tumor.models import GenotypeTumor
from iscc.sample import Resection
from iscc.data import bulkDNA, scDNA, scRNA, Visium

sc.settings.verbosity = 0
NB_DIR = os.getcwd()
CONFIG = os.path.join(NB_DIR, "example_config.yaml")
def coarse_type(v):                      # collapse a per-cell label to cancer / normal type
    return v if v in normal_names else "cancer"


print("iscc from:", iscc.__file__)

%matplotlib inline import os from collections import Counter import numpy as np import pandas as pd import matplotlib.pyplot as plt import scanpy as sc import anndata as ad import squidpy as sq import iscc from iscc.constants import normal_names from iscc.tumor.models import GenotypeTumor from iscc.sample import Resection from iscc.data import bulkDNA, scDNA, scRNA, Visium sc.settings.verbosity = 0 NB_DIR = os.getcwd() CONFIG = os.path.join(NB_DIR, "example_config.yaml") def coarse_type(v): # collapse a per-cell label to cancer / normal type return v if v in normal_names else "cancer" print("iscc from:", iscc.__file__)

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

Grow a tumor¶

We build a spatial GenotypeTumor from a YAML config and grow it. The config sets the spatial mode, genome size, the CINner-style selection parameters, and per-cell-type rates.

In [2]:

from iscc.tumor.models import GenotypeTumor

# One cm-scale tissue (~8.5 mm, ~10^6 cells): a single-clonal-founder DCIS -> confluent-IDC breast sample.
# The lesion spends a long opening phase shut inside its founding acinus, growing but going nowhere; that
# is where its clonal backbone is made, since whatever takes over the patch before a cell escapes ends up
# in every cell afterwards. Only then does it spread and invade, which is the part that takes the time.
tumor = GenotypeTumor(config=CONFIG, seed=2)
for _ in range(100):                      # grow to the confluent DCIS->IDC of the landing hero (~190k
    tumor.grow(n_steps=20, seed=2)        # cancer). Invasion is slow and breach-gated, so grow until it
    if tumor.get_cancer_size() >= 150_000:  # coalesces rather than for a fixed number of steps — most of
        break                             # the budget goes on the confined phase before anything moves.
print(f"Tumor: {tumor.get_cancer_size()} cancer / {tumor.get_tumor_size()} total cells "
      f"across {len(tumor.demes)} demes ({tumor.grid_size} x {tumor.grid_size} ~ 8.5 mm)")

from iscc.tumor.models import GenotypeTumor # One cm-scale tissue (~8.5 mm, ~10^6 cells): a single-clonal-founder DCIS -> confluent-IDC breast sample. # The lesion spends a long opening phase shut inside its founding acinus, growing but going nowhere; that # is where its clonal backbone is made, since whatever takes over the patch before a cell escapes ends up # in every cell afterwards. Only then does it spread and invade, which is the part that takes the time. tumor = GenotypeTumor(config=CONFIG, seed=2) for _ in range(100): # grow to the confluent DCIS->IDC of the landing hero (~190k tumor.grow(n_steps=20, seed=2) # cancer). Invasion is slow and breach-gated, so grow until it if tumor.get_cancer_size() >= 150_000: # coalesces rather than for a fixed number of steps — most of break # the budget goes on the confined phase before anything moves. print(f"Tumor: {tumor.get_cancer_size()} cancer / {tumor.get_tumor_size()} total cells " f"across {len(tumor.demes)} demes ({tumor.grid_size} x {tumor.grid_size} ~ 8.5 mm)")

Tumor: 225623 cancer / 491219 total cells across 28900 demes (170 x 170 ~ 8.5 mm)

What the config sets — and the physical scale¶

example_config.yaml defines this tissue. Three things are worth reading off it:

• A deme is not a cell. The engine's spatial unit is a deme (~50 µm) — a small 3-D column of tissue, not one cell. Its carrying capacity K is that column's cell population (K_duct 60 in a dense duct, K_stroma 30 in stroma); that is how a 2-D grid captures 3-D depth. A histology or Visium section samples only a thin slice of each column.
• The sizes are realistic. grid 170 × ~50 µm ≈ 8.5 mm of tissue (~0.5 M cells grown, ~10⁶ at capacity) — deliberately larger than a 6.5 mm Visium capture, so every assay samples a subset. A ductal carcinoma in situ (DCIS) duct is gland_radius 4 → 9 demes ≈ 450 µm; one Visium spot ≈ 55 µm ≈ one deme.
• The biology is one clonal founder, DCIS → invasive ductal carcinoma (IDC). selection_params/cell_params set CINner-style fitness (oncogene/TSG copy-number drivers) and rare, strong invasion: crossing the duct wall (breach) is rare, but the stroma is permissive, so an escaped clone invades and spreads along the connected ducts. n_steps sets how far invasion has coalesced — all the invasive tissue is one clone. Tau-leaping (update_mode: tau) + coarsen_passengers make cost scale with clones × demes, not cells, which is what makes ~10⁶ cells affordable.

Visualise growth and spatial structure¶

Three views of the same tumour, sharing one colouring so they line up: the tissue grid shows where clones sit (space), the Muller shows when they swept (time), and the phylogeny shows who descends from whom (ancestry). Clones are coloured by their driver clone — the distinct set of driver (oncogene/TSG) mutations they carry (by_drivers) — so a clone keeps its colour across all three.

The grid below also has a stage view: instead of driver identity, it colours each region by the stage-dominant trait its clone has reached (none → proliferation → invasion). That's the same scheme the landing animation uses to show the full metastatic cascade — here, with no metastasis or treatment, it just reads out the DCIS → IDC progression.

In [3]:

# plot_tissue draws the deme grid from the per-deme COUNTS, so it stays complete at cm-scale even
# though cell_data is a memory-safe subsample. "state" is categorical histology; "clone" colours each
# cancer deme by its dominant driver clone (matching the Muller + phylogeny); "stage" colours it by the
# stage-dominant trait (none / proliferation / invasion) — the SAME scheme the landing animation uses.
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
tumor.plot_tissue(ax=axes[0], color="state")
axes[0].set_title("Tissue (state): DCIS ducts + focal IDC")
tumor.plot_tissue(ax=axes[1], color="clone")
axes[1].set_title("Tissue (by driver clone)")
tumor.plot_tissue(ax=axes[2], color="stage")
axes[2].set_title("Tissue (by stage-dominant trait)")
plt.tight_layout()

# plot_tissue draws the deme grid from the per-deme COUNTS, so it stays complete at cm-scale even # though cell_data is a memory-safe subsample. "state" is categorical histology; "clone" colours each # cancer deme by its dominant driver clone (matching the Muller + phylogeny); "stage" colours it by the # stage-dominant trait (none / proliferation / invasion) — the SAME scheme the landing animation uses. fig, axes = plt.subplots(1, 3, figsize=(16, 5)) tumor.plot_tissue(ax=axes[0], color="state") axes[0].set_title("Tissue (state): DCIS ducts + focal IDC") tumor.plot_tissue(ax=axes[1], color="clone") axes[1].set_title("Tissue (by driver clone)") tumor.plot_tissue(ax=axes[2], color="stage") axes[2].set_title("Tissue (by stage-dominant trait)") plt.tight_layout()

The Muller (time) and the phylogeny (ancestry) complete the trio, in the same driver-clone colours. by_drivers collapses passenger-only diversity to a handful of driver clades; the Muller shows their frequencies over time and the phylogeny shows how they descend from the founder — the same clades you see tiling the grid.

In [4]:

fig, ax = plt.subplots(figsize=(9, 4))
tumor.plot_muller(ax=ax, by_drivers=True, min_freq=0.02)
ax.set_title("Clonal dynamics (Muller, by driver clone) — a single founder + its subclones")

fig, ax = plt.subplots(figsize=(9, 4)) tumor.plot_muller(ax=ax, by_drivers=True, min_freq=0.02) ax.set_title("Clonal dynamics (Muller, by driver clone) — a single founder + its subclones")

Out[4]:

Text(0.5, 1.0, 'Clonal dynamics (Muller, by driver clone) — a single founder + its subclones')

In [5]:

# Radial CLADE phylogeny, drawn straight to the axes (no image file). Rather than subsampling cells,
# the population is collapsed to driver clones (the same collapse the Muller uses) — one node per clade,
# sized by cell count, in the SAME driver-clone colours as the Muller and the grid.
fig, ax = plt.subplots(figsize=(6, 6))
tumor.plot_phylogeny(ax=ax)
ax.set_title("Clade phylogeny (by driver clone)")

# Radial CLADE phylogeny, drawn straight to the axes (no image file). Rather than subsampling cells, # the population is collapsed to driver clones (the same collapse the Muller uses) — one node per clade, # sized by cell count, in the SAME driver-clone colours as the Muller and the grid. fig, ax = plt.subplots(figsize=(6, 6)) tumor.plot_phylogeny(ax=ax) ax.set_title("Clade phylogeny (by driver clone)")

Out[5]:

Text(0.5, 1.0, 'Clade phylogeny (by driver clone)')

The ground truth¶

The per-cell ground truth lives in tumor.cell_data (a dict of DataFrames, also written to cell_data/ by tumor.write(path)). These matrices are what later stages sample and sequence.

In [6]:

cd = tumor.cell_data
print("cell_data matrices:", list(cd.keys()))
cd["cell_snv"].iloc[:5, :6]

cd = tumor.cell_data print("cell_data matrices:", list(cd.keys())) cd["cell_snv"].iloc[:5, :6]

cell_data matrices: ['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']

Out[6]:

	G_0_0	G_0_1	G_0_2	G_0_3	G_0_4	G_0_5
C0	0.0	0.0	0.0	0.0	0.0	0.0
C1	0.0	0.0	0.0	0.0	0.0	0.0
C2	0.0	0.0	0.0	0.0	0.0	0.0
C3	0.0	0.0	0.0	0.0	0.0	0.0
C4	0.0	0.0	0.0	0.0	0.0	0.0

In [7]:

from iscc.tumor.viz import CNV_CMAP, cnv_norm
fig, axes = plt.subplots(1, 3, figsize=(13, 3))
for ax, key, title in zip(axes, ["cell_snv", "cell_cnv", "cell_exp"], ["SNV", "Copy number", "Expression"]):
    if key == "cell_cnv":  # copy number: blue-white-red, white at diploid (2)
        im = ax.imshow(cd[key].values, aspect="auto", cmap=CNV_CMAP, norm=cnv_norm(cd[key].values))
    else:
        im = ax.imshow(cd[key].values, aspect="auto", cmap="viridis")
    ax.set_title(title); ax.set_xlabel("gene"); ax.set_ylabel("cell")
    plt.colorbar(im, ax=ax, fraction=0.046)
plt.tight_layout()

from iscc.tumor.viz import CNV_CMAP, cnv_norm fig, axes = plt.subplots(1, 3, figsize=(13, 3)) for ax, key, title in zip(axes, ["cell_snv", "cell_cnv", "cell_exp"], ["SNV", "Copy number", "Expression"]): if key == "cell_cnv": # copy number: blue-white-red, white at diploid (2) im = ax.imshow(cd[key].values, aspect="auto", cmap=CNV_CMAP, norm=cnv_norm(cd[key].values)) else: im = ax.imshow(cd[key].values, aspect="auto", cmap="viridis") ax.set_title(title); ax.set_xlabel("gene"); ax.set_ylabel("cell") plt.colorbar(im, ax=ax, fraction=0.046) plt.tight_layout()

Cutting the specimen¶

We sample the specimen with two orthogonal cuts (iscc.sample.Resection). An in-plane cut (bisect) splits the 2-D tissue: one part is dissociated for the sequencing assays, the remainder is kept for the Visium section. (The depth cut for Visium comes later.) Resection materialises only each piece, so this stays memory-safe at cm-scale.

In [8]:

spec = Resection(tumor)
CUT_FRAC = 0.42
cut, remainder = spec.bisect(frac=CUT_FRAC, axis="x")   # in-plane cut: dissociation | imaging
split = int(CUT_FRAC * tumor.grid_size)

fig, ax = plt.subplots(figsize=(5.6, 5.6))
tumor.plot_tissue(ax=ax, color="state")
ax.axvline(split - 0.5, color="k", lw=2.5, ls="--")
box = dict(boxstyle="round", fc="white", ec="0.5", alpha=0.85)
ax.text(split * 0.5, 6, "dissociation\ncut", ha="center", va="top", bbox=box)
ax.text((split + tumor.grid_size) * 0.5, 6, "Visium\nremainder", ha="center", va="top", bbox=box)
ax.set_title("resected specimen, cut in two")
print(f"in-plane cut: {len(cut)} demes dissociated | {len(remainder)} demes kept for imaging")

spec = Resection(tumor) CUT_FRAC = 0.42 cut, remainder = spec.bisect(frac=CUT_FRAC, axis="x") # in-plane cut: dissociation | imaging split = int(CUT_FRAC * tumor.grid_size) fig, ax = plt.subplots(figsize=(5.6, 5.6)) tumor.plot_tissue(ax=ax, color="state") ax.axvline(split - 0.5, color="k", lw=2.5, ls="--") box = dict(boxstyle="round", fc="white", ec="0.5", alpha=0.85) ax.text(split * 0.5, 6, "dissociation\ncut", ha="center", va="top", bbox=box) ax.text((split + tumor.grid_size) * 0.5, 6, "Visium\nremainder", ha="center", va="top", bbox=box) ax.set_title("resected specimen, cut in two") print(f"in-plane cut: {len(cut)} demes dissociated | {len(remainder)} demes kept for imaging")

in-plane cut: 12070 demes dissociated | 16830 demes kept for imaging

The dissociated sample and its ground truth¶

Dissociate the cut into a per-cell table — the sample the sequencing assays sequence. Because iscc simulated every cell, it knows the ground truth: each cell's true type, clone, copy number and single-nucleotide variants (SNVs). Here is the true SNV matrix of the dissociated cells (the assays below are noisy views of it).

In [9]:

cd = spec.dissociate(cut, max_cells=8000)          # dissociate the cut -> the sequencing sample
true_type = cd["cell_type"]["cell_id"].astype(str)  # ground-truth cell type per cell
drivers = tumor.get_gene_data()["driver_types"].iloc[:, 0]   # +1 oncogene, -1 TSG, 0 passenger
n_can = int(cd["cell_type"]["cell_id"].map(lambda g: tumor.genotypes[g].type == "cancer").sum())
print(f"dissociated {len(cd['cell_snv'])} cells ({n_can} cancer, {len(cd['cell_snv']) - n_can} "
      f"microenvironment) | {len(drivers)} genes, {int((drivers == 1).sum())} oncogenes, "
      f"{int((drivers == -1).sum())} tumour suppressors")

fig, ax = plt.subplots(figsize=(8, 4))
im = ax.imshow(cd["cell_snv"].values, aspect="auto", cmap="Greys", interpolation="none")
ax.set_xlabel("gene"); ax.set_ylabel("cell"); ax.set_title("True SNV matrix (dissociated cells)")
plt.colorbar(im, ax=ax, fraction=0.046)

cd = spec.dissociate(cut, max_cells=8000) # dissociate the cut -> the sequencing sample true_type = cd["cell_type"]["cell_id"].astype(str) # ground-truth cell type per cell drivers = tumor.get_gene_data()["driver_types"].iloc[:, 0] # +1 oncogene, -1 TSG, 0 passenger n_can = int(cd["cell_type"]["cell_id"].map(lambda g: tumor.genotypes[g].type == "cancer").sum()) print(f"dissociated {len(cd['cell_snv'])} cells ({n_can} cancer, {len(cd['cell_snv']) - n_can} " f"microenvironment) | {len(drivers)} genes, {int((drivers == 1).sum())} oncogenes, " f"{int((drivers == -1).sum())} tumour suppressors") fig, ax = plt.subplots(figsize=(8, 4)) im = ax.imshow(cd["cell_snv"].values, aspect="auto", cmap="Greys", interpolation="none") ax.set_xlabel("gene"); ax.set_ylabel("cell"); ax.set_title("True SNV matrix (dissociated cells)") plt.colorbar(im, ax=ax, fraction=0.046)

/Users/pedroferreira/projects/iscc/repo/src/iscc/sample/resection/resection.py:153: UserWarning: max_cells=8000 capped this section: 251,709 cells were selected by region/depth_frac but only ~8,000 are materialised (3.2%). The cap, not depth_frac, is setting the cell density. Pass a larger max_cells (~302,050) to let the physical section govern.
  return self.tumor.make_cell_data(**kwargs)

dissociated 8024 cells (4451 cancer, 3573 microenvironment) | 6000 genes, 130 oncogenes, 106 tumour suppressors

Out[9]:

<matplotlib.colorbar.Colorbar at 0x39257df00>

Bulk DNA-seq¶

Bulk DNA pools the dissociated cells into one library: per-gene coverage and alt-read counts, from which we read a variant-allele-frequency (VAF) spectrum.

In [10]:

bd = bulkDNA(n_reads=50000, data_mode="counts", fpr=0.02, fnr=0.02)
bd.run(cd)
vaf = bd.observed_data["vaf"][bd.observed_data["coverage"] > 0]   # drop zero-coverage loci
print("top variant genes:", list(vaf.sort_values(ascending=False).head(5).index))

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(vaf, bins=30); ax.set_xlabel("VAF"); ax.set_ylabel("genes")
ax.set_title("Bulk DNA VAF spectrum")

bd = bulkDNA(n_reads=50000, data_mode="counts", fpr=0.02, fnr=0.02) bd.run(cd) vaf = bd.observed_data["vaf"][bd.observed_data["coverage"] > 0] # drop zero-coverage loci print("top variant genes:", list(vaf.sort_values(ascending=False).head(5).index)) fig, ax = plt.subplots(figsize=(6, 3)) ax.hist(vaf, bins=30); ax.set_xlabel("VAF"); ax.set_ylabel("genes") ax.set_title("Bulk DNA VAF spectrum")

top variant genes: ['G_5_483', 'G_2_120', 'G_6_300', 'G_2_126', 'G_9_211']

Out[10]:

Text(0.5, 1.0, 'Bulk DNA VAF spectrum')

Single-cell DNA-seq¶

scDNA resolves SNVs per cell. We compute a per-cell VAF matrix (alt / coverage) — the substrate for clone reconstruction.

In [11]:

sd = scDNA(n_cells=200, breadth="wgs")
sd.run(cd)
sc_vaf = sd.vaf.fillna(0)
keep = sc_vaf.columns[sc_vaf.sum(0) > 0]   # genes with any signal

fig, ax = plt.subplots(figsize=(8, 4))
im = ax.imshow(sc_vaf[keep].values, aspect="auto", cmap="magma")
ax.set_xlabel("variant gene"); ax.set_ylabel("cell"); ax.set_title("scDNA per-cell VAF")
plt.colorbar(im, ax=ax, fraction=0.046)

sd = scDNA(n_cells=200, breadth="wgs") sd.run(cd) sc_vaf = sd.vaf.fillna(0) keep = sc_vaf.columns[sc_vaf.sum(0) > 0] # genes with any signal fig, ax = plt.subplots(figsize=(8, 4)) im = ax.imshow(sc_vaf[keep].values, aspect="auto", cmap="magma") ax.set_xlabel("variant gene"); ax.set_ylabel("cell"); ax.set_title("scDNA per-cell VAF") plt.colorbar(im, ax=ax, fraction=0.046)

Out[11]:

<matplotlib.colorbar.Colorbar at 0x392489120>

Single-cell RNA-seq¶

scRNA gives unique-molecular-identifier (UMI) counts per cell. We turn the assay into an AnnData, annotate each cell with its true type, then run a standard scanpy workflow (quality control, normalization, PCA, neighbours, Leiden, UMAP) — and colour the UMAP by the true cell type to see how well the unsupervised clusters recover it.

In [12]:

sr = scRNA(n_reads=4000, n_cells=400)
sr.run(cd)
adata = sr.to_anndata()
adata.obs["cell_type"] = pd.Categorical(true_type.reindex(adata.obs_names).map(coarse_type))
sc.pp.calculate_qc_metrics(adata, inplace=True, percent_top=None)
sc.pp.filter_genes(adata, min_cells=3)
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.pca(adata, n_comps=20)
sc.pp.neighbors(adata, n_neighbors=15)
sc.tl.leiden(adata, resolution=1.0, flavor="igraph", n_iterations=2, directed=False)
sc.tl.umap(adata)
print("cell-type annotations:", dict(adata.obs["cell_type"].value_counts()))

sc.pl.umap(adata, color=["cell_type", "leiden", "total_counts"], wspace=0.4)

sr = scRNA(n_reads=4000, n_cells=400) sr.run(cd) adata = sr.to_anndata() adata.obs["cell_type"] = pd.Categorical(true_type.reindex(adata.obs_names).map(coarse_type)) sc.pp.calculate_qc_metrics(adata, inplace=True, percent_top=None) sc.pp.filter_genes(adata, min_cells=3) adata.layers["counts"] = adata.X.copy() sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.pca(adata, n_comps=20) sc.pp.neighbors(adata, n_neighbors=15) sc.tl.leiden(adata, resolution=1.0, flavor="igraph", n_iterations=2, directed=False) sc.tl.umap(adata) print("cell-type annotations:", dict(adata.obs["cell_type"].value_counts())) sc.pl.umap(adata, color=["cell_type", "leiden", "total_counts"], wspace=0.4)

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

cell-type annotations: {'cancer': 215, 'stromal': 178, 'epithelial': 7}

In [13]:

sc.tl.rank_genes_groups(adata, groupby="cell_type", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=8, sharey=False)

sc.tl.rank_genes_groups(adata, groupby="cell_type", method="wilcoxon") sc.pl.rank_genes_groups(adata, n_genes=8, sharey=False)

Visium spatial transcriptomics¶

Visium is a spatial assay. From the imaging remainder we take a depth cut (spec.slice(remainder, depth_frac=0.2)) — a thin section keeping the whole 2-D field but only a fifth of each deme's 3-D column, which is about what a ~10 µm microtome section takes out of a ~50 µm-wide deme. First the H&E of the whole remainder (the tissue we image):

In [14]:

# The section is the WHOLE imaging remainder, not a patch of it: the v1 capture area is 156 x 111
# deme-widths and the remainder is 99 x 170, so essentially all of the tissue fits on one slide.
# max_cells has to be raised above what depth_frac implies, or the tumour's own 8,000-cell
# materialisation cap binds instead and every spot comes out near-empty.
section = spec.slice(remainder, depth_frac=0.2, max_cells=45000)

he, px = tumor.he_image(px=6)
fig, ax = plt.subplots(figsize=(4.4, 6.6))
ax.imshow(he[:, int(CUT_FRAC * tumor.grid_size * px):])      # crop to the remainder (right of the cut)
ax.set_xticks([]); ax.set_yticks([])
ax.set_title("H&E, imaging remainder\n(dense nuclei = invasive; pale = stroma)")

per_deme = Counter(section["cell_deme"]["deme_id"])
print(f"section: {len(section['cell_snv']):,} cells over {len(per_deme):,} of {len(remainder):,} demes"
      f" | median {int(np.median(list(per_deme.values())))} cells/deme, max {max(per_deme.values())}")

# The section is the WHOLE imaging remainder, not a patch of it: the v1 capture area is 156 x 111 # deme-widths and the remainder is 99 x 170, so essentially all of the tissue fits on one slide. # max_cells has to be raised above what depth_frac implies, or the tumour's own 8,000-cell # materialisation cap binds instead and every spot comes out near-empty. section = spec.slice(remainder, depth_frac=0.2, max_cells=45000) he, px = tumor.he_image(px=6) fig, ax = plt.subplots(figsize=(4.4, 6.6)) ax.imshow(he[:, int(CUT_FRAC * tumor.grid_size * px):]) # crop to the remainder (right of the cut) ax.set_xticks([]); ax.set_yticks([]) ax.set_title("H&E, imaging remainder\n(dense nuclei = invasive; pale = stroma)") per_deme = Counter(section["cell_deme"]["deme_id"]) print(f"section: {len(section['cell_snv']):,} cells over {len(per_deme):,} of {len(remainder):,} demes" f" | median {int(np.median(list(per_deme.values())))} cells/deme, max {max(per_deme.values())}")

/Users/pedroferreira/projects/iscc/repo/src/iscc/sample/resection/resection.py:153: UserWarning: max_cells=45000 capped this section: 47,985 cells were selected by region/depth_frac but only ~45,000 are materialised (93.8%). The cap, not depth_frac, is setting the cell density. Pass a larger max_cells (~57,582) to let the physical section govern.
  return self.tumor.make_cell_data(**kwargs)

section: 44,957 cells over 14,883 of 16,830 demes | median 2 cells/deme, max 28

The remainder is tall and narrow (99 × 170 deme-widths) and the v1 capture area is wide and short (78 × 64 spots = 156 × 111 at this pitch), so the section goes on the slide rotated 90° — that lays the tissue's long axis along the slide's long axis and puts nearly all of it under spots. vz.place_grid(section) lays the fixed grid down without assaying it, so we can check the placement first. Mind the rotation: the top of the H&E above is the left of every spot plot below.

In [15]:

vz = Visium(n_reads=20000, spot_pitch=2.0, spot_radius=0.55, section_frac=1.0,
            rotation=90, seed=0)              # placement defaults to the section's own centroid
vz.place_grid(section)
grid = vz.to_anndata()                        # spot grid + section H&E, X all zero (placement only)
lib = list(grid.uns["spatial"])[0]
sq.pl.spatial_scatter(grid, color="in_tissue", library_id=lib, img=True, shape="circle", size=1.0,
                      alpha=0.7, title="v1 spot grid on the section (placed, rotated 90 deg)")
print(f"{int((grid.obs['in_tissue'] == 1).sum())}/{grid.n_obs} spots land on tissue")

vz = Visium(n_reads=20000, spot_pitch=2.0, spot_radius=0.55, section_frac=1.0, rotation=90, seed=0) # placement defaults to the section's own centroid vz.place_grid(section) grid = vz.to_anndata() # spot grid + section H&E, X all zero (placement only) lib = list(grid.uns["spatial"])[0] sq.pl.spatial_scatter(grid, color="in_tissue", library_id=lib, img=True, shape="circle", size=1.0, alpha=0.7, title="v1 spot grid on the section (placed, rotated 90 deg)") print(f"{int((grid.obs['in_tissue'] == 1).sum())}/{grid.n_obs} spots land on tissue")

3872/4992 spots land on tissue

Now run the assay over the placed grid (vz.run() reuses the placement) and analyse the spots with squidpy and scanpy exactly as we would a real Visium slide.

In [16]:

vz.run()                                          # reuse the placed grid, now assay it
vis = vz.to_anndata()
vis = vis[vis.obs["in_tissue"] == 1].copy()          # keep the in-tissue spots
sc.pp.calculate_qc_metrics(vis, inplace=True, percent_top=None)
sc.pp.normalize_total(vis); sc.pp.log1p(vis)
sc.pp.pca(vis, n_comps=20); sc.pp.neighbors(vis)
sc.tl.leiden(vis, resolution=0.8, flavor="igraph", n_iterations=2, directed=False)
print(f"Visium: {vis.n_obs} in-tissue spots | median {int(np.median(vis.obs['n_cells']))} cells/spot "
      f"| {vis.obs['leiden'].nunique()} spatial niches")

# squidpy draws the spots straight onto the assay's own H&E image (img=True).
umi_vmax = float(np.percentile(vis.obs["total_counts"], 95))
fig, axes = plt.subplots(1, 2, figsize=(15, 6.2))
sq.pl.spatial_scatter(vis, color="total_counts", library_id=lib, ax=axes[0], fig=fig, img=True,
                      shape="circle", size=1.0, alpha=0.8, vmax=umi_vmax, cmap="magma",
                      title="UMI / spot on H&E")
sq.pl.spatial_scatter(vis, color="leiden", library_id=lib, ax=axes[1], fig=fig, img=True,
                      shape="circle", size=1.0, alpha=0.9, title="unsupervised niches (leiden)")
plt.tight_layout()

vz.run() # reuse the placed grid, now assay it vis = vz.to_anndata() vis = vis[vis.obs["in_tissue"] == 1].copy() # keep the in-tissue spots sc.pp.calculate_qc_metrics(vis, inplace=True, percent_top=None) sc.pp.normalize_total(vis); sc.pp.log1p(vis) sc.pp.pca(vis, n_comps=20); sc.pp.neighbors(vis) sc.tl.leiden(vis, resolution=0.8, flavor="igraph", n_iterations=2, directed=False) print(f"Visium: {vis.n_obs} in-tissue spots | median {int(np.median(vis.obs['n_cells']))} cells/spot " f"| {vis.obs['leiden'].nunique()} spatial niches") # squidpy draws the spots straight onto the assay's own H&E image (img=True). umi_vmax = float(np.percentile(vis.obs["total_counts"], 95)) fig, axes = plt.subplots(1, 2, figsize=(15, 6.2)) sq.pl.spatial_scatter(vis, color="total_counts", library_id=lib, ax=axes[0], fig=fig, img=True, shape="circle", size=1.0, alpha=0.8, vmax=umi_vmax, cmap="magma", title="UMI / spot on H&E") sq.pl.spatial_scatter(vis, color="leiden", library_id=lib, ax=axes[1], fig=fig, img=True, shape="circle", size=1.0, alpha=0.9, title="unsupervised niches (leiden)") plt.tight_layout()

Visium: 3872 in-tissue spots | median 2 cells/spot | 3 spatial niches

The same pipeline from the command line¶

The same grow → sample → assay steps as commands (the CLI takes a punch biopsy rather than the notebook's Resection cut, so the exact sampled cells differ):

isccsim    --sim-config example_config.yaml -r 2 -s 320 -o example_out/tumor
isccsample example_out/tumor --method biopsy --biopsy-type punch -o example_out/sample
isccdata   example_out/sample -a scrna --assay-config ../src/iscc/data/assayconfigs/scrna.yaml -o example_out/data_scrna

Next: Exploring the tumor.

Recap¶

From one resected tumour we sampled two ways: an in-plane cut, dissociated for bulk DNA / scDNA / scRNA; and a thin depth-slice of the remainder, imaged with Visium. Each is a noisy, partial measurement of the same ground truth, analysed with the same off-the-shelf tools (scanpy, squidpy) you would use on real data — yet iscc knows the ground truth (each cell's true type, clone, copy number and mutations), so these datasets are ideal benchmarks for tumour-evolution methods.

The assays model realistic technical structure — negative-binomial + dropout scRNA, allele-specific single-cell DNA with allelic dropout, a read/FASTQ mode, and spatial spot mixing with autocorrelation — and their parameters can be fit to real reference data. Further analysis-demo notebooks (mutual-hazard networks, MHN, from DNA; niche identification, batch effects, combining modalities, real-data comparison) build on these datasets.