slug = "human-colon"
release = "2026-10-28"
use_gpu = True
out_dir = "bundles"
smoke = 0
primers_ref = "feat-primers"
cohorts_path = None
catalog_path = None
representatives = 1
dev_note = ""

# Parameters
slug = "human-lung"
release = "2026-11-01"
out_dir = "bundles"
primers_ref = "4d7af340849c00077f9d21c5b54af9e689a4a070"
dev_note = ""

pip install "sddb-primers @ git+https://github.com/theislab/spatialdata-db-client@4d7af340849c00077f9d21c5b54af9e689a4a070#subdirectory=primers"

import inspect
import json
import platform
import time

import anndata as ad
import numpy as np
import pandas as pd
import scanpy as sc
import spatialdata_plot  # noqa: F401  (registers the .pl accessor)
import zarr
from primers.cohort import hd_table, load_cohort, smoke_spec
from primers.embed import CAVEAT, MAX_FLOOR_LOSS, MIN_SHARED_GENES, N_HVG, N_NEIGHBORS, N_PCS, SEED, embed_views
from primers.export import export_bundle, mixing_by_view, with_qc_status
from primers.markers import load_markers, panel_id
from primers.qc import MIN_GENES_FOR_QC, apply_qc, load_thresholds, qc_floors
from spatialdata.transformations import get_transformation

import sddb
from sddb.remote import storage_options

start = time.time()
out_dir = Path(out_dir)
if use_gpu:
    import cupy as cp
    import rapids_singlecell as rsc
    import rmm
    from rmm.allocators.cupy import rmm_cupy_allocator

    rmm.reinitialize(managed_memory=True)
    cp.cuda.set_allocator(rmm_cupy_allocator)


def show_call(func, what):
    """Show what a primers function does and link its source at ``primers_ref``."""
    _, line = inspect.getsourcelines(func)
    module = func.__module__.split(".")[-1]
    url = f"{REPO}/blob/{primers_ref}/primers/primers/{module}.py#L{line}"
    display(Markdown(f"`{func.__module__}.{func.__name__}`: {what} [Source at `{primers_ref}`]({url})."))


def show_table(df, digits=3):
    """Display ``df`` rounded, with missing values written as n/a, and long text columns in full."""
    shown = df.round(digits).astype(object)
    with pd.option_context("display.max_colwidth", None):
        display(shown.where(df.notna(), "n/a"))


catalog = sddb.Catalog.from_file(catalog_path) if catalog_path else sddb.Catalog(version=release)

if smoke:
    show_call(
        smoke_spec,
        "Keeps only the smallest kept member of each of a few papers, taken in turn from each family, for a "
        "quick end-to-end check.",
    )
    n_spec_excluded = len(spec.get("exclude_uids") or [])
    spec = smoke_spec(spec, catalog, smoke)
    print(
        f"Smoke run on {smoke} papers: {len(spec['exclude_uids']) - n_spec_excluded} more uids are excluded "
        f"({len(spec['exclude_uids'])} in all)."
    )
for key in ("slug", "query", "include_uids", "exclude_uids", "families", "marker_set", "rationale"):
    value = spec.get(key)
    if smoke and key == "exclude_uids":
        value = f"{len(value)} uids"
    print(f"{key}: {value}")

slug: human-lung
query: {'organism': 'human', 'tissue': 'lung', 'n_obs__lte': 500000}
include_uids: None
exclude_uids: ['ADgyHXEKFq7YjJTb0001']
families: ['spot', 'imaging']
marker_set: lung
rationale: Normal, fibrotic, sarcoid and cancer lung on Visium, Visium HD and Xenium; many independent Visium studies.

show_call(
    load_cohort,
    "Runs the query on the catalog, adds the included uids, removes the excluded ones, and loads each kept "
    "member's table keyed by gene symbol, joined within its technology family.",
)
adatas, cohort = load_cohort(spec, catalog=catalog)

columns = ["uid", "title", "technology", "family", "paper", "paper_source", "license", "n_obs", "status"]
show_table(cohort[columns].rename(columns={"license": "licence"}))
show_table(cohort["status"].value_counts().rename_axis("status").reset_index(name="members"))
print(f"{int((cohort['status'] == 'kept').sum())} kept, {int((cohort['status'] != 'kept').sum())} dropped.")

233 kept, 6 dropped.

show_call(load_thresholds, "Reads the QC threshold groups from thresholds.yaml.")
thresholds = load_thresholds()
show_table(
    pd.DataFrame(
        [
            {
                "group": name,
                "technologies": ", ".join(g["technologies"]),
                "min_counts": g["min_counts"],
                "min_genes": g["min_genes"],
                "max_pct_mt": g.get("max_pct_mt"),
                "flag: min median counts": g["flag"]["min_median_counts"],
                "flag: min share kept": g["flag"]["min_share_kept"],
                "reason": g["reason"],
            }
            for name, g in thresholds.items()
        ]
    ).set_index("group")
)
print(f"MIN_GENES_FOR_QC = {MIN_GENES_FOR_QC}: QC, and every view, needs at least this many genes.")

MIN_GENES_FOR_QC = 50: QC, and every view, needs at least this many genes.

show_call(
    apply_qc,
    "Leaves out control probes, computes counts, genes and the mitochondrial share per observation, keeps "
    "the observations that pass their technology's filters, and describes each dataset.",
)
qcs = []
for fam in list(adatas):
    adatas[fam], q = apply_qc(adatas[fam], thresholds)
    qcs.append(q.assign(family=fam))
qc = pd.concat(qcs, ignore_index=True)
show_table(qc)

MAX_VIOLIN_DATASETS = 40
for fam, adata in adatas.items():
    obs = adata.obs[["dataset_uid", "total_counts", "n_genes"]].copy()
    sizes = obs["dataset_uid"].value_counts()
    if len(sizes) > MAX_VIOLIN_DATASETS:
        print(f"{fam}: {len(sizes)} datasets, so the {MAX_VIOLIN_DATASETS} largest are shown.")
        obs = obs[obs["dataset_uid"].isin(sizes.index[:MAX_VIOLIN_DATASETS])]
    obs["dataset_uid"] = obs["dataset_uid"].astype(str).astype("category")
    display(Markdown(f"Counts and genes per observation after QC, family `{fam}`."))
    sc.pl.violin(ad.AnnData(obs=obs), ["total_counts", "n_genes"], groupby="dataset_uid", rotation=90, stripplot=False)

spot: 154 datasets, so the 40 largest are shown.

imaging_cell: 79 datasets, so the 40 largest are shown.

show_call(load_markers, "Reads the marker sets from markers.yaml, each with its citation and fallbacks.")
marker_set = load_markers()[spec["marker_set"]]
print(f"Marker set {spec['marker_set']!r}, cited from {', '.join(marker_set['citation'])}.")
show_table(
    pd.DataFrame({"markers": {c: ", ".join(g) for c, g in marker_set["compartments"].items()}}).rename_axis(
        "compartment"
    )
)
kept_after_qc = set(qc.loc[qc["status"] == "kept", "dataset_uid"].astype(str))
technology_of = cohort.set_index("uid")["technology"].astype(str)
panel_of = {
    uid: panel_id(technology_of[uid], genes)
    for adata in adatas.values()
    for uid, genes in adata.uns.get("member_genes", {}).items()
    if uid in kept_after_qc
}
used = [
    {
        "panel": pid,
        "compartment": comp,
        "genes": "not measured" if fb.get("not_measured") else ", ".join(fb.get("genes") or []),
        "note": fb.get("note", ""),
    }
    for pid, entries in (marker_set.get("fallback") or {}).items()
    if pid in set(panel_of.values())
    for comp, fb in entries.items()
]
if used:
    show_table(pd.DataFrame(used))
else:
    print("No fallback of this marker set applies to the panels of this cohort.")

Marker set 'lung', cited from 10.1038/s41586-020-2922-4.

show_call(with_qc_status, "Marks the members that have no cells left after QC as dropped.")
kept_members = with_qc_status(cohort, qc).query("status == 'kept'")
ENSEMBL = r"ENS[A-Z]*G\d{6}"
SYMBOL_COLUMNS = ("symbol", "gene_symbols", "gene_symbol", "gene_name", "feature_name")


def colour_gene(table, panel):
    """(var name, symbol) of the first gene of the first compartment the table measures.

    A compartment's genes are its markers, then the fallback genes for this section's own panel.
    """
    names = pd.Index(table.var_names.astype(str))
    symbols = names
    if names.str.match(ENSEMBL).mean() >= 0.5:
        column = next((c for c in SYMBOL_COLUMNS if c in table.var), None)
        if column is None:
            return None, None
        symbols = pd.Index(table.var[column].astype(str))
    by_symbol = {s.upper(): (n, s) for n, s in zip(names, symbols, strict=True)}
    fallback = (marker_set.get("fallback") or {}).get(panel) or {}
    for comp, genes in marker_set["compartments"].items():
        for gene in [*genes, *((fallback.get(comp) or {}).get("genes") or [])]:
            if gene.upper() in by_symbol:
                return by_symbol[gene.upper()]
    return None, None


def table_region(zarr_url, table_name):
    """The shapes element a table annotates, read from the table's zarr attributes without loading it."""
    group = zarr.open_group(zarr_url, mode="r", storage_options=storage_options(zarr_url) or None)
    attrs = dict(group[f"tables/{table_name}"].attrs)
    if "region" not in attrs and f"tables/{table_name}/uns/spatialdata_attrs" in group:
        attrs = ad.io.read_elem(group[f"tables/{table_name}/uns/spatialdata_attrs"])
    region = attrs.get("region")
    return region if isinstance(region, str) or region is None else (list(region) or [None])[0]


def representative_uids(members):
    """The ``representatives`` kept sections closest to the median number of observations."""
    n_obs = members.set_index("uid")["n_obs"].astype(float)
    return list((n_obs - n_obs.median()).abs().sort_values(kind="stable").index[:representatives])

for tech, members in kept_members.groupby("technology"):
    for uid in representative_uids(members):
        ds = catalog.dataset(uid)
        stored = ds.elements()
        table_name = hd_table(stored) if tech == "VisiumHD" else (getattr(ds, "default_table", None) or "table")
        region = table_region(ds.zarr_url, table_name)
        if region is None or f"shapes/{region}" not in stored:
            display(Markdown(f"{tech}, section `{uid}`: the table annotates no shapes element, so no plot."))
            continue
        # a small tissue image when there is one: lowres, else hires (never the full-resolution image)
        images = [k.split("/", 1)[1] for k in stored if k.startswith("images/")]
        lowres = next((k for k in images if "lowres" in k), None) or next((k for k in images if "hires" in k), None)
        wanted = [f"tables/{table_name}", f"shapes/{region}"] + ([f"images/{lowres}"] if lowres else [])
        sdata = ds.to_sdata(elements=wanted)
        gene, symbol = colour_gene(sdata[table_name], panel_of.get(uid))
        section = members.set_index("uid").loc[uid, "title"]
        title = f"{tech}, {section}: {symbol}" if gene else f"{tech}, {section}"
        display(Markdown(f"{tech}, section `{uid}` ({int(members.set_index('uid').loc[uid, 'n_obs'])} observations)."))
        systems = list(get_transformation(sdata.shapes[region], get_all=True))
        system = "global" if "global" in systems else systems[0]
        # spatialdata-plot 0.4.2 rasterises shapes with datashader in "global" only; elsewhere use matplotlib
        method = None if system == "global" else "matplotlib"
        # fmt: off
        if lowres:
            (
                sdata
                .pl.render_images(lowres)
                .pl.render_shapes(region, color=gene, table_name=table_name, method=method, colorbar_params={"label": symbol})
                .pl.show(
                    coordinate_systems=system,
                    title=title,
                )
            )
        else:
            (
                sdata
                .pl.render_shapes(region, color=gene, table_name=table_name, method=method, colorbar_params={"label": symbol})
                .pl.show(
                    coordinate_systems=system,
                    title=title,
                )
            )
        # fmt: on

2026-11-04 04:21:06 | [INFO] root_attr: version

2026-11-04 04:21:06 | [INFO] root_attr: multiscales

2026-11-04 04:21:06 | [INFO] root_attr: omero

2026-11-04 04:21:06 | [INFO] datasets [{'path': 's0', 'coordinateTransformations': [{'type': 'scale', 'scale': [1.0, 1.0, 1.0]}, {'type': 'translation', 'translation': [0.0, 0.0, 0.0]}]}]

2026-11-04 04:21:06 | [INFO] resolution: s0

2026-11-04 04:21:06 | [INFO]  - shape ('c', 'y', 'x') = (3, 421, 600)

2026-11-04 04:21:06 | [INFO]  - chunks =  ['3', '421', '600']

2026-11-04 04:21:06 | [INFO]  - dtype = uint8

2026-11-04 04:21:06 | [INFO] root_attr: version

2026-11-04 04:21:06 | [INFO] root_attr: multiscales

show_call(
    qc_floors, "Gives each technology's count floor, the min_counts of its QC group, for gate 3 and the joint view."
)
floors = qc_floors(thresholds)
show_table(
    pd.DataFrame(
        {
            "value": [str(x) for x in (MIN_SHARED_GENES, MAX_FLOOR_LOSS, N_HVG, N_PCS, N_NEIGHBORS, SEED)],
            "meaning": [
                "genes a view needs before a family is split per technology, and the joint view's gate 1",
                "largest share of a technology's observations that may fall below its floor on the shared genes (gate 3)",
                "highly variable genes kept when a view has more genes",
                "principal components",
                "neighbours per observation",
                "random seed",
            ],
        },
        index=["MIN_SHARED_GENES", "MAX_FLOOR_LOSS", "N_HVG", "N_PCS", "N_NEIGHBORS", "SEED"],
    )
)
show_table(pd.Series(floors, name="count floor").rename_axis("technology").to_frame())

show_call(
    embed_views,
    "Builds one view per family (or per technology when a family shares too few genes), plus the joint view "
    "when the gates pass, and embeds each with PCA, UMAP and Harmony on the section.",
)
views = embed_views(adatas, marker_set, floor=floors, use_gpu=use_gpu)
for name, v in views.items():
    p = v.uns["primer"]
    print(f"{name}: {p['n_obs']} observations, {p['n_sections']} sections; {p['genes']['rule']}. {p['note']}".strip())
first = next(iter(views.values())).uns["primer"]
for uid, reason in first["members_dropped_panel"].items():
    print(f"member {uid} {reason}")
for c in {tuple(c["uids"]): c for v in views.values() for c in v.uns["primer"]["panel_costs"]}.values():
    print(f"{c['technology']} panel of {', '.join(c['uids'])}: {c['reason']} shared genes")
for name, reason in first["views_skipped"].items():
    print(f"no view {name}: {reason}")
gates = next((v.uns["primer"]["gates"] for v in views.values() if "gates" in v.uns["primer"]), None)
if gates is not None:
    print(f"Joint view: {'yes' if gates['joint'] else 'no'}. {gates['reason']}".strip())
    for drop in gates["drops"]:
        print(f"  left out {drop['technology']} ({drop['step']}): {drop['reason']}")

spot: 579835 observations, 154 sections; 2000 highly variable genes of the 17190 measured, expressed genes (seurat_v3 on raw counts, without batch, because the per-section fit failed); its 7 gene panels have 17,922 to 37,082 genes each and share 17,191.
imaging_cell:Xenium: 8461782 observations, 71 sections; all 60 measured, expressed genes of the view; its 10 gene panels have 289 to 5,001 genes each and share 60. Xenium 60 ∩ Merscope 297 = 37 shared genes (under 200), so one view per technology
imaging_cell:Merscope: 318111 observations, 4 sections; all 297 measured, expressed genes of the view; its 2 gene panels have 499 to 500 genes each and share 297. Xenium 60 ∩ Merscope 297 = 37 shared genes (under 200), so one view per technology
member rW0mdZCUUxvVeaK20000 dropped: panel shares 39 genes with the rest of its technology
member yBAAm3DCafRcJ6zN0000 dropped: panel shares 39 genes with the rest of its technology
member vrrCn95ib2xLcGW40000 dropped: panel shares 39 genes with the rest of its technology
member 3RnmyFHYa8XiOHMf0000 dropped: panel shares 39 genes with the rest of its technology
Xenium panel of odTGrRWzlSkvdoF90000, CoPhJeIxAje3zWLd0000, Z1iFoCSHa9kIQlGC0000, DGf3EoxnUcn7FnAY0000: kept: shares 60 of 67 shared genes
Xenium panel of 1rooDV3jrEYzgsdz0000, tS2jpxnC9NS4FKBs0001: kept: shares 60 of 88 shared genes
Xenium panel of t0s2xlrQ6csETvQS0000: kept: shares 60 of 67 shared genes
Xenium panel of RdG9XtWiTzkUJb9g0000, zAhZa4MmArNE4xvy0000, QUo1w4FXBqDC9BAR0001, w28dHAmye37e87JR0000, q5sSLXkFXQjFqFXP0000, EWlXtev09gDq01CT0000, gUeCC5zUfJ89AOab0000, ILpAj8WcGQiH9IMQ0000, qQOOJICzrgozAwd30000, tfpzEbiNvSJLeLL40000, DEkH4oCkISbaSDB30000, kUw3Nw0otoLIWvFf0000, ig5lXDmSfZCU5dVh0000, cOjo8Fe6mVpB4Mu00000: kept: shares 60 of 66 shared genes
Joint view: no. fewer than two families keep a technology
  left out Visium (gate3): loses 100% of its observations below 500 counts on the 37 shared genes (limit 50%)
  left out VisiumHD (gate3): loses 100% of its observations below 200 counts on the 37 shared genes (limit 50%)
  left out Xenium (gate3): loses 81% of its observations below 10 counts on the 37 shared genes (limit 50%)
  left out Merscope (gate3): loses 98% of its observations below 10 counts on the 37 shared genes (limit 50%)

for name, v in views.items():
    colours = [c for c in ("paper", "technology", "tissue") if c in v.obs and v.obs[c].nunique() > 1]
    display(Markdown(f"View `{name}`, before Harmony."))
    sc.pl.umap(v, color=colours or None, ncols=1)
    if "X_umap_harmony" in v.obsm:
        display(Markdown(f"View `{name}`, after Harmony on the section."))
        sc.pl.embedding(v, basis="X_umap_harmony", color=colours or None, ncols=1)
    else:
        display(Markdown(f"View `{name}` has one section, so it has no Harmony embedding."))

show_call(
    mixing_by_view,
    "Measures how evenly the batches mix among each observation's nearest neighbours, per view, before and "
    "after Harmony; the bundle's manifest records the same numbers.",
)
mixing = mixing_by_view(views)
rows, definitions = [], {}
for name, by_batch in mixing.items():
    for label, by_basis in by_batch.items():
        for basis, m in by_basis.items():
            for stat in ("entropy", "other_share"):
                definitions.setdefault(stat, m[stat]["definition"].replace(m["batch_key"], "batch"))
                rows.append(
                    {
                        "view": name,
                        "batch": label,
                        "embedding": "after Harmony" if basis == "X_pca_harmony" else "before Harmony",
                        "statistic": stat,
                        "mean": m[stat]["mean"],
                        "random reference": m[stat]["random"],
                        "limit from proportions": m[stat]["random_limit"],
                        "note": m["note"],
                    }
                )
show_table(pd.DataFrame(rows))
for stat, text in definitions.items():
    print(f"{stat}: {text}")
print(CAVEAT)

entropy: How evenly the batch labels are spread among each observation's 30 nearest neighbours, divided by log(min(30, number of labels)): 0 when all neighbours share one label, 1 when the labels are as evenly spread as the neighbours allow.
other_share: The share of each observation's 30 nearest neighbours that carry another batch.
Well mixed is not the same as biologically right: a correction can also mix cells that differ for a real reason, such as tissue or disease.

fallback_genes = {}
for entries in (marker_set.get("fallback") or {}).values():
    for comp, fb in (entries or {}).items():
        fallback_genes.setdefault(comp, set()).update(str(g).upper() for g in (fb or {}).get("genes") or [])
compartment_genes = {
    comp: {str(g).upper() for g in genes} | fallback_genes.get(comp, set())
    for comp, genes in marker_set["compartments"].items()
}
readings = {}
for name, v in views.items():
    key = "nn_harmony" if "nn_harmony_distances" in v.obsp else "nn_pca"
    leiden = {"resolution": 0.5, "random_state": SEED, "neighbors_key": key, "key_added": "leiden"}
    if use_gpu:  # rapids-singlecell's n_iterations counts aggregation levels, so its default is kept
        rsc.tl.leiden(v, **leiden)
    else:
        sc.tl.leiden(v, flavor="igraph", n_iterations=2, **leiden)
    expr = v.obsm["marker_expr"]
    m = ad.AnnData(
        X=expr.to_numpy(dtype=np.float32),
        obs=v.obs[["leiden"]].copy(),
        var=pd.DataFrame(index=expr.columns.astype(str)),
    )
    display(Markdown(f"View `{name}`: {m.obs['leiden'].nunique()} clusters on `{key}`, {m.n_vars} marker genes."))
    if m.n_vars == 0 or m.obs["leiden"].nunique() < 2:
        display(Markdown("This view has fewer than two clusters or no marker gene, so there is no dotplot."))
        continue
    by_compartment = {
        comp: [g for g in expr.columns if g.upper() in genes] for comp, genes in compartment_genes.items()
    }
    sc.pl.dotplot(m, {k: g for k, g in by_compartment.items() if g}, groupby="leiden", standard_scale="var")
    scores = (
        pd.DataFrame(
            {
                comp: expr[[g for g in expr.columns if g.upper() in genes]].mean(axis=1)
                for comp, genes in compartment_genes.items()
                if any(g.upper() in genes for g in expr.columns)
            }
        )
        .groupby(m.obs["leiden"].to_numpy())
        .mean()
    )
    ranked = scores.apply(lambda r: r.sort_values(ascending=False), axis=1)
    table = pd.DataFrame(
        {
            "observations": m.obs["leiden"].value_counts().reindex(scores.index),
            "top compartment": scores.idxmax(axis=1),
            "top mean": scores.max(axis=1),
            "runner-up": [r.dropna().index[1] if r.notna().sum() > 1 else "" for _, r in ranked.iterrows()],
            "runner-up mean": [r.dropna().iloc[1] if r.notna().sum() > 1 else np.nan for _, r in ranked.iterrows()],
        }
    ).rename_axis("cluster")
    readings[name] = table
    show_table(table)

UserWarning: Over 500 categories found. Plot would be very large.

for name, table in readings.items():
    lines = [f"View `{name}`:"]
    for cluster, r in table.iterrows():
        text = (
            f"- Cluster {cluster} ({int(r['observations'])} observations) has its highest mean marker expression "
            f"for {r['top compartment']} ({r['top mean']:.2f})"
        )
        if r["runner-up"]:
            text += f", followed by {r['runner-up']} ({r['runner-up mean']:.2f})"
        lines.append(text + ".")
    display(Markdown("\n".join(lines)))

show_call(
    export_bundle,
    "Writes the manifest, umap.parquet, qc.parquet, cohort.parquet and embedding.h5ad for this cohort. A cohort whose embedding is over 1 GiB gets parts of whole datasets instead (embedding-1.h5ad, embedding-2.h5ad, ...), and anndata.concat over the parts gives the whole cohort.",
)
bundle = export_bundle(
    views,
    cohort,
    qc,
    spec,
    out_dir,
    release,
    runtime_s=time.time() - start,
    mixing=mixing,
    primers_ref=primers_ref,
)
for f in sorted(bundle.iterdir()):
    print(f"{f.relative_to(out_dir)}  {f.stat().st_size} bytes")
manifest = json.loads((bundle / "manifest.json").read_text())
for key in (
    "slug",
    "release",
    "primers_ref",
    "n_datasets",
    "n_papers",
    "technologies",
    "n_cells_before",
    "n_cells_after",
    "n_flagged",
    "path",
    "hardware",
    "runtime_s",
):
    print(f"{key}: {manifest[key]}")

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.

command = (
    f"papermill primers/notebooks/primer.ipynb primer.ipynb -p slug {slug} -p release {release} "
    f"-p primers_ref {primers_ref} -p use_gpu {use_gpu}" + (f" -p smoke {smoke}" if smoke else "")
)
versions = "\n".join(f"  - {k} {v}" for k, v in manifest["versions"].items())
display(
    Markdown(
        f"- Catalog release: `{release}`\n- primers reference: `{primers_ref}`\n- Versions:\n{versions}\n"
        f"- Hardware: {manifest['hardware']} ({platform.system()})\n- Path: {manifest['path']}\n\n"
        "To run this primer again:"
    )
)
show_command(command)

papermill primers/notebooks/primer.ipynb primer.ipynb -p slug human-lung -p release 2026-11-01 -p primers_ref 4d7af340849c00077f9d21c5b54af9e689a4a070 -p use_gpu True

	uid	title	technology	family	paper	paper_source	licence	n_obs	status
0	6JISVB3uFHGCxvmg0000	GSM8087031	Visium	spot	10.1183/13993003.01114-2024	doi	ncbi-geo-policy	4332	kept
1	k5Q3YIsLRHtzFqFR0000	visium-hd-cytassist-gene-expression-human-lung-cancer-post-xenium-expt	VisiumHD	spot	10x Genomics demo datasets	10x_demo	CC-BY-4.0	448109	kept
2	Wzknw691YzwJrG7m0000	human-lung-cancer-ffpe-2-standard	Visium	spot	10x Genomics demo datasets	10x_demo	CC-BY-4.0	3858	kept
3	PntNMJy70YLUP8eq0000	human-lung-cancer-11-mm-capture-area-ffpe-2-standard	Visium	spot	10x Genomics demo datasets	10x_demo	CC-BY-4.0	6195	kept
4	Dqvfg6twfJ6eEr7o0000	GSM8217625	Visium	spot	10.1038/s41467-024-50281-5	doi	ncbi-geo-policy	4002	kept
...	...	...	...	...	...	...	...	...	...
234	ig5lXDmSfZCU5dVh0000	GSM9509139	Xenium	imaging_cell	10.1038/s41592-026-03089-8	doi	ncbi-geo-policy	78955	kept
235	IlKfWuMs2VIqYXNL0000	GSM9714922	Xenium	imaging_cell	geo-gsm9714922	study_id	ncbi-geo-policy	89890	kept
236	Es8pHq6T7BY9N3GC0000	GSM9714927	Xenium	imaging_cell	geo-gsm9714927	study_id	ncbi-geo-policy	71688	kept
237	0qMt4mCUCzGXE16U0000	GSM9290279	Xenium	imaging_cell	geo-gsm9290279	study_id	ncbi-geo-policy	86009	kept
238	cOjo8Fe6mVpB4Mu00000	GSM9509137	Xenium	imaging_cell	10.1038/s41592-026-03089-8	doi	ncbi-geo-policy	64872	kept

	markers
compartment
at1	AGER, RTKN2, CLIC5
at2	SFTPC, SFTPB, NAPSA
airway	SCGB1A1, FOXJ1, KRT5
fibroblast	COL1A2, DCN, PDGFRA
smooth_muscle	ACTA2, DES, MYH11
endothelial	PECAM1, CLDN5, VWF
t_cell	CD3E, CD3D, CD2
b_plasma	MS4A1, CD79A, MZB1
myeloid	MARCO, CD68, LYZ

	panel	compartment	genes	note
0	Xenium-377-9193d582	at2	LAMP3	AT2 gene on this panel
1	Xenium-392-2248ed9e	at2	SFTPD, ETV5	AT2 genes on this panel; SFTPD is also in club cells
2	Xenium-289-b8cfbe0f	at2	SFTPD, ETV5	AT2 genes on this panel; SFTPD is also in club cells
3	Xenium-5001-079a3380	at2	LAMP3, ABCA3, LPCAT1, ETV5	AT2 genes on this panel
4	Xenium-5001-079a3380	smooth_muscle	MYLK, SYNPO2	smooth-muscle genes on this panel; MYLK is also in myofibroblasts

	count floor
technology
Visium	500.0
SpatialTranscriptomics	500.0
VisiumHD	200.0
Xenium	10.0
Merscope	10.0
MERFISH	10.0

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	16642	fibroblast	0.717	at1	0.099
1	34952	fibroblast	0.564	at1	0.102
10	44683	airway	1.062	at1	0.117
11	574	fibroblast	0.596	at1	0.12
12	188	fibroblast	0.67	at1	0.17
2	65818	fibroblast	0.964	at1	0.134
3	58103	fibroblast	1.223	at1	0.14
4	65417	at2	1.359	at1	0.545
5	137195	fibroblast	1.084	at1	0.286
6	45989	myeloid	1.279	at1	0.3
7	46055	fibroblast	1.185	at1	0.217
8	63613	at2	0.782	at1	0.114
9	606	fibroblast	1.073	at1	0.165

Cohort primer¶

Human lung in health, fibrosis and cancer¶

Setup¶

Query and cohort¶

Quality control¶

Marker set¶

Spatial overview¶

Embedding¶

Biology¶

Draft reading, to be reviewed by the Lead¶

Export and reproduce¶

	status	members
0	kept	233
1	dropped: no paper (no study_id, series or collection)	6

	technologies	min_counts	min_genes	max_pct_mt	flag: min median counts	flag: min share kept	reason
group
spot	Visium, SpatialTranscriptomics	500	250	30.0	1000	0.5	A Visium or ST spot holds about 1 to 10 cells, so a spot with under 500 UMIs or 250 genes is mostly background or tissue edge; solid tissues often reach 20 to 30% mitochondrial reads, so the cap is 30%.
binned_spot	VisiumHD	200	50	30.0	200	0.5	On the 5 Visium HD sections of oliveira-2024-colorectal, 16 um bins at 200 counts kept 74% of bins with the same marker signal as 8 um bins at Oliveira et al.'s 100 (Nat Genet 57:1512-1523), which kept 60%, and found T cells in all 5 sections against 3 (web tasks/2026-11-01-t094-hd-bin-size.md).
imaging_cell	Xenium, Merscope, MERFISH	10	3	n/a	30	0.5	Targeted panels give tens to hundreds of decoded transcripts per segmented cell, and a cell with under 10 transcripts or 3 genes cannot be given a type; the panels have no mitochondrial genes.

	dataset_uid	technology	n_before	n_after	median_counts	median_genes	pct_mt_median	flagged	flag_reason	status	family
0	0BThpfM4wTCZN9Z80000	Visium	4558	4465	29990.0	7807.0	0.0	False		kept	spot
1	0L44re6UyAsJJlNe0000	Visium	660	595	4723.5	2966.5	0.0	False		kept	spot
2	0Z2OCLl5KSlpIn1h0000	Visium	3525	3477	10110.0	4418.0	0.0	False		kept	spot
3	0vD8xbza2w1Eivbg0000	Visium	758	710	2805.5	1732.5	2.134	False		kept	spot
4	1Wqo5HzF9UddVxQu0000	Visium	4613	4409	7893.0	3386.0	0.0	False		kept	spot
...	...	...	...	...	...	...	...	...	...	...	...
228	wb27suQvtPl5ipru0000	Xenium	37429	18988	10.0	9.0	n/a	True	median counts 10 < 30	kept	imaging_cell
229	y763wU2uypqXfH4z0000	Xenium	47611	43641	46.0	27.0	n/a	False		kept	imaging_cell
230	yBAAm3DCafRcJ6zN0000	Xenium	37300	37132	119.0	52.0	n/a	False		kept	imaging_cell
231	zAhZa4MmArNE4xvy0000	Xenium	244374	239946	202.0	168.0	n/a	False		kept	imaging_cell
232	zU8xEuiIiOh0uPCY0001	Xenium	278659	268072	58.0	30.0	n/a	False		kept	imaging_cell

	value	meaning
MIN_SHARED_GENES	200	genes a view needs before a family is split per technology, and the joint view's gate 1
MAX_FLOOR_LOSS	0.5	largest share of a technology's observations that may fall below its floor on the shared genes (gate 3)
N_HVG	2000	highly variable genes kept when a view has more genes
N_PCS	30	principal components
N_NEIGHBORS	30	neighbours per observation
SEED	0	random seed

	view	batch	embedding	statistic	mean	random reference	limit from proportions
0	spot	section	before Harmony	entropy	0.14	0.868	1.288
1	spot	section	before Harmony	other_share	0.217	0.932	0.964
2	spot	section	after Harmony	entropy	0.559	0.868	1.288
3	spot	section	after Harmony	other_share	0.736	0.932	0.964
4	spot	paper	before Harmony	entropy	0.051	0.748	0.857
5	spot	paper	before Harmony	other_share	0.077	0.875	0.905
6	spot	paper	after Harmony	entropy	0.334	0.748	0.857
7	spot	paper	after Harmony	other_share	0.516	0.875	0.905
8	imaging_cell:Xenium	section	before Harmony	entropy	0.374	0.883	1.16
9	imaging_cell:Xenium	section	before Harmony	other_share	0.894	0.96	0.976
10	imaging_cell:Xenium	section	after Harmony	entropy	0.346	0.88	1.16
11	imaging_cell:Xenium	section	after Harmony	other_share	0.459	0.956	0.976
12	imaging_cell:Xenium	paper	before Harmony	entropy	0.301	0.746	0.911
13	imaging_cell:Xenium	paper	before Harmony	other_share	0.815	0.918	0.935
14	imaging_cell:Xenium	paper	after Harmony	entropy	0.282	0.753	0.911
15	imaging_cell:Xenium	paper	after Harmony	other_share	0.41	0.915	0.935
16	imaging_cell:Merscope	section	before Harmony	entropy	0.464	0.881	0.919
17	imaging_cell:Merscope	section	before Harmony	other_share	0.347	0.673	0.696
18	imaging_cell:Merscope	section	after Harmony	entropy	0.607	0.881	0.919
19	imaging_cell:Merscope	section	after Harmony	other_share	0.467	0.673	0.696
20	imaging_cell:Merscope	paper	before Harmony	entropy	0.392	0.97	0.994
21	imaging_cell:Merscope	paper	before Harmony	other_share	0.157	0.479	0.496
22	imaging_cell:Merscope	paper	after Harmony	entropy	0.582	0.97	0.994
23	imaging_cell:Merscope	paper	after Harmony	other_share	0.268	0.48	0.496

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	2240306	myeloid	0.237	fibroblast	0.141
1	606901	t_cell	0.231	fibroblast	0.049
10	32447	b_plasma	0.195	fibroblast	0.033
100	3141	myeloid	0.002	fibroblast	0.0
1000	190	myeloid	0.004	fibroblast	0.0
...	...	...	...	...	...
995	23	myeloid	0.49	fibroblast	0.0
996	40	fibroblast	0.007	fibroblast	0.007
997	200	t_cell	0.895	fibroblast	0.0
998	149	myeloid	0.004	fibroblast	0.004
999	650	myeloid	0.899	fibroblast	0.0

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	17882	b_plasma	0.265	b_plasma	0.265
1	4502	smooth_muscle	0.045	b_plasma	0.016
10	33228	endothelial	0.021	b_plasma	0.003
11	71349	endothelial	0.019	b_plasma	0.004
12	12981	at2	0.396	b_plasma	0.004
13	9087	smooth_muscle	0.073	b_plasma	0.006
14	21568	endothelial	0.444	b_plasma	0.006
15	10137	smooth_muscle	0.749	b_plasma	0.005
16	15121	b_plasma	0.205	b_plasma	0.205
17	546	endothelial	0.002	b_plasma	0.0
18	422	smooth_muscle	0.001	b_plasma	0.001
19	1515	smooth_muscle	0.036	b_plasma	0.006
2	11208	myeloid	0.016	b_plasma	0.007
20	1729	myeloid	0.01	b_plasma	0.007
21	287	smooth_muscle	0.004	b_plasma	0.001
22	541	smooth_muscle	0.001	b_plasma	0.0
3	5012	smooth_muscle	0.042	b_plasma	0.011
4	20620	t_cell	0.103	b_plasma	0.021
5	28962	myeloid	0.066	b_plasma	0.01
6	7775	myeloid	0.229	b_plasma	0.005
7	3436	smooth_muscle	0.163	b_plasma	0.009
8	32359	smooth_muscle	0.027	b_plasma	0.005
9	7844	smooth_muscle	0.038	b_plasma	0.009