Skip to content

Commit 6bb5236

Browse files
authored
Merge pull request #76 from theislab/feat/atera-primer
docs(primers): Atera docs notebook on catalog-2026-11-12
2 parents 66dc3b3 + 28ac7c8 commit 6bb5236

15 files changed

Lines changed: 10947 additions & 1 deletion

‎docs/conf.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ def _write_catalog_columns() -> None:
3737
"**/.pixi",
3838
"reference/_catalog_columns.md",
3939
"reference/_cli.md",
40+
"primers/*/STORY_DRAFT.md", # the page text and its sources, read by export_hero.py, not a docs page
4041
]
4142
smartquotes = False # keep straight quotes, and `--flag` in the CLI reference instead of an en dash
4243

‎docs/execute_notebooks.sh‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,8 @@
44
# bash docs/execute_notebooks.sh
55
# 01_find_and_filter downloads the gene index (about 900 MB) and needs about 40 GB of RAM to query it, so run
66
# this on an HPC node rather than a laptop. The Visium atlas is a GPU run with its own environment and is not
7-
# executed here (see docs/use_cases/visium_atlas/pixi.toml).
7+
# executed here (see docs/use_cases/visium_atlas/pixi.toml). Primers run one at a time in their own environments
8+
# with docs/primers/run_primer.sh <slug>.
89
set -euo pipefail
910
cd "$(dirname "$0")/.."
1011
export SDDB_CACHE_DIR="${SDDB_CACHE_DIR:-$(mktemp -d)}"

‎docs/index.md‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,14 @@ use_cases/training_cohort
142142
use_cases/visium_atlas/visium-atlas
143143
```
144144

145+
```{toctree}
146+
:hidden:
147+
:maxdepth: 2
148+
:caption: Primers
149+
150+
primers/index
151+
```
152+
145153
```{toctree}
146154
:hidden:
147155
:maxdepth: 2

‎docs/primers/atera/.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
atera_out/
2+
.pixi/
3+
web_hero/

‎docs/primers/atera/STORY_DRAFT.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# Atera primer: page story
2+
3+
Numbers from the notebook's run on catalog release 2026-11-12 (outputs in atera_out_release, not in git). One paragraph per figure, caveats first. The story list for primer.json is these 5 paragraphs in order (export_hero.py reads them).
4+
5+
1. This is a first look at six Atera sections, not an analysis. The same rules apply to every section, and three settings were chosen on a 100,000-cell test: a quality floor, scaling before the principal components, and a threshold for cell types. Each tissue has one section, and breast has two (one fresh-frozen, FF, and one formalin-fixed, FFPE), so differences between tissues also contain differences between sections and preservation. The cell types are automated and cursory.
6+
7+
2. The six sections hold 2.31 million cells. 275,586 of them (11.9%) are marked low quality and left out of the analysis: cells with very few counts, or with less than a quarter of their section's median counts, a floor that is our own judgement. The kept FF cells have a median of 5,648 to 6,277 counts per cell and the FFPE cells 1,014 to 2,377, with the cervical cancer section lowest.
8+
9+
3. The sections are integrated with Harmony, and Leiden finds 34 clusters among the 2.04 million kept cells. When the UMAP is coloured by section, immune and stromal cells from different sections share regions, while epithelial cells stay apart by section, as expected for tumour and organ epithelium. One cluster of 41,589 cells, mostly from the cervix, still has low depth (median 449 counts) and is mostly left unassigned.
10+
11+
4. Each cell takes the broad class whose marker genes score clearly highest, or stays unassigned. The brain section is mostly oligodendrocytes (31%), neurons (31%) and astrocytes (20%). In the other five sections, epithelial cells are the largest class (20% to 65%), the cervical cancer section has the most T and NK cells (16%), and the kidney and cervix sections have the most unassigned cells (19% and 17%), partly because the marker list has no kidney tubule markers.
12+
13+
5. On the 4 sections that come with 10x's own broad labels, 52% to 85% of cells get the cursory type that matches their 10x label, with unassigned cells counted as disagreements. This is a consistency check, not a benchmark. The executed notebook has every step and the full figures, including each section drawn on the tissue.
14+
15+
## Sources for each statement
16+
- Paragraph 1: the quality, scaling and label rules, and their choice on the 100,000-cell tests, are in primer.py (quality control, embedding and cell type sections) and in qc_sample_check.py and core_qc.py. One section per tissue and two for breast: atera_out_release/cohort.csv. FF and FFPE from preservation_method ("Fresh Frozen" is shown as FF) and the canonical keys.
17+
- Paragraph 2: atera_out_release/qc_per_dataset.csv. cells sum to 2,310,946, low_quality to 275,586 and cells_kept to 2,035,360. Kept median counts: FF 6,277 (breast), 5,648 (brain), 6,186 (colon); FFPE 2,377 (breast), 1,014 (cervix), 1,845 (kidney). Rule: primer.py qc_mark.
18+
- Paragraph 3: the notebook's cluster table (after the Leiden UMAPs): 34 clusters; the shallowest has 41,589 cells, median 449 counts, mostly from the cervix, and its most common label is "unassigned". The section statement comes from looking at fig_umap_dataset.png (after Harmony).
19+
- Paragraph 4: the notebook's crosstab of cell_type by section. Brain: oligodendrocytes 24,494 of 77,867 (31.46%), neurons 31.1%, astrocytes 20.2%. Epithelial 20.0% (cervix) to 64.9% (colon), the largest class in each non-brain section. T and NK cells: cervix 15.9%, the highest. Unassigned: kidney 18.7%, cervix 17.2%.
20+
- Paragraph 5: atera_out_release/cell_types_vs_10x.csv, 0.519 (kidney) to 0.848 (colon) on 4 sections. Each 10x label is matched to the type most of its assigned cells get, so an unassigned cell never counts as agreeing (primer.py, "Our cursory types next to 10x's labels"). The page shows no tissue view, so the page text has no tissue sentences; the notebook keeps its tissue figure (fig_cell_types_space.png).

‎docs/primers/atera/core_qc.py‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
"""Segmentation and control fields of the low-count core vs the rest of the cells, per section.
2+
3+
Core = Leiden 2 of the 100,000-cell "scale" run of qc_sample_check.py (its obs_scale.parquet). Reads the obs of all six Atera
4+
tables of catalog-2026-11-12 through the client and writes core_vs_rest.csv (100k sample) and section_fields.csv (every
5+
cell of each table: the distributions a rule on the full data would act on).
6+
7+
python -I core_qc.py QCT_DIR (CPU, from docs/primers/atera; SDDB_CACHE_DIR set)
8+
"""
9+
import sys
10+
from pathlib import Path
11+
12+
import numpy as np
13+
import pandas as pd
14+
import sddb
15+
16+
D = Path(sys.argv[1])
17+
FIELDS = ["cell_area", "nucleus_area", "nucleus_count", "transcript_counts", "control_probe_counts", "control_codeword_counts",
18+
"genomic_control_counts", "unassigned_codeword_counts", "deprecated_codeword_counts"]
19+
o = pd.read_parquet(D / "obs_scale.parquet")
20+
core = o.index[o["leiden"].astype(str) == "2"]
21+
22+
atera = sddb.Catalog(version="2026-11-12").query(technology="10x_atera")
23+
m = atera.members
24+
parts = []
25+
for r in m.itertuples():
26+
ds = next(iter(atera.where(m["uid"] == r.uid)))
27+
ob = ds.to_sdata(elements=["tables/table"]).tables["table"].obs
28+
keep = [c for c in FIELDS + ["segmentation_method", "z_level"] if c in ob]
29+
x = ob[keep].copy()
30+
x.index = ob.index.astype(str) + "-" + r.uid
31+
x["section"] = o["section"].astype(str).reindex(x.index)
32+
x["uid"] = r.uid
33+
print(r.uid, r.tissue, len(x), "missing:", sorted(set(FIELDS + ["segmentation_method", "z_level"]) - set(keep)), flush=True)
34+
parts.append(x)
35+
allc = pd.concat(parts)
36+
ctrl = [c for c in FIELDS if c.endswith("_counts") and c != "transcript_counts" and c in allc]
37+
allc["control_share"] = allc[ctrl].sum(axis=1) / (allc["transcript_counts"] + allc[ctrl].sum(axis=1)).replace(0, np.nan)
38+
allc["unassigned_share"] = allc["unassigned_codeword_counts"] / (allc["transcript_counts"] + allc[ctrl].sum(axis=1)).replace(0, np.nan)
39+
allc["counts_per_um2"] = allc["transcript_counts"] / allc["cell_area"].replace(0, np.nan)
40+
allc["no_nucleus"] = allc["nucleus_count"].fillna(0) == 0
41+
if "z_level" in allc:
42+
allc["z_level"] = pd.to_numeric(allc["z_level"], errors="coerce")
43+
44+
samp = allc.loc[allc.index.intersection(o.index)].copy()
45+
samp["section"] = o.loc[samp.index, "section"].astype(str)
46+
samp["group"] = np.where(samp.index.isin(core), "core", "rest")
47+
num = ["cell_area", "nucleus_area", "transcript_counts", "counts_per_um2", "control_share", "unassigned_share", "z_level"]
48+
t = samp.groupby(["section", "group"]).agg(n=("cell_area", "size"), **{f"median_{c}": (c, "median") for c in num if c in samp},
49+
share_no_nucleus=("no_nucleus", "mean"))
50+
if "segmentation_method" in samp:
51+
seg = pd.crosstab([samp["section"], samp["group"]], samp["segmentation_method"].astype(str), normalize="index").round(3)
52+
t = t.join(seg.add_prefix("seg_"))
53+
t.round(3).to_csv(D / "core_vs_rest.csv")
54+
print(t.round(3).to_string(), flush=True)
55+
56+
sec = allc.groupby("uid").agg(n=("cell_area", "size"), **{f"median_{c}": (c, "median") for c in num if c in allc},
57+
share_no_nucleus=("no_nucleus", "mean"))
58+
if "segmentation_method" in allc:
59+
sec = sec.join(pd.crosstab(allc["uid"], allc["segmentation_method"].astype(str), normalize="index").round(3).add_prefix("seg_"))
60+
sec.round(3).to_csv(D / "section_fields.csv")
61+
print(sec.round(3).to_string(), flush=True)
62+
allc[["uid", "section", "transcript_counts", "cell_area", "nucleus_count", "control_share", "unassigned_share", "counts_per_um2"]].to_parquet(D / "all_cells_fields.parquet")
63+
print("done", flush=True)

‎docs/primers/atera/export_hero.py‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
"""Web hero figures for the Atera collection page, from the primer's outputs.
2+
3+
python export_hero.py ATERA_OUT DEST -> DEST/*.webp (<= 1600 px wide, < 400 KB, white background) and DEST/primer.json
4+
"""
5+
import io
6+
import csv
7+
import json
8+
import re
9+
import sys
10+
from pathlib import Path
11+
12+
import numpy as np
13+
from PIL import Image
14+
15+
CREDIT = "Data: 10x Genomics, CC BY 4.0."
16+
FIGS = [ # (source png, web name, alt text, caption), in the page's order
17+
("fig_qc.png", "atera-qc.webp", "Violin plots of counts and genes per cell for each Atera section",
18+
f"Counts and detected genes per cell in each of the six Atera sections, after the same quality filters. {CREDIT}"),
19+
("fig_umap_dataset.png", "atera-umap-dataset.webp", "UMAP of all cells coloured by section, before and after Harmony",
20+
f"All cells from the six sections in one UMAP, before and after Harmony integration across sections. {CREDIT}"),
21+
("fig_cell_types.png", "atera-cell-types.webp", "UMAP coloured by broad cell type, with the cell-type share per tissue",
22+
f"Broad cell types from marker-gene scores per cluster, an automated and cursory annotation, and their share in each tissue. {CREDIT}"),
23+
("fig_cell_types_space.png", "atera-cell-types-space.webp", "Each Atera section with its cells coloured by broad cell type",
24+
f"The same automated, cursory cell types drawn on each section with spatialdata-plot. {CREDIT}"),
25+
]
26+
27+
28+
def squeeze_blank_rows(im, keep=30):
29+
"""Shorten every run of all-white rows (e.g. between the last panel row and the credit footer) to `keep` px."""
30+
a = np.asarray(im)
31+
blank = (a > 250).all(axis=(1, 2))
32+
rows, run = [], 0
33+
for i, b in enumerate(blank):
34+
run = run + 1 if b else 0
35+
if run <= keep:
36+
rows.append(i)
37+
return Image.fromarray(a[rows])
38+
39+
40+
def save_webp(src, dest, max_w=1600, max_kb=400):
41+
im = Image.open(src).convert("RGBA")
42+
bg = Image.new("RGBA", im.size, "white")
43+
im = Image.alpha_composite(bg, im).convert("RGB")
44+
im = squeeze_blank_rows(im)
45+
if im.width > max_w:
46+
im = im.resize((max_w, round(im.height * max_w / im.width)), Image.LANCZOS)
47+
for q in (90, 80, 70, 60, 50, 40): # highest quality under the size cap
48+
buf = io.BytesIO()
49+
im.save(buf, "WEBP", quality=q, method=6)
50+
if buf.tell() < max_kb * 1024:
51+
break
52+
else:
53+
raise ValueError(f"{src}: still {buf.tell() // 1024} KB at quality 40")
54+
dest.write_bytes(buf.getvalue())
55+
return im.width, im.height, buf.tell()
56+
57+
58+
def main(src_dir, dest):
59+
src_dir, dest = Path(src_dir), Path(dest)
60+
dest.mkdir(parents=True, exist_ok=True)
61+
figures = []
62+
for src, name, alt, cap in FIGS:
63+
w, h, size = save_webp(src_dir / src, dest / name)
64+
figures.append({"src": name, "alt": alt, "caption": cap, "width": w, "height": h})
65+
print(name, w, h, size // 1024, "KB")
66+
summary = ("spatialdata-db holds all six 10x Atera whole-transcriptome sections, about 2.3 million cells, in one format. "
67+
"Here the same quality control, Harmony integration and a cursory cell-type annotation run on all of them on one GPU. "
68+
+ CREDIT)
69+
# quality control per section (qc_per_dataset.csv of the same run): low-quality cells are counted on the page
70+
qc = list(csv.DictReader((src_dir / "qc_per_dataset.csv").open()))
71+
qc = [{"section": f"{r['tissue']} ({r['preservation']})", "cells": int(r["cells"]), "low_quality": int(r["low_quality"]),
72+
"cells_kept": int(r["cells_kept"]), "median_counts": float(r["median_counts"])} for r in qc]
73+
# the page story: the numbered paragraphs of STORY_DRAFT.md next to this script, in order
74+
story = [m.group(1) for m in re.finditer(r"^\d+\. (.+)$", (Path(__file__).with_name("STORY_DRAFT.md")).read_text(), re.M)]
75+
(dest / "primer.json").write_text(json.dumps({"summary": summary, "figures": figures, "qc": qc, "story": story}, indent=1))
76+
77+
78+
if __name__ == "__main__":
79+
main(sys.argv[1], sys.argv[2])

0 commit comments

Comments
 (0)