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-skin"
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-skin
query: {'organism': 'human', 'tissue': ['skin of body', 'zone of skin']}
include_uids: None
exclude_uids: None
families: ['spot', 'imaging']
marker_set: skin
rationale: Inflammatory skin disease, melanoma and squamous cell carcinoma on ST, Visium, Visium HD and Xenium; many studies, small total size.

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)

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

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.")

328 kept, 48 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: 196 datasets, so the 40 largest are shown.

imaging_cell: 132 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 'skin', cited from 10.1126/science.aba6500.

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

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.

2026-11-04 05:58:31 | [INFO] root_attr: version

2026-11-04 05:58:31 | [INFO] root_attr: multiscales

2026-11-04 05:58:31 | [INFO] root_attr: omero

2026-11-04 05:58:31 | [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 05:58:31 | [INFO] resolution: s0

2026-11-04 05:58:31 | [INFO]  - shape ('c', 'y', 'x') = (3, 600, 595)

2026-11-04 05:58:31 | [INFO]  - chunks =  ['3', '600', '595']

2026-11-04 05:58:31 | [INFO]  - dtype = uint8

2026-11-04 05:58:31 | [INFO] root_attr: version

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: 338120 observations, 196 sections; 2000 highly variable genes of the 11476 measured, expressed genes (seurat_v3 on raw counts, without batch, because the per-section fit failed); its 22 gene panels have 14,789 to 62,710 genes each and share 11,476.
imaging_cell:Xenium: 284733 observations, 13 sections; all 155 measured, expressed genes of the view; its 3 gene panels have 360 to 5,101 genes each and share 155. Xenium 155 ∩ Merscope 375 = 64 shared genes (under 200), so one view per technology
imaging_cell:Merscope: 1213541 observations, 115 sections; all 375 measured, expressed genes of the view; its 29 gene panels have 480 to 500 genes each and share 375. Xenium 155 ∩ Merscope 375 = 64 shared genes (under 200), so one view per technology
member LCAXoD0ctWEYb92Y0001 dropped: panel shares 56 genes with the rest of its technology
member yY3XV9HxnHqU3ZC30001 dropped: panel shares 56 genes with the rest of its technology
member 26Fo6iDsILpxQ0MX0000 dropped: panel would reduce 155 shared genes to 61
member gHOpZTOsGPHECkuK0000 dropped: panel would reduce 155 shared genes to 61
Xenium panel of 26Fo6iDsILpxQ0MX0000, gHOpZTOsGPHECkuK0000: dropped: panel would reduce 155 shared genes to 61 shared genes
Xenium panel of nceJKOdo07khx69z0000: kept: shares 155 of 215 shared genes
Xenium panel of zJI63TcbTlkIJsKY0000, pBLPikFm2hQ3Qgde0000: kept: shares 155 of 229 shared genes
Joint view: no. fewer than two families keep a technology
  left out Visium (gate3): loses 97% of its observations below 500 counts on the 48 shared genes (limit 50%)
  left out SpatialTranscriptomics (gate3): loses 100% of its observations below 500 counts on the 48 shared genes (limit 50%)
  left out VisiumHD (gate3): loses 100% of its observations below 200 counts on the 48 shared genes (limit 50%)
  left out Merscope (gate3): loses 76% of its observations below 10 counts on the 48 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)

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-skin -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	dW1m7GzIvN6c1cmz0000	GSM7657718	Visium	spot	10.1016/j.crmeth.2024.100759	doi	ncbi-geo-policy	4992	kept
1	XgW5WDDDqXSgZoxG0000	GSM9483977	Visium	spot	10.3389/fimmu.2026.1870524	doi	ncbi-geo-policy	870	kept
2	tmXsTSERDNOqO2OM0000	NCBI469	Visium	spot	10.1038/s41467-022-35319-w	doi	CC-BY-NC-SA-4.0	1504	kept
3	LqSr9Rjv3YT3FlvO0000	human-melanoma-if-stained-ffpe-2-standard	Visium	spot	10x Genomics demo datasets	10x_demo	CC-BY-4.0	3458	kept
4	4CtcXsNU1uQrEWuy0000	GSM9338175	Visium	spot	10.1172/jci.insight.198954	doi	CC-BY-4.0	4833	kept
...	...	...	...	...	...	...	...	...	...
371	LScuFf9BHmZETwiU0000	S-BIAD2376:202405111530_D077--6-12--05112024_VMSC02301:region_6	Merscope	imaging_cell	10.1038/s41588-026-02552-8	doi	CC-BY-4.0	7868	kept
372	Ovk3fZisAhLAi0Rz0000	S-BIAD2376:202405081011_D082--4-12--05082024_VMSC02301:region_5	Merscope	imaging_cell	10.1038/s41588-026-02552-8	doi	CC-BY-4.0	9260	kept
373	6TECNrLJiJWOcyS00000	S-BIAD2376:202308311115_D041-43--lastprobes-083123_VMSC02301:region_1	Merscope	imaging_cell	10.1038/s41588-026-02552-8	doi	CC-BY-4.0	78640	kept
374	QPaAs1CwrK5ANxeQ0000	S-BIAD2376:202405081011_D082--4-12--05082024_VMSC02301:region_3	Merscope	imaging_cell	10.1038/s41588-026-02552-8	doi	CC-BY-4.0	9653	kept
375	acHEWjjLFbpxXZWP0000	S-BIAD2376:202404241001_D165--7-12--04242024_VMSC02301:region_1	Merscope	imaging_cell	10.1038/s41588-026-02552-8	doi	CC-BY-4.0	22397	kept

	dataset_uid	technology	n_before	n_after	median_counts	median_genes	pct_mt_median	flagged	flag_reason	status	family
0	02Og8NVAa41JDFhi0000	Visium	1906	1717	3398.5	1582.5	5.244	False		kept	spot
1	1Haj9pIFZeaHFPvN0000	Visium	14335	10909	1886.0	1010.0	0.554	False		kept	spot
2	1Lehzb5WV60DegjI0000	Visium	1117	364	127.0	105.0	4.989	True	median counts 127 < 1000; 33% kept < 50%	kept	spot
3	1Y79jorcVL1O2Onf0000	Visium	598	87	38.0	36.0	4.974	True	median counts 38 < 1000; 15% kept < 50%	kept	spot
4	2KBVVl4aVjRqy8Ik0000	Visium	825	684	4854.0	1816.0	4.411	False		kept	spot
...	...	...	...	...	...	...	...	...	...	...	...
323	yG8s7Y0LO6Clk4pt0000	Merscope	7154	4880	17.0	12.0	n/a	True	median counts 17 < 30	kept	imaging_cell
324	ySHxtZGcFaEpbWoE0000	Merscope	10569	6319	12.0	10.0	n/a	True	median counts 12 < 30	kept	imaging_cell
325	yY3XV9HxnHqU3ZC30001	Xenium	68476	46019	21.0	13.0	n/a	True	median counts 21 < 30	kept	imaging_cell
326	yr83FMrgCggrRR5W0000	Merscope	10783	7288	18.0	13.0	n/a	True	median counts 18 < 30	kept	imaging_cell
327	zJI63TcbTlkIJsKY0000	Xenium	64091	61027	90.0	42.0	n/a	False		kept	imaging_cell

	markers
compartment
basal_keratinocyte	KRT14, KRT5
differentiated_keratinocyte	KRT10, KRT1
melanocyte	PMEL, TYRP1, DCT, MLANA
fibroblast	COL1A1, DCN, PDGFRA
smooth_muscle	ACTA2, DES, MYH11
endothelial	PECAM1, VWF
t_cell	CD3E, CD3D, CD2
myeloid	CD68, CD14, CSF1R
langerhans	CD207, CD1A

	panel	compartment	genes	note
0	Xenium-5101-9ddba058	basal_keratinocyte	KRT15, COL17A1, TP63	basal keratinocyte genes on this panel
1	Xenium-377-9193d582	basal_keratinocyte	COL17A1	basal keratinocyte gene on this panel
2	Xenium-377-9193d582	differentiated_keratinocyte	LY6D, SERPINB3	suprabasal keratinocyte genes on this panel
3	Xenium-360-97551f14	differentiated_keratinocyte	KRTDAP, KRT2	differentiated keratinocyte genes on this panel
4	Xenium-480-18bfc240	differentiated_keratinocyte	not measured	only stress keratins (KRT6A, KRT16) on this panel, which do not mark differentiation
5	Xenium-480-18bfc240	fibroblast	LUM, MMP2	fibroblast genes on this panel

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	41180	melanocyte	1.483	differentiated_keratinocyte	0.134
1	60988	fibroblast	1.847	differentiated_keratinocyte	0.401
10	10327	fibroblast	0.943	differentiated_keratinocyte	0.04
11	7224	basal_keratinocyte	1.211	differentiated_keratinocyte	0.648
12	777	melanocyte	1.214	differentiated_keratinocyte	0.164
13	7042	melanocyte	1.719	differentiated_keratinocyte	0.132
14	914	melanocyte	0.273	differentiated_keratinocyte	0.056
2	52844	differentiated_keratinocyte	2.228	differentiated_keratinocyte	2.228
3	15229	differentiated_keratinocyte	1.731	differentiated_keratinocyte	1.731
4	39318	fibroblast	1.075	differentiated_keratinocyte	0.439
5	46816	fibroblast	1.203	differentiated_keratinocyte	0.678
6	18747	basal_keratinocyte	2.205	differentiated_keratinocyte	1.611
7	20347	basal_keratinocyte	1.601	differentiated_keratinocyte	0.472
8	10825	basal_keratinocyte	1.45	differentiated_keratinocyte	0.623
9	5542	melanocyte	2.179	differentiated_keratinocyte	0.719

Cohort primer¶

Human skin¶

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	328
1	dropped: no paper (no study_id, series or collection)	48

	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.

	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

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

	view	batch	embedding	statistic	mean	random reference	limit from proportions	note
0	spot	section	before Harmony	entropy	0.3	0.934	1.404
1	spot	section	before Harmony	other_share	0.444	0.955	0.988
2	spot	section	after Harmony	entropy	0.586	0.935	1.404
3	spot	section	after Harmony	other_share	0.742	0.955	0.988
4	spot	paper	before Harmony	entropy	0.084	0.766	0.89
5	spot	paper	before Harmony	other_share	0.126	0.893	0.924
6	spot	paper	after Harmony	entropy	0.352	0.766	0.89
7	spot	paper	after Harmony	other_share	0.51	0.893	0.924
8	imaging_cell:Xenium	section	before Harmony	entropy	0.369	0.739	0.823
9	imaging_cell:Xenium	section	before Harmony	other_share	0.502	0.807	0.833
10	imaging_cell:Xenium	section	after Harmony	entropy	0.478	0.737	0.823
11	imaging_cell:Xenium	section	after Harmony	other_share	0.557	0.806	0.833
12	imaging_cell:Xenium	paper	before Harmony	entropy	0.2	0.882	0.912
13	imaging_cell:Xenium	paper	before Harmony	other_share	0.121	0.587	0.606
14	imaging_cell:Xenium	paper	after Harmony	entropy	0.461	0.883	0.912
15	imaging_cell:Xenium	paper	after Harmony	other_share	0.292	0.587	0.606
16	imaging_cell:Merscope	section	before Harmony	entropy	0.631	0.889	1.258
17	imaging_cell:Merscope	section	before Harmony	other_share	0.767	0.943	0.976
18	imaging_cell:Merscope	section	after Harmony	entropy	0.756	0.889	1.258
19	imaging_cell:Merscope	section	after Harmony	other_share	0.878	0.943	0.976
20	imaging_cell:Merscope	paper	before Harmony	entropy	n/a	n/a	n/a	one paper, so neither statistic is defined
21	imaging_cell:Merscope	paper	before Harmony	other_share	n/a	n/a	n/a	one paper, so neither statistic is defined
22	imaging_cell:Merscope	paper	after Harmony	entropy	n/a	n/a	n/a	one paper, so neither statistic is defined
23	imaging_cell:Merscope	paper	after Harmony	other_share	n/a	n/a	n/a	one paper, so neither statistic is defined

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	29717	fibroblast	0.636	differentiated_keratinocyte	0.005
1	56226	differentiated_keratinocyte	0.733	differentiated_keratinocyte	0.733
10	801	fibroblast	0.001	differentiated_keratinocyte	0.0
11	135	fibroblast	0.026	differentiated_keratinocyte	0.0
12	4042	fibroblast	0.162	differentiated_keratinocyte	0.006
13	4266	fibroblast	0.186	differentiated_keratinocyte	0.007
14	5612	differentiated_keratinocyte	0.923	differentiated_keratinocyte	0.923
15	242	myeloid	0.009	differentiated_keratinocyte	0.0
16	312	basal_keratinocyte	0.023	differentiated_keratinocyte	0.0
17	229	fibroblast	0.011	differentiated_keratinocyte	0.0
18	242	fibroblast	0.024	differentiated_keratinocyte	0.0
19	2695	endothelial	0.175	differentiated_keratinocyte	0.002
2	18529	t_cell	0.753	differentiated_keratinocyte	0.008
20	140	basal_keratinocyte	0.005	differentiated_keratinocyte	0.0
21	678	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
22	1849	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
23	882	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
24	149	basal_keratinocyte	0.005	differentiated_keratinocyte	0.0
25	1704	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
26	180	basal_keratinocyte	0.03	differentiated_keratinocyte	0.015
27	134	differentiated_keratinocyte	0.014	differentiated_keratinocyte	0.014
28	527	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
29	1582	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
3	17534	smooth_muscle	0.501	differentiated_keratinocyte	0.005
30	161	differentiated_keratinocyte	0.012	differentiated_keratinocyte	0.012
31	1425	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
32	46	fibroblast	0.033	differentiated_keratinocyte	0.0
33	18	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
34	52	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
35	928	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
36	137	differentiated_keratinocyte	0.089	differentiated_keratinocyte	0.089
37	55	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
38	10	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
39	41	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
4	7620	basal_keratinocyte	0.678	differentiated_keratinocyte	0.035
40	44	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
41	13	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
42	94	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
43	28	fibroblast	0.076	differentiated_keratinocyte	0.0
44	23	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
45	51	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
46	8	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
47	93	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
48	43	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
49	32	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
5	18542	myeloid	0.62	differentiated_keratinocyte	0.045
50	32	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
51	40	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
52	13	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
53	57	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
54	103	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
55	458	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
56	48	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
57	11	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
58	87	basal_keratinocyte	0.0	differentiated_keratinocyte	0.0
6	22745	endothelial	0.559	differentiated_keratinocyte	0.004
7	42051	fibroblast	0.731	differentiated_keratinocyte	0.003
8	1334	fibroblast	0.27	differentiated_keratinocyte	0.005
9	39883	basal_keratinocyte	1.191	differentiated_keratinocyte	0.029

	observations	top compartment	top mean	runner-up	runner-up mean
cluster
0	31058	t_cell	0.1	endothelial	0.032
1	112399	basal_keratinocyte	0.449	endothelial	0.013
10	129909	basal_keratinocyte	0.114	endothelial	0.012
11	153941	basal_keratinocyte	0.104	endothelial	0.016
12	43320	smooth_muscle	0.165	endothelial	0.028
13	58932	smooth_muscle	0.756	endothelial	0.093
14	5301	fibroblast	0.114	endothelial	0.036
2	21369	endothelial	0.287	endothelial	0.287
3	18011	fibroblast	0.089	endothelial	0.046
4	76415	myeloid	0.245	endothelial	0.031
5	56117	basal_keratinocyte	0.027	endothelial	0.012
6	337870	fibroblast	0.284	endothelial	0.019
7	21338	fibroblast	0.079	endothelial	0.022
8	34437	fibroblast	0.043	endothelial	0.04
9	113124	endothelial	0.883	endothelial	0.883