diff --git a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh index b695b0f08..4d4a3aebe 100644 --- a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh +++ b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh @@ -19,7 +19,7 @@ # SCRATCH_DIR=/path/to/big/scratch ./mirror_bruker_to_s3.sh set -euo pipefail - +SCRATCH_DIR="/Volumes/SeagateHHD" # --- Config ----------------------------------------------------------------- SRC_BASE="https://smi-public.objects.liquidweb.services" @@ -33,11 +33,11 @@ SCRATCH_DIR="${SCRATCH_DIR:-$PWD/bruker_mirror_scratch}" # - SOURCE is either a full https URL, or a name relative to SRC_BASE (the mouse/liver host). # The URL-encoded names are what the liquidweb server serves. # - LOCAL_NAME is the (decoded) name to store under on S3. -FILES=( - "HalfBrain.zip|HalfBrain.zip" - "Half%20%20Brain%20simple%20%20files%20.zip|Half Brain simple files.zip" - "NormalLiverFiles.zip|NormalLiverFiles.zip" -) +#FILES=( +# "HalfBrain.zip|HalfBrain.zip" +# "Half%20%20Brain%20simple%20%20files%20.zip|Half Brain simple files.zip" +# "NormalLiverFiles.zip|NormalLiverFiles.zip" +#) # NSCLC lung-cancer samples: each ships a flat-files+cell-labels archive and a # raw-morphology-images archive. The bruker_cosmx_nsclc loader streams both from S3. @@ -132,7 +132,7 @@ for entry in "${FILES[@]}"; do # Upload to S3 echo " Uploading to s3://$bucket/$key ..." - aws s3 cp "$local_path" "s3://$bucket/$key" + aws s3 cp "$local_path" "s3://$bucket/$key" --profile op # Free scratch space before the next (much larger) file echo " Removing local copy to free space" diff --git a/scripts/create_resources/spatial/upload_vizgen_merscope_2d.sh b/scripts/create_resources/spatial/upload_vizgen_merscope_2d.sh new file mode 100644 index 000000000..3fff704fe --- /dev/null +++ b/scripts/create_resources/spatial/upload_vizgen_merscope_2d.sh @@ -0,0 +1,185 @@ +#!/bin/bash + +# Mirror a z3-only (2D) copy of the Vizgen MERSCOPE FFPE showcase raw data from Google +# Cloud Storage to S3, so the Nebius/Seqera benchmark can stage its inputs without any +# Google Cloud credentials. +# +# WHY THIS EXISTS +# process_vizgen_merscope_nebius.sh used to point `input:` directly at +# gs://vz-ffpe-showcase/. That bucket is Vizgen's access-controlled +# data-release bucket (NOT public: an unauthenticated GET returns +# "Anonymous caller does not have storage.objects.get access ..."). The Nebius compute +# env has no GCP credentials, so Nextflow staged the input as an anonymous caller and +# the run died before the loader ran. Mirroring the data to s3://openproblems-data +# (which the Nebius env already reads for every other spatial dataset) removes the +# GCP-auth dependency from every future run. +# +# WHAT IS KEPT (lossless for this 2D pipeline) +# The vizgen_merscope loader calls spatialdata_io.merscope(..., z_layers=3), i.e. it +# only ever reads z-plane 3 of the mosaic images. So we mirror ONLY *_z3.tif and drop +# every *_z{0,1,2,4,5,6}.tif (7 focal planes -> 1). Everything else is copied verbatim: +# - cell_boundaries/*.hdf5 (OLD showcase format; the loader rebuilds +# cell_boundaries.parquet from these at run time) +# - images/mosaic_DAPI_z3.tif (DAPI only, by default — the loader keeps ONLY the DAPI +# channel: `sdata["morphology_mip"].sel(c=["DAPI"])`, so +# dropping PolyT/Cellbound stains is lossless for this +# pipeline and cuts the image bytes to ~1/5. If the +# keep-more-stains TODO in the loader is ever done, re-run +# with KEEP_ONLY_DAPI=0 to mirror all stains at z3.) +# - images/micron_to_mosaic_pixel_transform.csv, images/manifest.json +# - cell_by_gene.csv, cell_metadata.csv, detected_transcripts.csv +# KEEP_ONLY_DAPI defaults to 1 (DAPI-only). Set KEEP_ONLY_DAPI=0 to keep all stains at z3. +# +# PARALLEL +# Each patient ships ~1100-2500 tiny cell_boundaries/*.hdf5 files (~15k total across the 7 +# patients). Streaming them one-at-a-time is dominated by per-file spawn/handshake overhead +# (~15 s/file => ~2-3 DAYS). This script runs the transfers through a pool of $JOBS parallel +# workers (default 16), which hides that overhead and cuts the whole mirror to a few hours. +# Big z3 tifs / detected_transcripts.csv are bandwidth-bound, so a handful running +# concurrently just keeps the uplink saturated. Tune with JOBS= (e.g. 32 for the +# mostly-small-file tail). +# +# REQUIREMENTS (run this on a machine that is authenticated to BOTH clouds) +# - gcloud CLI, authenticated with an account granted access to gs://vz-ffpe-showcase +# (Vizgen data-release-program access — request it via https://info.vizgen.com/ffpe-showcase). +# - aws CLI, with a profile that can write s3://openproblems-data (default: profile "op"). +# The transfer STREAMS each object gcloud->aws (no local disk staging). +# +# IDEMPOTENT + RESUMABLE +# Objects already in S3 with a matching byte size are skipped (one bulk `aws s3 ls` up +# front, so resuming does not cost a HEAD per object). A partially-copied file is aborted +# by aws and reads back as absent, so it is simply re-sent on the next run. Just re-run. +# +# Usage: +# AWS_PROFILE=op ./upload_vizgen_merscope_2d.sh +# # more workers for the small-file tail: +# JOBS=32 AWS_PROFILE=op ./upload_vizgen_merscope_2d.sh +# # or override anything: +# GCS_BASE=gs://vz-ffpe-showcase \ +# S3_DEST=s3://openproblems-data/resources/raw_data/txSim_custom/vizgen_merscope \ +# SAMPLES="HumanBreastCancerPatient1 HumanLungCancerPatient1" \ +# JOBS=24 AWS_PROFILE=op ./upload_vizgen_merscope_2d.sh + +set -euo pipefail + +GCS_BASE="${GCS_BASE:-gs://vz-ffpe-showcase}" +S3_DEST="${S3_DEST:-s3://openproblems-data/resources/raw_data/txSim_custom/vizgen_merscope}" +AWS_PROFILE="${AWS_PROFILE:-op}" +Z="${Z:-3}" # which single z-plane to keep +KEEP_ONLY_DAPI="${KEEP_ONLY_DAPI:-1}" # 1 = DAPI-only (default, lossless for this loader); 0 = all stains at z3 +JOBS="${JOBS:-12}" # number of concurrent gcloud->aws transfers (each = 2 TLS + # conns; >~16 tends to trigger transient TLS/read-timeout + # errors, which the worker retries with backoff anyway) + +# The 7 patients currently active (uncommented) in process_vizgen_merscope_nebius.sh. +# Additional patients (melanoma/ovarian/prostate/uterine) are commented-out there and in +# the s3 sibling process_vizgen_merscope.sh — add their folder names here to mirror them. +SAMPLES=(${SAMPLES:-\ + HumanBreastCancerPatient1 \ + HumanLiverCancerPatient1 \ + HumanLiverCancerPatient2 \ + HumanLungCancerPatient1 \ + HumanLungCancerPatient2 \ + HumanColonCancerPatient1 \ + HumanColonCancerPatient2}) + +export AWS_PROFILE +command -v gcloud >/dev/null || { echo "ERROR: gcloud CLI not found (needed to read gs://vz-ffpe-showcase)" >&2; exit 1; } +command -v aws >/dev/null || { echo "ERROR: aws CLI not found (needed to write $S3_DEST)" >&2; exit 1; } + +gcs_bucket="$(echo "$GCS_BASE" | sed -E 's#^gs://([^/]+).*#\1#')" +s3_bucket="$(echo "$S3_DEST" | sed -E 's#^s3://([^/]+)/.*#\1#')" +s3_prefix="$(echo "$S3_DEST" | sed -E 's#^s3://[^/]+/(.*)#\1#')" +export gcs_bucket s3_bucket s3_prefix + +# Should this object be dropped from the mirror? +drop_object() { + local rel="$1" base="${1##*/}" + # keep only z-plane $Z: drop any *_z.tif (matches mosaic_*_z*.tif and boundaries_z*.tif) + if [[ "$base" =~ _z([0-9]+)\.tif$ && "${BASH_REMATCH[1]}" != "$Z" ]]; then return 0; fi + # optional: drop non-DAPI stains entirely + if [[ "$KEEP_ONLY_DAPI" == "1" && "$base" =~ ^mosaic_ && ! "$base" =~ ^mosaic_DAPI_ ]]; then return 0; fi + return 1 +} + +# Worker: stream one "|" record gcloud->aws. Exported so `xargs -P` can run it in +# parallel `bash -c` subshells. Kept free of bash-4 features so it works on the stock macOS +# /bin/bash 3.2. The already-uploaded files are filtered out up front (see below), so the +# worker just transfers. A failed transfer is logged, not fatal: the batch keeps going and the +# file (absent, since aws aborts a partial multipart) is retried on the next resumable run. +transfer_one() { + local rec="$1" + local size="${rec%%|*}" + local rel="${rec#*|}" + local attempt=0 max=4 + # Retry transient GCS/S3 errors (read timeouts, SSL/TLS handshake failures) that show up under + # high concurrency. pipefail (set by the caller) makes a gcloud-side failure fail the pipe too. + while :; do + attempt=$((attempt+1)) + if gcloud storage cat "gs://$gcs_bucket/$rel" 2>/dev/null \ + | aws s3 cp - "s3://$s3_bucket/$s3_prefix/$rel" --expected-size "$size" --only-show-errors 2>/dev/null; then + printf ' OK %s (%s B)\n' "$rel" "$size" + return 0 + fi + if [ "$attempt" -ge "$max" ]; then + printf ' FAIL %s (after %s attempts; retry on re-run)\n' "$rel" "$attempt" >&2 + return 0 + fi + sleep $((attempt * 3)) # linear backoff: 3s, 6s, 9s + done +} +export -f transfer_one + +present="$(mktemp -t viz_present)" +wanted_raw="$(mktemp -t viz_wanted_raw)" +wanted="$(mktemp -t viz_wanted)" +worklist="$(mktemp -t viz_worklist)" +trap 'rm -f "$present" "$wanted_raw" "$wanted" "$worklist"' EXIT + +# 1) One bulk listing of what is already in S3 -> "\t" (sorted). One API call for the +# whole resume state, instead of a HEAD per object. +echo "Listing existing objects under $S3_DEST ..." +aws s3 ls "s3://$s3_bucket/$s3_prefix/" --recursive 2>/dev/null \ + | awk -v p="$s3_prefix/" 'BEGIN{OFS="\t"} {k=$4; sub("^"p,"",k); if(k!="") print k,$3}' \ + | sort > "$present" +echo " $(wc -l < "$present" | tr -d ' ') objects already in S3." + +# 2) Enumerate GCS, apply the keep-rules -> "\t". Written to a file (not piped) so +# the counters stay in this shell and the progress echoes go to stderr, not into the data. +total=0; dropped=0 +for sample in "${SAMPLES[@]}"; do + echo "Enumerating $GCS_BASE/$sample ..." >&2 + while IFS= read -r line; do + [[ "$line" == *"gs://"* ]] || continue # skip the TOTAL: summary line + url="$(awk '{print $NF}' <<<"$line")" + size="$(awk '{print $1}' <<<"$line")" + [[ "$url" == */ ]] && continue # skip directory placeholders + [[ "$size" =~ ^[0-9]+$ ]] || continue # skip anything without a numeric size + rel="${url#gs://$gcs_bucket/}" + total=$((total+1)) + if drop_object "$rel"; then dropped=$((dropped+1)); continue; fi + printf '%s\t%s\n' "$rel" "$size" >> "$wanted_raw" + done < <(gcloud storage ls -l "$GCS_BASE/$sample/**") +done +sort "$wanted_raw" > "$wanted" + +# 3) Work list = wanted objects not already in S3 at the same size (comm on the sorted files). +comm -23 "$wanted" "$present" | awk -F'\t' '{print $2"|"$1}' > "$worklist" +queued="$(wc -l < "$worklist" | tr -d ' ')" +echo "================================================================" +echo "Enumerated $total objects: $dropped dropped, $(( $(wc -l < "$wanted" | tr -d ' ') - queued )) already in S3, $queued to transfer." +echo "Transferring with $JOBS parallel workers ..." + +# 4) Run the transfers in parallel. -I REC feeds one record per worker; -P runs $JOBS at once. +# Use the SAME bash that is running this script ($BASH) so the exported transfer_one +# imports correctly regardless of which bash is first on PATH. +if [[ "$queued" -gt 0 ]]; then + xargs -P "$JOBS" -I REC "${BASH:-bash}" -c 'set -o pipefail; transfer_one "$@"' _ REC < "$worklist" +fi + +echo "================================================================" +echo "Done. objects=$total dropped=$dropped queued=$queued (JOBS=$JOBS)" +echo "Loader --input paths (use these in process_vizgen_merscope_nebius.sh):" +for sample in "${SAMPLES[@]}"; do + echo " $S3_DEST/$sample" +done diff --git a/src/base/labels_nebius.config b/src/base/labels_nebius.config index 93ab8d178..5dc68c20e 100644 --- a/src/base/labels_nebius.config +++ b/src/base/labels_nebius.config @@ -193,6 +193,7 @@ process { withLabel: hightime { time = 12.h } withLabel: veryhightime { time = 24.h } +withLabel: veryveryhightime { time = 48.h } // 2 days (e.g. allen_brain_cell_atlas_merfish loader) // make sure publishstates gets enough disk space and memory withName:'.*publishStatesProc' { diff --git a/src/base/setup_spatialdata_partial.yaml b/src/base/setup_spatialdata_partial.yaml index 7e3f4be48..c887a67f0 100644 --- a/src/base/setup_spatialdata_partial.yaml +++ b/src/base/setup_spatialdata_partial.yaml @@ -1,3 +1,3 @@ setup: - type: python - pypi: ["spatialdata>=0.7.3", "anndata>=0.12.0", "zarr>=3.0.0"] + pypi: ["spatialdata>=0.7.3", "anndata>=0.12.0,<0.13", "zarr>=3.0.0"] diff --git a/src/datasets/loaders/allen_brain_cell_atlas_merfish/config.vsh.yaml b/src/datasets/loaders/allen_brain_cell_atlas_merfish/config.vsh.yaml index 700616235..0fe590ecf 100644 --- a/src/datasets/loaders/allen_brain_cell_atlas_merfish/config.vsh.yaml +++ b/src/datasets/loaders/allen_brain_cell_atlas_merfish/config.vsh.yaml @@ -91,4 +91,4 @@ runners: - type: executable - type: nextflow directives: - label: [highmem, midcpu, hightime] + label: [veryhighmem, midcpu, veryveryhightime] diff --git a/src/datasets/loaders/bruker_cosmx/NOTES.md b/src/datasets/loaders/bruker_cosmx/NOTES.md new file mode 100644 index 000000000..e970957dc --- /dev/null +++ b/src/datasets/loaders/bruker_cosmx/NOTES.md @@ -0,0 +1,255 @@ +# bruker_cosmx loader — notes + +## What this component is + +Dataset loader (`datasets/loaders` namespace) that converts a **Bruker / NanoString +CosMx** export into the standardized `raw_ist` SpatialData zarr consumed by the rest of the +pipeline (API: `/src/api/file_common_ist.yaml`, merged as `--output` in `config.vsh.yaml`). + +It is a thin wrapper around **`sopa.io.cosmx`** (the sopa dependency is unpinned in +`config.vsh.yaml` → currently sopa ~2.2.x). sopa does the heavy lifting: it stitches the +per-FOV morphology images and cell-label tifs into one global image/label, reads the +transcripts flat file, and builds the cells table. Our script's job is (1) massage the +on-disk directory layout into the shape sopa expects, (2) rename sopa's element keys to our +API names, (3) drop non-DAPI image channels, (4) attach dataset metadata, (5) write. + +Some CosMx exports ship **no AtoMx flat files at all** (the liver "Raw Data Files" release — +see Version 4 below). For those the script first **reconstructs** the flat files sopa needs +from the raw per-FOV decoding CSVs, then feeds the same `sopa.io.cosmx` call. This is +detected automatically (no config knob): if `input_flat_files` is unset AND the raw zip has +no `*_tx_file.csv`, reconstruction mode is on. + +Sibling loader `bruker_cosmx_nsclc` handles the CosMx lung-cancer layout (per-sample +subdirs, per-z-plane images) — that's a *different* export structure (see the docstring +"Version 3"), not covered here. + +Datasets driven through this loader (see `scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh`): +- `bruker_cosmx/bruker_mouse_brain_cosmx/rep1` — full mouse-brain hemisphere (`HalfBrain.zip` + + separate `Half Brain simple files.zip`). **The big one** — the source of the OOM this + NOTES documents. +- `bruker_cosmx/bruker_human_liver_cosmx` — human liver (`NormalLiverFiles.zip`). **No flat + files exist for this dataset** (NanoString only released raw images + TileDB + Seurat), so + it runs through the **reconstruction** path (Version 4 / see section below). `input_flat_files` + is left unset in the run script. + +## Supported input layouts + +The script's module docstring (`script.py:3-80`) is the authoritative catalogue of the three +CosMx export shapes seen in the wild. The two this loader handles: + +- **Flat files bundled in the raw zip** (liver): `input_flat_files` is left unset. +- **Flat files in a separate zip** (mouse brain): `input_flat_files` is set, and the CSVs are + symlinked into `CellStatsDir` (`script.py:170-184`). + +In both, sopa wants a `CellStatsDir/CellLabels/` folder, but many exports only ship the +per-FOV `CellLabels_F###.tif` inside `FOV*/` folders. `script.py:186-200` gathers those into +a `CellLabels/` folder via symlinks (see sopa issue #285). + +## script.py — step by step + +1. **Extract** `input_raw` (strip the single wrapper dir) → find `CellStatsDir` + (`script.py:155-168`). +2. **Extract flat files** if provided, symlink the `*.csv` into `CellStatsDir` + (`script.py:170-184`). +3. **Assemble `CellLabels/`** from per-FOV tifs if needed (`script.py:186-200`). +4. **Lean transcript-read shim** installed here (`script.py:202-253`) — see next section. +5. **`sopa.io.cosmx(...)`** with `cells_labels/table/polygons=True`, `read_proteins=False` + (`script.py:261`). Stitched image/labels come back as **lazy dask** — they are *not* the + memory peak; they're written chunked. +6. **Rename element keys** to API names (`script.py:277-288`): + `stitched_image→morphology_mip`, `stitched_labels→cell_labels`, `points→transcripts`, + `cells_polygons→cell_boundaries`. +7. **Add API transcript columns** (`script.py:290-292`): `cell_id = global_cell_id`, + `feature_name = target` (dask-lazy copies; cheap now that `target` is categorical). +8. **Keep only the DNA channel** of the morphology image (`script.py:300`, + `.sel(c=["DNA"])`) — we treat the "DNA" stain as DAPI. TODO in code: keep PolyT / + Cellbound channels one day (currently saving/plotting fails when they're kept). +9. **Attach dataset metadata** to `table.uns` (`script.py:303-307`). +10. **`sdata.write(par["output"])`** (`script.py:315`). + +## Flat-file reconstruction for raw-only exports (Version 4, liver — 2026-07-29) + +`NormalLiverFiles.zip` contains only imaging + segmentation + raw decoding output; the AtoMx +"flat files" were never released for the liver dataset (verified: no `*_tx_file.csv` / +`*_fov_positions_file.csv` anywhere in the zip, and no companion "simple files" zip on S3 — +only TileDB/Seurat). `sopa.io.cosmx` therefore died at `_infer_dataset_id` +(`ValueError: Could not infer dataset_id ...`). The raw material to rebuild the flat files is +present, so the loader now reconstructs them. + +**Detection** (`script.py`, `RECONSTRUCT`): on when `input_flat_files` is unset *and* the raw +zip has no `*_tx_file.csv` (peeked cheaply from the zip central directory via +`_zip_has_flat_files`). Off for mouse brain (separate flat-files zip) and any export that +ships flat files inside the raw zip. + +**Extraction** switches from the denylist to an *allowlist* (`_member_wanted_reconstruct` / +`_RECON_KEEP`) that keeps only what reconstruction + sopa need: per-FOV `CellLabels_F*.tif`, +`*Cell_Stats_F*.csv`, `Morphology2D/*.TIF`, `AnalysisResults/**/*complete_code_cell_target_call_coord.csv`, +and `RunSummary/*.fovs.csv`. Footprint for the liver hemisphere ≈ **61 GB** (vs 440 GB of +Morphology3D alone). Crucially this keeps `AnalysisResults`/`RunSummary`, which the flat-file +denylist drops. + +**`reconstruct_flat_files()`** writes four flat files into `CellStatsDir` (`DATA_DIR`): +- `{id}_fov_positions_file.csv` ← `RunSummary/latest.fovs.csv` (headerless; `x_mm`=col 1, + `y_mm`=col 2, `fov`=col 6 — matches the Metrics4D `Xpos_um`/`Ypos_um`). +- `{id}_tx_file.csv` ← concat of the 304 `complete_code_cell_target_call_coord.csv` + (`fov, cell_ID=CellId, target, x_global_px, y_global_px, z`), written FOV-by-FOV so only + one FOV is in RAM (≈229 M transcripts total for liver). +- `{id}_metadata_file.csv` ← `Cell_Stats_F*.csv` (`fov, cell_ID, CenterX/Y_global_px, Area`). +- `{id}_exprMat_file.csv` ← per-cell gene counts crosstab'd from the assigned transcripts, + reindexed to the Cell_Stats cell list so metadata and counts share the exact same cells/order. +- No `-polygons.csv`: the raw export has no cell boundaries → `cells_polygons=False`, + `cell_boundaries` absent (optional in the API; downstream methods segment themselves). + +**Coordinate formula** (the load-bearing part — validated end-to-end at **0.998** of assigned +transcripts landing on their own cell in sopa's stitched labels, on a 3-FOV subset run through +the real `sopa.io.cosmx`): +``` +scale = 1e3 / COSMX_PIXEL_SIZE # 0.120280945 µm/px; == sopa's own mm→px factor +x_global_px = x_mm * scale + x_local +y_global_px = y_mm * scale + (H - 1 - y_local) # H=4256; sopa flips each label tile along y +``` +The transcript local `(x, y)` map DIRECTLY to `CellLabels[row=y, col=x]` (top-origin, no flip; +measured 0.993 raw), but sopa's stitcher does `da.flip(tile, axis=y)`, so the `(H-1-y)` term +undoes that flip. The same transform is applied to the Cell_Stats centroids so the table's +`obsm['spatial']` stays consistent with transcripts + labels. `fov_shift` is left at sopa's +inferred default (False — no polygons to infer from). + +**Global-cell-id assert workaround:** sopa derives `max(cell_ID)` separately from each of tx / +metadata / exprMat and asserts they agree; a segmented cell with zero transcripts makes the +tx-file max smaller and trips the assert. In reconstruction mode we monkeypatch +`_CosMXReader._get_global_cell_id` to use one fixed max (the segmentation's true max from +Cell_Stats), keeping ids consistent across all three files. + +**Not yet run end-to-end on the full 304-FOV liver dataset on Nebius** — validated on a 3-FOV +subset. Watch memory (≈229 M-row tx read, mitigated by the lean-read shim) and wall-time +(reconstruction + stitch + zarr write) against the loader's `hightime` (8 h) cap; bump to +`veryhightime` if it times out. + +## The real OOM: scratch exhaustion during extraction (the actual fix, 2026-07-27) + +**Symptom:** `bruker_mouse_brain_cosmx` (full hemisphere) is OOM-killed (**exit 137**) on +every retry at the *very first* step — the `log("Extract zip of raw files")` line, 58 µs in — +**before any sopa / transcript / image code runs**. + +**Root cause (measured — NOT the transcript read):** the raw `HalfBrain.zip` is **208 GB +uncompressed** (3990 files), but sopa only ever reads **Morphology2D (13.6 GB) + +`FOV*/CellLabels` (0.12 GB)**. ~190 GB of the zip is `Morphology3D` (per-z-plane 3D stacks — +sopa uses the 2D composite), `AnalysisResults` (decoding CSVs), and `RunSummary` — none of +which sopa touches. The loader extracted all of it. And because `--input_raw` was a Viash +`file`, Nextflow first **staged the 187.7 GB zip** into the task's scratch, then the script +extracted ~208 GB *alongside* it. The `veryhighmem` node has only ~220 GB ephemeral disk and a +**RAM-backed scratch** (which is why a *disk* operation OOMs as *memory*): staged zip (187) + +extraction blow the 400→480 GB memory limit → 137. Retries can't help — `disk` is pinned at +200 GB in `labels_nebius.config` and memory only reaches the 480 GB cap. + +**Fix (2026-07-27):** +1. `--input_raw` is now `type: string` in **both** the loader and the `process_bruker_cosmx` + workflow config, so Nextflow no longer stages the 187 GB zip (verified: `main.nf` passes it + through `fromState`). +2. `extract_zip` (`script.py`) streams the zip **in place** — local path via `open()`, + `s3://` via **s3fs** (`_open_zip_source`, 32 MB blocks) — and extracts only members **not** + under `SKIP_DIRS = {Morphology3D, AnalysisResults, RunSummary}` (`_member_wanted`). + Footprint drops **208 GB → 14.6 GB** (916 of 3990 files); S3 transfer drops ~14×. Verified + end-to-end against the real S3 zip (dry-run filter count + streamed a Morphology2D and a + CellLabels tif → valid TIFFs, `strip_root` correct). New dep **s3fs** in `config.vsh.yaml` + pypi ⇒ **the container must be rebuilt** before a run picks this up. + +Denylist (not allowlist) on purpose: it drops only the three known-huge unused dirs and keeps +every small file, so it can't accidentally drop the flat CSVs the **liver** dataset ships +*inside* its raw zip. + +## Secondary optimization: lean transcript read (NOT the OOM — kept anyway) + +The *earlier* hypothesis (recorded here before 2026-07-27) was that the OOM was sopa's +`_CosMXReader.read_transcripts` doing one `pd.read_csv` of the whole `*_tx_file.csv`. That read +is genuinely wasteful (the mouse-brain tx file is 7 GB / ~100M rows, object-dtype strings) and +the lean shim below is worth keeping — but it was **never the OOM**: the job died at +extraction, long before transcripts are read, and the panel is only ~960 genes (so sopa's +table densification is a non-issue too). The shim is only now *reached* for the first time, +once extraction succeeds. + +**Lean transcript read (`script.py`):** we replace `sopa.io.reader.cosmx.pd` with a thin proxy +(`_PandasReadCsvProxy`) that forwards every attribute to real pandas but intercepts +`read_csv`. For the transcripts file only (name contains `_tx_file.csv`), it: +- reads `target` as **categorical** (object → int codes; the single biggest saver), and +- restricts `usecols` to the columns sopa's stitched-FOV path + `_get_global_cell_id` + the + downstream schema actually use: `fov, cell_ID, target, x_global_px, y_global_px, z` + (dropping `cell`, `CellComp`, `x_local_px`, `y_local_px`). + +It pre-reads the header (`nrows=0`) to intersect the keep-list with the real columns and only +leans the read if sopa's required columns are present — otherwise it falls through to a normal +read, so it degrades gracefully on an unexpected export. Every other sopa read (fov positions, +metadata, counts, polygons) is left untouched because their filenames don't match +`_tx_file.csv`. **All of sopa's coordinate/stitching math runs unchanged** on the leaner frame. + +Why the proxy (not patching `pandas.read_csv` globally, and not a private-API reimplement): +sopa is unpinned, so we avoid touching its internals; the proxy only rebinds the module-level +`pd` name *inside sopa's cosmx module*, so no global pandas mutation and no dependence on +sopa's private `read_transcripts` signature. + +Measured ~5.6× smaller (82%) transcript frame on a synthetic CosMx-shaped file; larger on the +real data. Safe to drop the vendor-specific columns because the pipeline is vendor-agnostic — +downstream methods consume the standardized API columns (`x, y, feature_name, cell_id, z`); +nothing depends on CosMx-only `cell`/`CellComp` (verified: the `df["cell"]` refs in +`methods_transcript_assignment/baysor/script_no_sopa.py` are Baysor's *own* output, not our +input). + +## Gotcha: obs `cell_id` vs vendor `cell_ID` (spatialdata write, 2026-07-30) + +`sdata.write()` failed for the mouse brain with `ValidationError: SpatialData contains elements +with invalid names`. spatialdata's writer forbids two keys in a table attribute that differ only +in **case**. sopa's table already carries the CosMx per-FOV-local **`cell_ID`** obs column; the +loader then adds **`cell_id`** (the global unique index — required by `file_common_ist.yaml`'s +`obs.cell_id`). `cell_ID` vs `cell_id` collide → invalid. + +Fix (`script.py`, "Add info to metadata table"): rename the vendor local id to `cell_ID_local` +*before* adding `cell_id`. Note this is a **case-insensitive-key** rule, not a bad-character one +— the message's "invalid names" is a catch-all. (Aside: the mouse panel also has `/` gene names +like `Tuba1a/b/c`; those live in the var **index**, which spatialdata does **not** validate — so +they are harmless for `write` and are left untouched. Only obs/var *column* names, obsm/uns +*keys*, and element names are checked.) + +## Arguments + +| Argument | Required | Notes | +|----------|----------|-------| +| `--input_raw` | yes | **`type: string`** (not a staged file) — URL/path to the raw zip. Streamed in place (s3:// via s3fs) and selectively extracted; see the OOM section. | +| `--input_flat_files` | no | Second zip with the `*_.csv` flat files, when not in the raw zip (mouse brain needs it) | +| `--segmentation_id` | default `["cell"]` | Must be exactly `["cell"]` — asserted at `script.py`; CosMx ships only the cell segmentation | +| `--dataset_*` metadata | mixed | Written into `table.uns` | + +## Setup / Docker + +`config.vsh.yaml`: `openproblems/base_python:1` + `__merge__` of +`/src/base/setup_spatialdata_partial.yaml` + `pypi: [sopa, s3fs]`. **sopa is unpinned** — a +load-bearing risk given the transcript shim reaches into `sopa.io.reader.cosmx` (see risk +points). **s3fs** was added 2026-07-27 for streaming the raw zip from `s3://` — adding it means +the container must be rebuilt/pushed before a run sees the extraction fix. + +## Wiring + +- Registered as a Nextflow module dep of the per-dataset workflow + `src/datasets/workflows/process_bruker_cosmx/` (`main.nf` runs loader → optional + `crop_region` → `setState`). +- Resource label: **`[veryhighmem, midcpu, hightime]`** in `config.vsh.yaml`. +- Run script: `scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh` (launches via + Seqera `tw launch` against `build/main`). +- The workflow's `crop_region` step (gated on `crop_region_min_x`) is the escape hatch to + tile/crop the section if it still won't fit — the mouse-brain params don't set it today. + +## Risk points / gotchas + +- **sopa is unpinned** and the transcript fix depends on `sopa.io.reader.cosmx` existing and + reading the tx file via the module-level `pd.read_csv`. If a future sopa refactors the + reader (renames the module, switches to `dask`/`pyarrow`, changes the tx filename match), + the shim silently stops leaning (falls through to a normal read → OOM returns). Consider + pinning sopa if this bites. +- **Secondary spike, not yet fixed:** sopa's `read_tables` densifies the cell×gene matrix via + `csr_matrix(counts.values)` (`io/reader/cosmx.py`). Negligible for a ~1k panel; large if the + mouse-brain export is **WTx (~19k genes)**. Can't fix without reimplementing sopa's table + read — next lever if OOM persists after the transcript fix. +- **Script-only change:** no image rebuild needed, but it must reach `origin/main` and be + regenerated into `build/main` before a Nebius run picks it up (use the `check-component` + skill to confirm deploy-freshness). +- **Not yet confirmed end-to-end:** the fix passes local syntax + logic tests; it has **not** + yet been validated on a real mouse-brain run on Nebius. diff --git a/src/methods_cell_type_annotation/moscot/NOTES.md b/src/methods_cell_type_annotation/moscot/NOTES.md new file mode 100644 index 000000000..5b3c92caa --- /dev/null +++ b/src/methods_cell_type_annotation/moscot/NOTES.md @@ -0,0 +1,155 @@ +# moscot — cell-type annotation (MOSCOT / optimal transport) + +## What this component is + +Cell-type-annotation stage method (API `src/api/comp_method_cell_type_annotation.yaml`, +subtype `method_cell_type_annotation`). It labels each spatial cell by **optimal-transport +mapping** of an scRNA-seq reference onto the spatial cells and then transferring the +reference's cell-type labels through the learned transport plan. + +- Docs: https://moscot.readthedocs.io — Repo: https://github.com/theislab/moscot +- Paper (config `references.doi`, **verified correct**): Klein, Palla, Lange, Klein et al., + "Mapping cells through time and space with moscot", *Nature* (2025), + DOI `10.1038/s41586-024-08453-2`. This is the moscot method paper — **no citation caveat**. +- **GPU-only** in practice: the engine installs `jax[cuda12]` + `moscot` + `flax` + `diffrax` + on `openproblems/base_pytorch_nvidia:1.1.0`, and the Nextflow directives carry `gpuh100`. + +Under the hood it uses `moscot.problems.space.MappingProblem`, a **fused Gromov-Wasserstein +(FGW)** problem solved with OTT-jax: a quadratic (Gromov-Wasserstein) term matches +intra-domain structure (the SC feature space vs the spatial coordinates) and a linear term +matches shared-gene expression across the two domains; `alpha` interpolates between them. + +## script.py — step by step + +1. **`:13-15` — pop `LD_LIBRARY_PATH` before importing jax.** JAX's `jax[cuda12]` wheel ships + its own CUDA/cuDNN; a system `LD_LIBRARY_PATH` pointing at an older cuDNN otherwise wins and + crashes with "Loaded runtime CuDNN library: 9.1.0 but source was compiled with: 9.8.0". + **Load-bearing** — do not remove. +2. `:18-25` — jax GPU sanity prints (version, backend, devices). +3. `:58-59` — read `input_scrnaseq_reference` and `input_spatial_normalized_counts` (h5ad). +4. **`:62-66` — SMALL-DATA GUARD (the key gotcha).** If `adata_sp.n_obs < 10000` it + **overrides** `par['rank'] = -1` (full rank) and `par['tau'] = 1.0` (balanced OT), + regardless of what was passed in. See "Risk points". +5. `:69-71` — assert a `"normalized"` layer exists in both AnnDatas and `centroid_x`/`centroid_y` + exist in the spatial `obs`. +6. `:74-76` — set `X = layers["normalized"]` for both; build `adata_sp.obsm["spatial"]` from the + centroid columns. +7. `:78` — `sc.pp.pca(adata_sc, n_comps=30)` — the SC quadratic cost is built from a **hardcoded + 30-dim PCA**, not raw genes. +8. `:81-86` — `MappingProblem(adata_sc, adata_sp).prepare(sc_attr={obsm, X_pca}, + xy_callback="local-pca")`: SC intra-domain cost = the 30-dim PCA; spatial intra-domain cost = + a local-PCA callback over the spatial coordinates; the linear (fused) term links shared genes. +9. `:92-99` — `mp.solve(alpha=, epsilon=, tau_a=tau, tau_b=tau, rank=, batch_size=)`. +10. `:102-108` — `mp.annotation_mapping(mapping_mode=, annotation_label=celltype_key, + source="src", forward=False)`; the returned per-cell label is written to + `adata_sp.obs[celltype_key]`. (`forward=False`/`source="src"` are hardcoded.) +11. `:111` — write the annotated spatial AnnData. + +## Arguments + +| Arg | Type | Config default | moscot tool default | Effective on the <10k test data | Maps to | +|-----|------|----------------|---------------------|---------------------------------|---------| +| `--alpha` | double | **0.8** | **0.5** | used as-is | `solve(alpha=)` — FGW quadratic↔linear weight | +| `--epsilon` | double | 0.01 | 0.01 | used as-is | `solve(epsilon=)` — entropic regularization | +| `--tau` | double | **0.3** | **1.0** | **forced to 1.0** (`:65-66`) | `solve(tau_a=tau_b=)` — marginal relaxation / unbalancedness | +| `--rank` | integer | **500** | **-1** | **forced to -1** (`:63-64`) | `solve(rank=)` — low-rank OT (−1 = full rank) | +| `--batch_size` | integer | 1024 | `None` | used (but n_obs10k-cell) dataset**: for a real-data sweep, `tau ∈ [0.1, 0.2, 0.3, 1.0]` (the config + `# TODO` notes it "seems only to work with tau=1 on our data") and `rank ∈ [500, 1000, 2000, + 5000, -1]` (the `# TODO` scales rank with cell count, ~5000 for 300k cells). + +**Tier 2 — label read-out (swept):** +- **`mapping_mode`** *(exposed; config `max`, moscot has no tool default — required arg).* + `max` = label of the single highest-mass source cell; `sum` = aggregate plan mass per cell type + then argmax (uses the full coupling, smoother). One meaningful non-default value: `[sum]`. + +**Tier 3 — high-value knobs NOT exposed today (future work; would need config expose + script +wiring + a container rebuild, so NOT in the submittable sweep):** +- **`n_comps`** — the SC PCA dimensionality is hardcoded at 30 (`:78`); it sets the resolution of + the quadratic cost. Highest-value un-exposed knob. +- **`scale_cost`** *(tool default "mean")* — hardcoded (unset ⇒ tool default); how the cost + matrices are normalized before the solve, affects the effective `epsilon`. +- **`initializer`** — for low-rank solves (`rank>0`) the initializer (e.g. `"rank2"`/`"k-means"`) + materially changes convergence; only relevant once `rank` is swept on a large dataset. +- **`scale_by_marginals`** *(annotation_mapping, tool default True)* and **`forward`/`source`** + (hardcoded `forward=False`, `source="src"`) — the label-projection direction/weighting. +- **`threshold` / `max_iterations`** *(convergence: tool `threshold=1e-3`)* — Sinkhorn/LR-GW + convergence controls (quality↔runtime). + +To sweep any Tier-3 knob: add it to `config.vsh.yaml` `arguments:` (default = the moscot default), +thread it into the `solve`/`prepare`/`annotation_mapping` call in `script.py`, then +`viash ns build` + rebuild the container (see `check-component`). A sweep launched against a stale +`build/main` container silently ignores a freshly-added arg — which is why these stay out of the +submittable-now sweep. + +## Risk points / gotchas + +- **The `<10k` small-data guard silently overrides `tau` and `rank`** (`:62-66`). On the 306-cell + test resource this means those two args do nothing — any sweep over them is a no-op. This is + *intended* behavior (unbalanced/low-rank OT is only for large references), not a bug. +- **GPU-only + not installable/verifiable locally.** moscot/jax[cuda12] is not in the local env; + all tool defaults and signatures in this NOTES were read from moscot's upstream source and docs, + not executed here. The sweep has **not** been run end-to-end yet. +- **A hard argmax label is emitted** — the transport plan / per-cell confidence is not exported; + downstream metrics see labels only. +- **`batch_size` is a pure memory knob** (online GW-cost batching to bound GPU memory on large + references); it does not change the output and is held fixed (`n_obs=306 < 1024` ⇒ no batching + on the test data anyway). diff --git a/src/methods_cell_type_annotation/rctd/NOTES.md b/src/methods_cell_type_annotation/rctd/NOTES.md new file mode 100644 index 000000000..568cbb843 --- /dev/null +++ b/src/methods_cell_type_annotation/rctd/NOTES.md @@ -0,0 +1,142 @@ +# RCTD — cell type annotation (NOTES) + +Authoritative how-it-works / why-the-setup / where-it-breaks reference for the RCTD +component. The running iteration log lives in the memory file +`rctd-split-zero-umi-reference-nan.md` (repo auto-memory) — that pointer + log, this +file the depth. + +## What this component is + +- **Stage / API:** cell type annotation (`src/api/comp_method_cell_type_annotation.yaml`). + Inputs `input_spatial_normalized_counts` (the aggregated cell x gene `.h5ad` with + centroids) and `input_scrnaseq_reference` (labelled scRNA-seq atlas); outputs the same + spatial object with a `cell_type` column in `.obs`/`colData`. +- **Unusual:** this is an **R** component (`script.R`, `openproblems/base_r:1`) that wraps + **spacexr::RCTD** (Robust Cell Type Decomposition). It installs `spacexr` from GitHub + (`dmcable/spacexr`) at build time — an un-pinned `install_github`, so the exact spacexr + version floats with the build date (last observed live: 2.2.1). +- **Links:** docs/repo `https://github.com/dmcable/spacexr`. +- **Publication:** Cable et al., *Robust decomposition of cell type mixtures in spatial + transcriptomics*, Nat Biotechnol 2022 (DOI `10.1038/s41587-021-00830-w`, PMID 33603203). + **The config DOI is correct** (it is the RCTD method paper, not a generic framework + paper) — no citation caveat. The paper's headline contribution is the `doublet_mode` + decomposition (each pixel assigned at most two cell types), and platform-effect / DE-gene + selection is the machinery the exposed thresholds control. + +## script.R — step by step + +1. **L27** read the spatial `.h5ad` as a `SingleCellExperiment`. +2. **L30-35** build a `SpatialRNA` "puck" from `centroid_x/centroid_y` and the raw + `counts` assay (RCTD works on **raw integer counts**, not the normalized layer — despite + the input being named `*_normalized_counts`, the raw `counts` assay is what is read). +3. **L38** read the scRNA-seq reference as a `SingleCellExperiment`. +4. **L45-49 — zero-UMI reference fix (load-bearing).** After `process_dataset` subsets the + reference to the shared spatial panel (~hundreds of genes), some reference cells express + none of those genes (`nUMI == 0`). spacexr's `get_cell_type_info` divides each cell by + its `nUMI` → `NaN`; a single `NaN` column poisons `other_mean` (rowMeans) for **every** + cell type in `get_de_genes` → all logFC `NaN` → 0 DE genes → `create.RCTD` aborts with + "fewer than 10 regression differentially expressed genes". Dropping zero-count cells here + is the fix (see memory `rctd-split-zero-umi-reference-nan.md`). **This — not the + thresholds — was the true root cause of the historical abort.** +5. **L52-53 — CELL_MIN_INSTANCE, applied by hand.** Keep only reference cell types with + >=25 cells. This replicates RCTD's `CELL_MIN_INSTANCE=25` default but is a **hardcoded + pre-filter, not an exposed arg** (so it never reaches `create.RCTD`). +6. **L82-87** sanitize cell-type factor levels containing `/` (spacexr's + `check_cell_types` rejects them); keep a safe→original name map to restore labels later. +7. **L89** `Reference(ref_counts, cell_types, min_UMI = 0)` — `min_UMI` hardcoded to 0 so + the (now zero-UMI-free) cells all pass. +8. **L100-105** `create.RCTD(...)` — the **6 exposed threshold args** are forwarded here: + `gene_cutoff`, `fc_cutoff`, `gene_cutoff_reg`, `fc_cutoff_reg`, `UMI_min`, + `UMI_min_sigma`. `max_cores` = `meta$cpus` (performance only). +9. **L106** `run.RCTD(myRCTD, doublet_mode = "doublet")` — `doublet_mode` is **hardcoded**, + not exposed. +10. **L109-114** take `results_df$first_type`; pixels with `spot_class == "reject"` are + relabelled `None_sp`. +11. **L117-128** write predictions back into `colData(sce)$cell_type` (restoring the + original `/`-containing names) and `write_h5ad` the result. + +## Arguments + +Six exposed args, all `create.RCTD` DE-gene / UMI thresholds. **Every config default is +deliberately relaxed away from the spacexr default** — RCTD's defaults are tuned for +whole-transcriptome references (~20k genes) with high per-cell UMIs, whereas iST uses small +curated panels (~100-500 genes) with compressed fold-changes and low per-cell counts, so the +stock defaults strand too few DE genes / drop too many cells. + +| Arg (config) | Config default | spacexr default | Maps to | Meaning | +|---|---|---|---|---| +| `--gene_cutoff` | 0.0 | 0.000125 | `create.RCTD(gene_cutoff)` | min normalized mean expr for a platform-effect DE gene | +| `--fc_cutoff` | 0.1 | 0.5 | `create.RCTD(fc_cutoff)` | min log-FC for a platform-effect DE gene | +| `--gene_cutoff_reg` | 0.0 | 0.0002 | `create.RCTD(gene_cutoff_reg)` | min normalized mean expr for a regression DE gene | +| `--fc_cutoff_reg` | 0.1 | 0.75 | `create.RCTD(fc_cutoff_reg)` | min log-FC for a regression DE gene | +| `--umi_min` | 20 | 100 | `create.RCTD(UMI_min)` | min total UMI per spatial cell to annotate it | +| `--umi_min_sigma` | 20 | 300 | `create.RCTD(UMI_min_sigma)` | min UMI for cells used to fit platform-effect variance | + +Hardcoded (NOT exposed): `doublet_mode="doublet"` (L106), the `>=25` cells/type filter +(≈`CELL_MIN_INSTANCE`, L52), `Reference(min_UMI=0)` (L89). + +## Setup / Docker + +Base `openproblems/base_r:1`; installs `SingleCellExperiment`, `anndataR`, `rhdf5`, +`devtools` via Bioc, then `devtools::install_github('dmcable/spacexr', build_vignettes=FALSE)`. +The `SingleCellExperiment` reinstall comment (config L63-66) is a workaround for a +`SpatialExperiment`/Seurat install-order bug. **The `install_github` is not commit-pinned**, +so the spacexr version floats — a risk point if upstream changes DE-gene defaults or the +`create.RCTD`/`run.RCTD` signatures. + +## Wiring + +- Registered as `rctd` in the `celltype_annotation_methods` stage of the run_benchmark + workflow. Stage default is `tacco`. +- Nextflow labels (config L75): `hightime, midcpu, highmem`. +- Sweep scripts: `scripts/run_benchmark/param_sweep/rctd_params.yaml` + + `run_test_rctd_nebius.sh` (see below). + +## Risk points / gotchas + +- **Small-panel DE-gene floor.** Even with the zero-UMI fix, raising the fold-change + thresholds toward the spacexr defaults shrinks the DE-gene set; on a very small panel this + can drop back under the 10-gene floor and re-trigger the `create.RCTD` abort. The relaxed + defaults exist precisely to stay clear of that floor. **Consequence for the sweep: the + high-threshold variants (values approaching the spacexr defaults) may legitimately fail on + small-panel datasets** — that failure boundary is part of what the sweep measures. +- Un-pinned spacexr (version floats with build date). +- `doublet_mode` and the `CELL_MIN_INSTANCE`-equivalent filter are hardcoded — see + Optimization / tuning below. + +## Optimization / tuning + +Impact tiers (grounding: Cable et al. 2022 + spacexr `create.RCTD`/`run.RCTD` docs). The +sweep in `param_sweep/rctd_params.yaml` varies only **already-exposed** args (so it is +submittable against the current build/main container without a rebuild). + +- **Tier 0 — input, not a parameter.** The scRNA-seq **reference**: its cell-type + granularity and, above all, its **overlap with the spatial gene panel**. RCTD learns + profiles only on the shared-panel genes, so panel size/quality dominates everything below. + Not a sweep axis. +- **Tier 1 — highest impact on quality (EXPOSED, all six on the sweep).** The DE-gene / + UMI thresholds. `fc_cutoff` / `fc_cutoff_reg` are the sharpest levers (they set which genes + are "informative"); `umi_min` sets which cells get annotated at all vs dropped; + `gene_cutoff` / `gene_cutoff_reg` / `umi_min_sigma` are companion filters. **All six are + set to a non-default (relaxed) value, so per the non-default rule all six are on a sweep + axis, each range straddling the relaxed config default and the spacexr default.** Because + the relaxed default already sits at/near the permissive extreme, the meaningful sweep + direction is *toward* the stricter spacexr default. +- **Tier 1 but NOT exposed → Tier 3 (deferred).** `doublet_mode` (`run.RCTD`, hardcoded + `"doublet"`). This is RCTD's single biggest behavioural lever: `"doublet"` assigns ≤2 + types/pixel (the paper's headline mode, aimed at mixed Slide-seq pixels), `"full"` does + unrestricted deconvolution, `"multi"` is a greedy multi-type extension. For one-cell-per- + object segmented iST, `doublet` + `first_type` is a reasonable default, but `full`/`multi` + could change results materially. **Worth exposing** as `--doublet_mode` (thread into the + `run.RCTD` call). Not in this sweep — exposing it needs `viash ns build` + a container + rebuild, so a stale build/main container would silently ignore it. +- **Tier 3 (deferred) — other un-exposed knobs.** `CELL_MIN_INSTANCE` (currently the + hardcoded `>=25` cells/type pre-filter, L52) — exposing it would let the sweep trade rare- + cell-type coverage against profile stability. `Reference(min_UMI=0)` (L89) is intentional + given the zero-UMI fix and is best left fixed. +- **Performance only (fixed, never swept):** `max_cores` (= `meta$cpus`). + +**To promote any Tier-3 knob into a real sweep:** add the arg to `config.vsh.yaml` +(`type`, `default` = the spacexr default), thread it into the `run.RCTD`/`create.RCTD` call +in `script.R`, `viash config view` to validate, then `viash ns build` + rebuild the +container (see `check-component`) before it can be swept on the cloud. diff --git a/src/methods_cell_type_annotation/singler/NOTES.md b/src/methods_cell_type_annotation/singler/NOTES.md new file mode 100644 index 000000000..23de849d5 --- /dev/null +++ b/src/methods_cell_type_annotation/singler/NOTES.md @@ -0,0 +1,147 @@ +# singler — developer notes + +## What this component is + +Cell-type annotation method (stage API `src/api/comp_method_cell_type_annotation.yaml`, +subtype `method_cell_type_annotation`). It is a **thin Python wrapper around +[singler-py](https://github.com/SingleR-inc/singler-py)** (the Python port of the +Bioconductor `SingleR`), a **reference-correlation labeller**: it builds per-label marker +sets from an annotated scRNA-seq reference, then assigns each spatial cell the label whose +reference profile its expression correlates with best (Spearman, per-label score at a +correlation quantile), optionally refined by fine-tuning. + +- CPU-only (no GPU); resource label `[ midtime, midcpu, midmem ]`. +- Docker: `openproblems/base_python:1` + `pypi: [singler]`, merging + `/src/base/setup_spatialdata_partial.yaml`. Merges the base setup and needs no version + gymnastics — nothing load-bearing to record here. +- Links: docs/repo `https://github.com/SingleR-inc/singler-py`. + +### Citation note +`references.doi: 10.1038/s41590-018-0276-y` is **correct** — Aran et al., *Nat Immunol* +2019, the paper that introduced SingleR. (Unlike some components in this repo, the DOI is +the method-specific paper, not a generic framework paper, so no fix is needed.) + +## script.py — step by step + +1. `sce.read_h5ad(input_spatial_normalized_counts)` and a parallel `ad.read_h5ad` of the + same file (L27-28). The `SingleCellExperiment` is used for the matrix; the `AnnData` is + the object we write labels back onto. +2. `sce_ref = sce.read_h5ad(input_scrnaseq_reference)` (L30). +3. Test matrix = the spatial **`counts`** assay, `mat.sorted_indices()` (L34-35). Reference + matrix = the reference **`normalized`** assay, `sorted_indices()` (L37-38). + `sorted_indices()` is the "magic line" that puts the CSR/CSC into the layout singler + expects. +4. `singler.train_single(ref_data=mat_ref, ref_labels=, + ref_features=..., test_features=...)` (L41-44) — builds the prebuilt reference. +5. `singler.classify_single(mat, ref_prebuilt=built)` (L47) — classifies. +6. `adata_sp.obs["cell_type"] = output["best"]` (L49), then `adata_sp.write(output)` (L53). + +### Two things that bite (read before touching this) + +- **The two config-exposed non-IO args are DEAD.** `--celltype_key` (default `cell_type`) + and `--labels_key` (default `cell_labels`) are declared in the merged config but the + script **never reads `par['celltype_key']` or `par['labels_key']`**. `ref_labels` is + hardcoded to `sce_ref.get_column_data().column("cell_type")` (L42), and the output column + is hardcoded to `"cell_type"` (L49). So changing either arg on the command line (or via a + sweep) currently has **no effect whatsoever**. See "Optimization / tuning" — wiring + `celltype_key` is the single highest-leverage fix and the basis of the prepared sweep. +- **Test matrix is raw `counts`, reference is `normalized`.** SingleR scores with a + rank-based (Spearman) correlation, which is invariant to per-cell monotonic scaling, so a + raw-counts test side is mostly harmless; but marker detection on the reference does expect + log-normalized input, which is why the reference uses `normalized`. Not a tuning axis, + just context. + +## Arguments + +| arg | default | what it maps to upstream | status | +|-----|---------|--------------------------|--------| +| `--input_spatial_normalized_counts` | — | test matrix (`counts` assay) + AnnData to annotate | used | +| `--input_scrnaseq_reference` | — | `train_single(ref_data=…)` (`normalized` assay) | used | +| `--input_transcript_assignments` | — | (optional, unused by script) | ignored | +| `--celltype_key` | `cell_type` | *should* select the reference label column for `ref_labels` | **declared but NOT wired** (hardcoded `"cell_type"`) | +| `--labels_key` | `cell_labels` | (spatial label key) | **declared but NOT wired** (never read) | +| `--output` | — | `adata_sp.write(...)` | used | + +None of singler-py's algorithm knobs (`marker_method`, `num_de`, `quantile`, +`use_fine_tune`, `fine_tune_threshold`, `aggregate`) are exposed — `train_single` / +`classify_single` are called entirely at their upstream defaults. + +## Wiring + +- Registered in `src/workflows/run_benchmark/config.vsh.yaml`: dependency list + (`methods_cell_type_annotation/singler`) and the `--celltype_annotation_methods` default + string `ssam:tacco:moscot:mapmycells:tangram:singler:rctd`; fan-out in + `src/workflows/run_benchmark/main.nf` (~L383). Present on `build/main` (config + `target/`). +- Sweep scripts: `scripts/run_benchmark/param_sweep/singler_params.yaml` + + `run_test_singler_nebius.sh`. + +## Optimization / tuning + +**Grounded in the singler-py source** (`_train_single.py`, `_classify_single.py`, +fetched from `SingleR-inc/singler-py`), whose real defaults are: + +`train_single(..., marker_method="classic", num_de=None, aggregate=False, ...)` and +`classify_single(..., quantile=0.8, use_fine_tune=True, fine_tune_threshold=0.05, ...)`. + +### Tier 0 — the input, not a parameter +The **reference** dominates SingleR accuracy: which scRNA-seq atlas, how well its panel +overlaps the iST gene set, and **which annotation granularity** it is labelled at. The +mouse-brain test reference carries several nested label columns — `cell_type`, +`cell_type_level2`, `cell_type_level3`, `cell_type_level4` (coarse→fine). Selecting among +them is exactly what `--celltype_key` is *meant* to do (see Tier 1). + +### Tier 1 — highest impact on output quality +- **`celltype_key` (annotation granularity)** — ALREADY EXPOSED but **UNWIRED** (see the + "dead args" note). Coarser labels → higher per-class accuracy but less resolution; finer + labels (`level3/4`) → more classes, harder correlation separation on a small iST panel. + This is the one meaningful axis for a reference labeller and the basis of the prepared + sweep — **but it only varies output after this one-line fix + a container rebuild**: + + ```python + # L42, replace: + ref_labels = sce_ref.get_column_data().column("cell_type") + # with: + ref_labels = sce_ref.get_column_data().column(par["celltype_key"]) + ``` + +- **`marker_method`** `{"classic","auc","cohens_d"}`, default `"classic"` — chooses how + per-label markers are ranked. On small, low-count iST panels `cohens_d`/`auc` can pick + more robust discriminative genes than classic pairwise log-FC. **Not exposed.** +- **`use_fine_tune`** (default `True`) / **`fine_tune_threshold`** (default `0.05`) — + fine-tuning re-scores the top candidate labels using only their mutual markers; the + single biggest accuracy lever between similar cell types. Currently left at the + (good) default `True`. **Not exposed** (would need exposing to sweep the threshold or to + measure the cost of turning it off). + +### Tier 2 — quality/speed trade-offs +- **`num_de`** (default: auto for classic, else 10) — markers per pairwise comparison. **Not exposed.** +- **`quantile`** (default `0.8`) — quantile of the correlation distribution used for each + label's score; lower is more robust to a few high-correlation outlier genes. **Not exposed.** +- **`aggregate`** (default `False`) — pseudo-bulk the reference for speed on large atlases, + small accuracy change. **Not exposed.** + +### Tier 3 — not exposed by the component (future work) +`marker_method`, `num_de`, `quantile`, `use_fine_tune`, `fine_tune_threshold`, `aggregate` +are all reachable in singler-py but **not surfaced** in `config.vsh.yaml`, and the script +passes none of them. Exposing any of these means: add the arg to `config.vsh.yaml` +`arguments:`, thread it into the `train_single`/`classify_single` call in `script.py`, then +`viash ns build` + rebuild the container (see the `check-component` skill) — a +**stale `build/main` container silently ignores a newly-added arg**, so a sweep over any of +these is NOT submittable until the rebuild lands. + +### The prepared sweep (and its caveat) +`singler_params.yaml` sweeps **`celltype_key`** over `[cell_type_level2, cell_type_level3, +cell_type_level4]` (default `cell_type` covered by the default variant) → 4 variants total. +`celltype_key` already exists in `build/main`'s config, so the sweep **launches without a +rebuild**. HOWEVER, because of the dead-arg bug above, the deployed `build/main` script +hardcodes the `"cell_type"` column and ignores the value — so **until the one-line wiring +fix is applied and the container is rebuilt, the four variants produce identical output** +(the sweep is a no-op that measures nothing). Apply the fix + rebuild before using this +sweep to draw conclusions. + +### Non-default audit +No exposed argument sits at a non-default (non-performance) value that deviates from an +upstream tool default: `celltype_key`/`labels_key` are data-key strings (no upstream tool +default to deviate from), and every singler-py algorithm knob is left at its library +default. Nothing was promoted onto the sweep axis by the non-default audit; `celltype_key` +is on the axis by design choice (Tier-1 granularity), not because it was mis-defaulted. diff --git a/src/methods_cell_type_annotation/ssam/NOTES.md b/src/methods_cell_type_annotation/ssam/NOTES.md new file mode 100644 index 000000000..776cc98f1 --- /dev/null +++ b/src/methods_cell_type_annotation/ssam/NOTES.md @@ -0,0 +1,108 @@ +# ssam — cell-type annotation (SSAM) + +## What this component is + +Cell-type-annotation stage method (API `src/api/comp_method_cell_type_annotation.yaml`, +subtype `method_cell_type_annotation`). It labels each spatial cell with a cell type by the +**SSAM** approach: signature-based, cell-segmentation-free inference of cell types from an +mRNA-density map, then a per-cell majority vote over the transcripts assigned to that cell. + +- Docs: https://ssam.readthedocs.io — Repo: https://github.com/HiDiHlabs/ssam +- Paper (config `references.doi`, verified correct): Park, Choi, Tiesmeyer et al., + "Cell segmentation-free inference of cell types from in situ transcriptomics data", + *Nat Commun* 12, 3545 (2021), DOI `10.1038/s41467-021-23807-4` (PMID 34112806). + No citation caveat — the config DOI is the SSAM method paper. + +**Important implementation caveat.** Despite the name/DOI, this component does **not** run +the original `ssam` package. It calls `txsim.preprocessing.run_ssam` (txsim from +`theislab/txsim@dev`), which internally uses the **`plankton`/`planktonspace`** +re-implementation: `from plankton.utils import ssam`. The Docker setup installs +`planktonspace` + `matplotlib<3.9` on top of the spatialdata/txsim base. So the algorithm +is plankton's `ssam()`, with defaults that differ from the paper (see below). + +## script.py — step by step + +1. Assert `input_transcript_assignments` and `input_scrnaseq_reference` are provided + (both required for this method). +2. Read `input_spatial_normalized_counts` (h5ad), the `transcripts` element of the + `transcript_assignments.zarr`, and the SC reference; set `adata_sc.X` to its + `layers["normalized"]`. +3. Subset the spatial AnnData to genes shared with the SC reference. +4. Call `tx.preprocessing.run_ssam(adata_sp, transcripts.compute(), adata_sc, + um_p_px=par['um_per_pixel'], cell_id_col='cell_id', gene_col='feature_name', + sc_ct_key=par['celltype_key'])`. +5. Copy `obs['ct_ssam']` -> `obs['cell_type']` (string) and write output. + +**What `run_ssam` actually does** (txsim `preprocessing/_ctannotation.py::run_ssam`): +builds a `plankton.SpatialData` from the transcript gene column and `x*um_p_px`, +`y*um_p_px`; derives mean expression signatures per SC cell type; calls +`ssam(sdata, signatures=..., kernel_bandwidth=4, patch_length=1500, threshold_cor=0.2, +threshold_exp=0.1)` to produce a cell-type map; samples the map at every molecule; then +majority-votes a cell type per `cell_id`. `ct_ssam_cert` records the fraction of a cell's +spots that agree with the winning label. + +Two live `# TODO`s in the script flag a known correctness risk: transcripts are passed in +**physical (µm) space**, not pixel space, so `um_p_px` scaling interacts with the fixed +kernel in a way the author suspects yields poor results ("ssam most likely outputs bad +results since the transcripts are provided in physical space instead of pixel space"). This +is exactly why `um_per_pixel` is the meaningful thing to sweep. + +## Arguments + +| Arg | Type | Config default | txsim/tool default | Maps to | +|-----|------|----------------|--------------------|---------| +| `--um_per_pixel` | double | **0.5** | `um_p_px=0.325` | scales transcript x/y before the density map (only exposed tuning knob) | +| `--celltype_key` | string | `cell_type` (from API) | `sc_ct_key='celltype'` | which SC `obs` column holds cell types (schema arg, not a tuning knob) | +| `--input_*` / `--output` | file | — | — | standard stage I/O | + +Note: the config `um_per_pixel` default (0.5) **deviates** from txsim's `um_p_px` default +(0.325). The `# TODO` on the arg ("Should be able to infer this from transcripts") indicates +it is a placeholder rather than a tuned value. + +## Setup / Docker + +Merges `setup_spatialdata_partial.yaml` + `setup_txsim_partial.yaml` on +`openproblems/base_python:1`, then `pypi: [planktonspace, "matplotlib<3.9"]`. The +`matplotlib<3.9` pin is load-bearing (see commit `80c18d704 "ssam matplotlib"`): plankton's +plotting import breaks against matplotlib >= 3.9. No further Docker drama recorded. + +## Wiring + +- Registered as `celltype_annotation_methods` in the run_benchmark workflow config; the + stage default is `tacco`. `main.nf` fans out one variant per enabled annotation method. +- Param sweep: `scripts/run_benchmark/param_sweep/ssam_params.yaml` + + `run_test_ssam_nebius.sh` (CPU compute env, no `gpu` label; enables `tacco` + `ssam`). + +## Optimization / tuning + +**Tier 0 — input / coordinate space (biggest lever, not a sweepable arg).** Transcripts +arrive in µm (physical) space, contradicting SSAM's pixel-grid assumption. Fixing the +coordinate convention (or deriving the true µm/px from the transcript table, per the arg's +TODO) would likely matter more than any single-knob sweep. + +**Tier 1 — the only exposed knob: `um_per_pixel`.** With the SSAM Gaussian +`kernel_bandwidth` fixed at 4 (µm-equivalent) inside txsim, `um_per_pixel` is the *only* +reachable control over the effective KDE smoothing scale: it rescales the point cloud +relative to the fixed kernel. Smaller -> coarser smoothing; larger -> finer. Swept over +`[0.1, 0.2, 0.325, 1.0]` (straddling the txsim default 0.325 and the shipped 0.5, which the +default variant already covers). + +**Tier 3 — high-value knobs NOT exposed today (future work; NOT in this sweep).** These are +the true SSAM quality dials, but they are **hardcoded inside txsim's `run_ssam`** call to +`ssam(...)`, so the component cannot forward them. Exposing them means editing +`txsim/preprocessing/_ctannotation.py` (add params to the `run_ssam` signature + the `ssam()` +call), then adding matching `config.vsh.yaml` args and rebuilding the image +(`viash ns build` + container rebuild on `build/main` — see the `check-component` skill). +Until then they must NOT go in `ssam_params.yaml`: + +- `kernel_bandwidth` (txsim hardcodes **4**; SSAM's own default is 2.5 µm) — the KDE + bandwidth; the single biggest quality dial. Lower = sharper/noisier, higher = smoother. +- `threshold_cor` (hardcoded **0.2**) — minimum correlation between a pixel's local + expression vector and a signature for a cell type to be called; SSAM typically uses ~0.6. + Directly trades recall vs precision of assignments. +- `threshold_exp` (hardcoded **0.1**) — minimum total expression for a pixel to be + classified at all (a foreground/vector-field threshold). +- `patch_length` (hardcoded **1500**) — tiling size; primarily a memory/speed knob (would + stay fixed even if exposed). +- `no_ct_assigned_value` (`run_ssam` default `'None_sp'`) — label for unassigned cells; + behavioural, not a quality dial. diff --git a/src/methods_cell_type_annotation/tangram/NOTES.md b/src/methods_cell_type_annotation/tangram/NOTES.md new file mode 100644 index 000000000..40b5c6b53 --- /dev/null +++ b/src/methods_cell_type_annotation/tangram/NOTES.md @@ -0,0 +1,190 @@ +# Tangram cell-type annotation — implementation notes + +Reference for anyone (human or Claude) working on this component. Captures how it +works, which of Tangram's many knobs are actually wired up, why the two shipped +defaults were chosen, and where to tune for quality vs. time. + +Links: docs https://tangram-sc.readthedocs.io · repo +https://github.com/broadinstitute/Tangram. + +**Citation:** the config's `references.doi` is `10.1038/s41592-021-01264-7` — +Biancalani et al., "Deep learning and alignment of spatially resolved single-cell +transcriptomes with Tangram", *Nature Methods* 18, 1352–1362 (2021), PMID +34711971. Verified via PubMed: this **is** the method-specific paper (not a generic +framework reference), so no citation caveat is needed. + +## What this component is + +One of the methods at the **cell-type-annotation** stage of the iST preprocessing +benchmark (`src/methods_cell_type_annotation/`). It implements the standard API +`src/api/comp_method_cell_type_annotation.yaml`: + +- in: `input_spatial_normalized_counts` (`.h5ad`, cell × gene with a `normalized` + layer) + `input_scrnaseq_reference` (`.h5ad`, reference atlas with a `normalized` + layer and a `celltype_key` obs column, default `cell_type`). +- out: the spatial AnnData with a per-cell label written to `obs[celltype_key]`. + +It is a **GPU deep-learning method**: Tangram (torch) learns a soft mapping matrix +between the reference and the spatial cells by gradient descent, then projects the +reference labels through that mapping. Nextflow directives are +`[ midtime, midcpu, midmem, gpuh100 ]` — the `gpuh100` label is what schedules it +onto a GPU node. + +## script.py — step by step + +1. **Device pick** (`:19-23`) — `"cuda:0"` if `torch.cuda.is_available()` else + `"cpu"`. Tangram runs on CPU but is much slower there. +2. **Read inputs** (`:30-31`) — spatial + reference AnnData. +3. **Use the log1p-normalized layer as `X`** (`:33-35`) — both adatas' `X` is + replaced by their `normalized` layer before mapping. (That layer is produced + upstream by the normalization stage / the SC process workflow, not here.) + `adata_sp_orig` is stashed so the final output carries the original layers, not + Tangram's intermediate additions. +4. **Markers = the full spatial panel** (`:39-40`) — every `adata_sp.var_name` is + used as a training gene. iST panels are small (hundreds of genes) so "all genes" + is reasonable; there is no marker-selection step. +5. **`tg.pp_adatas`** (`:44`) — intersects genes across the two adatas, drops + all-zero genes, and stores density priors. Mutates both adatas in place. +6. **Build `map_kwargs` and map** (`:56-66`) — calls `tg.map_cells_to_space` with + `mode=par['mode']`, `num_epochs=par['num_epochs']`, `density_prior='uniform'` + (hardcoded), `device`. In `clusters` mode it also passes + `cluster_label=par['celltype_key']` (`:64-65`) — required so Tangram averages + reference cells within each cluster. +7. **Project labels** (`:69-73`) — `tg.project_cell_annotations` writes a + spot × cell-type score frame to `adata_sp.obsm['tangram_ct_pred']`. +8. **Argmax → label** (`:76-80`) — restore `adata_sp_orig`, then set + `obs[celltype_key] = tangram_ct_pred.idxmax(axis=1)` (hard label = highest-scoring + type). The commented block (`:83-87`) would additionally store a per-cell + confidence score and the full score matrix — currently disabled. +9. **Write** (`:90`). + +## Arguments + +Only **two** knobs are exposed; everything else in `map_cells_to_space` is left at +the Tangram default or hardcoded in the script. + +| config arg | default | Tangram default | forwarded to | note | +|---|---|---|---|---| +| `--mode` | `clusters` | `cells` | `map_cells_to_space(mode=…)` | **deviates** from tool default; see below | +| `--num_epochs` | `1000` | `1000` | `map_cells_to_space(num_epochs=…)` | matches tool default | + +**Why `mode: clusters` (a deliberate non-default).** In `cells` mode Tangram learns +an `(n_sc_cells × n_spatial_cells)` mapping matrix on the GPU, which exhausts VRAM +for large references (the config comment cites the 101k-cell MPII lung reference). +`clusters` mode maps per-cell-type *clusters* instead, shrinking the matrix by +~cells-per-type so it fits comfortably. It is also the natural mode for **label +projection** (we only want the cell-type label, not an individual-cell +correspondence). So the deviation is a VRAM/scale choice — but note it *does* change +the output (cluster-level vs. cell-level mapping), so it is treated as a tunable +axis, not a silent baseline (see "Optimization / tuning"). + +**Hardcoded (not a config arg):** `density_prior='uniform'` (`:62`). Tangram's tool +default is `'rna_count_based'`. The docstring guidance (`:48-50`) is: use `'uniform'` +when the spatial voxels are at **single-cell resolution** (MERFISH/Xenium-style), +and `'rna_count_based'` when density is expected to track RNA counts. Since iST here +is single-cell-resolution, `'uniform'` is the appropriate choice — but it is a +deviation from the tool default and is **not exposed**, so it can only be swept +after being added to the config (see Tier 3). + +## config.vsh.yaml — the Docker setup + +- Base image `openproblems/base_python:1` + `pypi: [tangram-sc]`. +- The config carries a **history comment**: `openproblems/base_pytorch_nvidia:1` + (the CUDA/torch base used by other GPU methods) was tried first and "leads to + dependency issues", so it was reverted to the plain python base and lets + `tangram-sc` pull its own torch. A `TODO` notes a different pytorch+CUDA base + might be worth trying. If you touch the image, that dependency clash is the thing + to re-verify. +- `native` engine is also declared (for local/CI without Docker). + +## Optimization / tuning + +Two framing facts before tuning: +1. Unlike the segmentation methods, the tangram defaults are **not** aggressively + speed-floored — `num_epochs=1000` is Tangram's own recommended default. The one + real deviation is `mode` (chosen for VRAM/scale, see above). +2. All ranges below are grounded in the Tangram docs signature (verified against + `tangram/mapping_utils.py`: `map_cells_to_space(mode="cells", + learning_rate=0.1, num_epochs=1000, lambda_g1=1, lambda_d=0, + density_prior='rna_count_based', …)`). + +**Tier 0 — the input, not a parameter.** The biggest lever is the *reference*: which +SC atlas, how well its cell types match the tissue, and the shared gene panel. +`pp_adatas` restricts training to genes present in both — a poorly matched reference +caps achievable accuracy regardless of any knob. Also note markers = the entire +spatial panel (no marker curation). + +**Tier 1 — highest impact on quality (quality vs. time): `num_epochs`** +*(exposed; default 1000 = tool default).* Tangram optimizes the mapping by gradient +descent for `num_epochs` steps; more epochs = better-converged mapping (to a point), +fewer = faster but under-converged and noisier labels. This is the primary +accuracy↔runtime dial. Sweep below and around 1000 (e.g. 100 / 500 / 2000 / 3000) to +find the knee — where extra epochs stop improving the argmax labels. + +**Tier 2 — `mode`** *(exposed; default `clusters`, tool default `cells`).* Not just a +resource knob — it changes the mapping granularity and therefore the labels. +`cells` maps individual reference cells (finer, can be more accurate on **small** +references) but its matrix scales with `n_sc_cells` and OOMs on large ones; the +config picked `clusters` for that reason. The one meaningful non-default value to +sweep is `cells` (feasible on the small test reference; expect OOM on large +references such as MPII lung on limited VRAM). Tangram also offers `constrained`, +but that mode is designed for **cell filtering** and needs `target_count` + +`lambda_count`/`lambda_f_reg` to be meaningful, so it is not a drop-in axis for pure +label projection — left out of the sweep. + +**Tier 3 — high-value knobs NOT exposed today (future work; would need config +expose + container rebuild, so NOT in the submittable sweep):** +- **`density_prior`** — hardcoded `'uniform'`; tool default `'rna_count_based'`. + Directly shapes how mass is distributed across spatial cells. Expose it (choices + `uniform` / `rna_count_based`) to test whether the single-cell-resolution + assumption actually helps on each dataset. Highest-value un-exposed knob. +- **`learning_rate`** *(tool default 0.1)* — step size of the Adam optimizer; + interacts with `num_epochs` (lower LR needs more epochs). Worth a small sweep + (0.05 / 0.1 / 0.5) once epochs are set. +- **`lambda_g1` / `lambda_d` / `lambda_r`** *(1 / 0 / 0)* — the loss weights + (gene-expression term, density term, entropy regularizer). Advanced; only worth + touching after epochs + density_prior are settled. + +To sweep any Tier-3 knob you must add it to `config.vsh.yaml`'s `arguments:` (type + +default = the Tangram default), thread it into the `map_kwargs` dict in `script.py`, +then `viash ns build` + rebuild the container (see the `check-component` skill). A +sweep launched against a stale `build/main` container would silently ignore a +freshly-added arg — which is exactly why these stay out of the current +submittable-now sweep. + +## Risk points / gotchas + +- **`cells` mode OOMs on large references.** The whole reason `clusters` is the + default. Sweeping `mode: cells` is safe on the small test reference but expect + VRAM blow-ups on 100k-cell references at limited VRAM. +- **`density_prior` is hardcoded, not a config arg.** Changing it requires editing + `script.py`, not the sweep file. +- **Only a hard argmax label is emitted** (`:80`); the per-cell confidence score + block is commented out (`:83-87`). Downstream metrics see labels only, no scores. +- **GPU base-image fragility.** The plain `base_python:1` + `tangram-sc` combo is + deliberate; `base_pytorch_nvidia:1` caused dependency issues (see config comment). +- **No end-to-end tuning run has been done yet** — this NOTES.md documents the setup + and the sweep design; results are not yet in. + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: dependency + `methods_cell_type_annotation/tangram` (`:181`); part of the annotation default + string `ssam:tacco:moscot:mapmycells:tangram:singler:rctd` (`:101`). +- `src/workflows/run_benchmark/main.nf:382`: imported in the annotation fan-out. +- Run scripts: enabled (uncommented) in `run_gpu_nebius.sh`; commented in the plain + `run_test_nebius.sh` / `run_test_local.sh` (GPU + time cost). +- **Parameter sweep:** `scripts/run_benchmark/param_sweep/tangram_params.yaml` + (committed default + sweep) and `run_test_tangram_nebius.sh` (Nebius GPU launch, + reads the params file from GitHub raw). The stage default `tacco` is kept alongside + `tangram` so the run has a baseline to compare against. + +## References + +- Docs: https://tangram-sc.readthedocs.io (the + `tangram.mapping_utils.map_cells_to_space` page has the full parameter reference). +- Repo: https://github.com/broadinstitute/Tangram +- **Paper (the DOI in `config.vsh.yaml`, verified correct):** Biancalani T. et al., + "Deep learning and alignment of spatially resolved single-cell transcriptomes with + Tangram", *Nature Methods* 18, 1352–1362 (2021), DOI 10.1038/s41592-021-01264-7, + PMID 34711971. diff --git a/src/methods_expression_correction/resolvi_correction/NOTES.md b/src/methods_expression_correction/resolvi_correction/NOTES.md new file mode 100644 index 000000000..67d24e88c --- /dev/null +++ b/src/methods_expression_correction/resolvi_correction/NOTES.md @@ -0,0 +1,126 @@ +# resolvi_correction — NOTES + +Authoritative how-it-works / why-the-setup / where-it-breaks reference for the +`resolvi_correction` expression-correction method. Committed with the code. + +## What this component is + +- **Stage / API:** `methods_expression_correction` (merges + `/src/api/comp_method_expression_correction.yaml`). Runs late in the iST pipeline: it takes + `spatial_with_cell_types.h5ad` (a cell x gene matrix that already has cell-type labels, + centroids, and a `normalized` layer) and rewrites its `counts` / `normalized` layers with a + denoised, spillover-/background-corrected version. +- **Underlying tool:** **resolVI** (`scvi.external.RESOLVI`) from scvi-tools — a variational + autoencoder that models each observed cell's expression as a mix of (α0) true self-expression, + (α1) diffusion/spillover from spatial neighbours, and (α2) unspecific background, then samples + the posterior to reassign misplaced molecules. +- **GPU:** deep VAE model → tagged `gpu` in the nextflow directives. NB the docker `image:` is + currently `openproblems/base_python:1` (CPU base) with an explicit TODO to move to a GPU + pytorch image (`base_pytorch_nvidia:1`) — see the commented line in `config.vsh.yaml`. So the + component requests a GPU node but the container is a CPU-first python image with `scvi-tools` + pip-installed; whether it actually lights up CUDA depends on the torch build in that image. +- **Links:** docs https://docs.scvi-tools.org/en/latest/user_guide/models/resolvi.html , + repo https://github.com/scverse/scvi-tools . + +## Citation + +`references.doi: 10.1101/2025.01.20.634005` is **correct** (verified via bioRxiv): Ergen & Yosef, +2025, *"ResolVI - addressing noise and bias in spatial transcriptomics"* — the method-specific +paper, not a generic framework citation. No caveat needed. + +## script.py — step by step + +1. Read `input_spatial_with_cell_types`; stash the incoming `normalized` layer as + `normalized_uncorrected` (script.py L32-33). +2. `sc.pp.filter_cells(min_genes=5)` — drop near-empty cells (L36). +3. Build the spatial neighbour input: stack `obs.centroid_x/centroid_y` into + `obsm["X_spatial"]` (L38-39) — this is what resolVI uses to find neighbours. +4. `RESOLVI.setup_anndata(labels_key=celltype_key, layer="counts")` (L45). +5. Construct `RESOLVI(..., semisupervised=True, n_hidden=, encode_covariates=, downsample_counts=)` + (L48-51). `semisupervised=True` is **hard-coded** (tutorials recommend the supervised variant + when cell-type labels are available, which they are at this stage). +6. `train(max_epochs=100)` — **hard-coded** (L52). +7. Two `sample_posterior` calls (**`num_samples=20`, `batch_size=4000` both hard-coded**): + - `model_corrected` → `px_rate` summarised at the q50/median (L55-59). + - `model_residuals` → `mixture_proportions, mean_poisson, per_gene_background, + per_neighbor_diffusion, px_r_inv` (L63-70). The exact `return_sites` set matters — see the + scvi-tools issue #3252 link in the code. +8. Corrected counts = `counts * px_rate_q50 / (1 + px_rate_q50 + mean_poisson)` (L80-81), then a + size-factor renormalisation into the `normalized` layer (L84-86). The renormalisation is + acknowledged as non-ideal in a long code comment (L88-97): the pipeline runs normalization + *before* correction, so the corrected counts are re-normalised here by a crude size factor + rather than re-running the real normalization method. + +**Input-contract note:** the code comment at L28 claims `par['input_sc']` is required, but the +script never reads `input_scrnaseq_reference` — only `input_spatial_with_cell_types`. The SC +reference is optional in the API and unused here. + +## Arguments (exposed in config.vsh.yaml) + +| arg | type | config default | upstream (scvi RESOLVI) default | notes | +|-----|------|----------------|--------------------------------|-------| +| `--celltype_key` | string | `cell_type` | n/a | `labels_key` for `setup_anndata`; input mapping, not a tuning knob (Tier 0) | +| `--n_hidden` | integer | `32` | **32** (matches) | VAE hidden-layer width; model capacity | +| `--encode_covariates` | boolean | `false` | `False` (matches; via `**model_kwargs`) | feed batch/covariates through the encoder | +| `--downsample_counts` | boolean | `true` | `True` (matches) | subsample per-cell counts to equalise depth in training | + +**Hard-coded (NOT exposed) knobs that matter:** `max_epochs=100`, `num_samples=20` (upstream +`sample_posterior` default is **1000**), `mixture_k` (upstream default **50**, unset here), +`batch_size=4000`, `semisupervised=True`, `n_latent=10`, `n_layers=2`, `dropout_rate=0.05`. + +## Optimization / tuning + +Verified upstream defaults (scvi-tools RESOLVI): `n_hidden=32`, `n_hidden_encoder=128`, +`n_latent=10`, `n_layers=2`, `dropout_rate=0.05`, `mixture_k=50`, `downsample_counts=True`; +`sample_posterior(num_samples=1000)`. + +**Non-default audit result:** all three exposed *tunable* args (`n_hidden`, `encode_covariates`, +`downsample_counts`) sit at their scvi upstream defaults — nothing was silently deviated, so +nothing is *forced* onto a sweep axis. They are swept because they are the only exposed levers. + +Tiers: + +- **Tier 0 — input, not a parameter.** `celltype_key` (which obs column carries the labels), + the spatial neighbour graph built from `centroid_x/y`, and the `semisupervised=True` choice. + Biggest lever overall is segmentation/assignment quality upstream, not resolVI's own knobs. +- **Tier 1 — highest posterior-quality impact, but NOT EXPOSED (see below).** `num_samples` + (posterior samples; script uses **20** vs upstream **1000** — the single biggest + quality-vs-time dial for denoising stability) and `max_epochs` (**100**, hard-coded). +- **Tier 2 — exposed quality/capacity trade-offs (what the sweep varies).** + - `n_hidden`: 32 (default) → **64, 128**. More capacity to model diffusion + background, + slower. (script.py's own grid-search note lists exactly 32/64/128.) + - `encode_covariates`: false → **true** (boolean flip). + - `downsample_counts`: true → **false** (boolean flip). +- **Tier 3 — high-value knobs NOT surfaced by the component (future work; would need a config + arg + container rebuild — NOT in the submittable sweep):** + - `num_samples` (posterior samples): expose as e.g. `--num_samples`, sweep `[50, 100, 200]` + straddling the script's 20 toward upstream's 1000. Highest expected quality gain. + - `max_epochs`: expose as `--n_epochs`, sweep around 100 (e.g. `[50, 100, 200, 400]`). + - `mixture_k`: number of background/diffusion mixture components (upstream 50); expose and + sweep e.g. `[10, 50, 100]`. + + Exposing any of these needs: add to `config.vsh.yaml arguments:`, thread `par[...]` into the + matching `RESOLVI(...)` / `train(...)` / `sample_posterior(...)` call in `script.py`, then + `viash ns build` + a container rebuild (see `check-component`) before a build/main Nebius run + would honour them. + +## Wiring + +- Enabled at the `expression_correction` stage. Present (commented) in the stock run scripts + (`run_test_nebius.sh`), active in `run_gpu_nebius.sh` alongside `no_correction`. +- Sweep: `scripts/run_benchmark/param_sweep/resolvi_correction_params.yaml` (committed params) + + `scripts/run_benchmark/param_sweep/run_test_resolvi_correction_nebius.sh` (Nebius launcher, + GPU compute env `5hfmdCBxMRd4nHZaJKYEQZ`, labels include `gpu`). Total variants = + 1 default + 2 (`n_hidden`) + 1 (`encode_covariates`) + 1 (`downsample_counts`) = **5**. + +## Risk points / gotchas + +- **Submittable-now constraint:** the Nebius run uses `--revision build/main`, so the sweep can + only vary args already in the build/main container (the four above). Any newly-exposed Tier-3 + arg is silently ignored until rebuilt. +- **`return_sites` set is load-bearing** (script.py L63-70, scvi issue #3252) — changing it can + break the residual sampling. +- **CPU image / gpu label mismatch** — see "What this component is". Confirm the image actually + provides a CUDA torch before relying on GPU acceleration. +- **Renormalisation is crude** (script.py L84-97) — the corrected `normalized` layer is a + size-factor rescale, not a re-run of the pipeline's real normalization step. diff --git a/src/methods_expression_correction/resolvi_correction/config.vsh.yaml b/src/methods_expression_correction/resolvi_correction/config.vsh.yaml index d494ebcc0..543c379aa 100644 --- a/src/methods_expression_correction/resolvi_correction/config.vsh.yaml +++ b/src/methods_expression_correction/resolvi_correction/config.vsh.yaml @@ -48,7 +48,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0", scvi-tools] + pypi: ["anndata>=0.12.0,<0.13", scvi-tools] - type: native runners: diff --git a/src/methods_expression_correction/split/NOTES.md b/src/methods_expression_correction/split/NOTES.md new file mode 100644 index 000000000..de8090f55 --- /dev/null +++ b/src/methods_expression_correction/split/NOTES.md @@ -0,0 +1,146 @@ +# split — NOTES + +Authoritative how-it-works / why-the-setup / where-it-breaks reference for the `split` +expression-correction method. Committed with the code. + +## What this component is + +- **Stage / API:** `methods_expression_correction` (merges + `/src/api/comp_method_expression_correction.yaml`). Runs late in the iST pipeline: it takes + `spatial_with_cell_types.h5ad` (cell x gene matrix with cell-type labels, centroids, a + `counts` and a `normalized` layer) plus the `scrnaseq_reference.h5ad`, and rewrites the + spatial object's `counts` / `normalized` layers with purified (contamination-corrected) ones. + Unlike resolvi_correction, split **requires the SC reference** (it deconvolves against it). +- **Underlying tools:** **SPLIT** (`bdsc-tds/SPLIT`, R) layered on top of **RCTD** (spacexr). + RCTD deconvolves each spatial cell against the SC reference in `doublet_mode`; SPLIT then + "purifies" the counts — for a doublet it splits the mixed transcripts to the primary cell + type, removing spillover/contamination. R method (`script.R`), CPU-only. +- **GPU:** none. Nextflow directives `[ hightime, highcpu, highmem ]`; no `gpu` label. +- **Links:** docs/repo https://github.com/bdsc-tds/SPLIT . + +## Citation + +`references.doi: 10.1101/2025.04.23.649965` is **correct** (verified via bioRxiv): Bilous, +Buszta, et al., 2025, *"From Transcripts to Cells: Dissecting Sensitivity, Signal +Contamination, and Specificity in Xenium Spatial Transcriptomics"* — the paper that introduces +SPLIT (Spatial Purification of Layered Intracellular Transcripts). Method-specific, not a +generic framework citation. No caveat needed. + +## script.R — step by step + +1. Read `input_spatial_with_cell_types` twice — as `SingleCellExperiment` (`sce`) and as a + `Seurat` object (`xe`) (L30-31). `sce` carries the assays/coords; `xe` supplies the raw + `counts` layer handed to `SPLIT::purify`. +2. If `!keep_all_cells`: drop zero-count cells from both objects (L34-38). +3. Build the `SpatialRNA` puck from `centroid_x/centroid_y` + the `counts` assay (L41-46). +4. Read the SC reference; **drop reference cells with zero panel-gene counts** (L56-60) — after + subsetting to the shared ~few-hundred-gene panel, some cells have `nUMI==0`, which makes + spacexr's per-cell normalisation `NaN` and poisons every cell type's DE-gene means (see the + memory note `rctd-split-zero-umi-reference-nan`). This guard is load-bearing. +5. Keep only cell types with >=25 cells (RCTD minimum) (L62-64); sanitise `/` out of cell-type + names (spacexr's `check_cell_types` rejects them) (L76); build the `Reference(min_UMI=0)`. +6. `create.RCTD(puck, reference, max_cores, gene_cutoff, fc_cutoff, gene_cutoff_reg, + fc_cutoff_reg, UMI_min, UMI_min_sigma)` (L89-94) — **this is where all six exposed threshold + args are forwarded.** Then `run.RCTD(doublet_mode = "doublet")` (L95) — `doublet_mode` is + **hard-coded**. +7. `SPLIT::run_post_process_RCTD(myRCTD)` (L106) then `SPLIT::purify(counts=, + rctd=RCTD, DO_purify_singlets=TRUE)` (L110-114) — `DO_purify_singlets` is **hard-coded TRUE** + and `DO_remove_residual_contamination` is **left at its package default (not set)**. +8. Stash the incoming `normalized` as `normalized_uncorrected` (L121); write purified counts + into `corrected_counts` (copy counts, then overwrite the updated cells) (L124-127); + `logNormCounts` with library-size factors -> new `normalized` (L130-131); `write_h5ad` (L136). + +## Arguments (exposed in config.vsh.yaml) + +All six thresholds map straight into the `create.RCTD(...)` call; `keep_all_cells` is a +component-local zero-cell guard (no upstream equivalent). + +| arg | type | config default | spacexr (RCTD) default | notes | +|-----|------|----------------|------------------------|-------| +| `--keep_all_cells` | boolean | `false` | n/a | drop zero-count cells; description warns TRUE "may cause errors". Robustness flag, not a tuning lever. | +| `--gene_cutoff` | double | `0.0` | **0.000125** | min normalized mean expr, platform-effect DE gene | +| `--fc_cutoff` | double | `0.1` | **0.5** | min log-FC, platform-effect DE gene | +| `--gene_cutoff_reg` | double | `0.0` | **0.0002** | min normalized mean expr, regression DE gene | +| `--fc_cutoff_reg` | double | `0.1` | **0.75** | min log-FC, regression DE gene | +| `--umi_min` | integer | `20` | **100** | min total UMI per spatial cell to include | +| `--umi_min_sigma` | integer | `20` | **300** | min UMI for cells fitting the platform-effect variance | + +**Why the six defaults are relaxed (deliberate, not accidental):** RCTD's DE-gene / UMI +thresholds are tuned for whole-transcriptome references (~20k genes). iST panels are small +(~100-500 genes) and low-count, so the stock spacexr thresholds strand too few DE genes and +`create.RCTD` aborts with *"fewer than 10 regression differentially expressed genes"*. These +defaults mirror the standalone `rctd` cell-type-annotation component exactly. + +**Hard-coded (NOT exposed) SPLIT/RCTD knobs that matter:** `doublet_mode="doublet"` (L95), +`DO_purify_singlets=TRUE` (L113), `DO_remove_residual_contamination` (unset -> package default, +L110-114), `Reference(min_UMI=0)` (L78), the `>=25`-cells-per-type filter (L63). + +## Optimization / tuning + +**Verified upstream defaults** — spacexr `create.RCTD`: `gene_cutoff=0.000125`, `fc_cutoff=0.5`, +`gene_cutoff_reg=0.0002`, `fc_cutoff_reg=0.75`, `UMI_min=100`, `UMI_min_sigma=300`. SPLIT +`purify` (from the repo README): `DO_purify_singlets` and `DO_remove_residual_contamination` +(the latter "removes an additional 2-5% of counts" for improved specificity, v0.3.0+). + +**Non-default audit result:** all six exposed thresholds sit at RELAXED, non-default values +(walked away from the spacexr defaults). None is a performance/resource knob, so per the +phase-3 audit **each is forced onto a sweep axis** — the deviation is itself evidence the knob +matters. `keep_all_cells` is the one exposed arg NOT swept: it has no upstream tool default to +straddle and its own description warns that flipping it to TRUE "may cause errors", so it stays +fixed at `false`. (This is the "do not sweep every knob" — 6 of the 7 exposed args are swept.) + +Tiers: + +- **Tier 0 — input, not a parameter.** The SC reference quality/annotation and the shared panel + gene set (subsetting to a small panel is what forces the relaxed thresholds in the first + place); upstream segmentation/assignment quality. Biggest levers overall. +- **Tier 1 — the exposed thresholds (what the sweep varies).** Each relaxed default is walked + back TOWARD its spacexr default (the meaningful direction; the relaxed default already sits at + the permissive extreme): + - `fc_cutoff`: 0.1 -> **0.25, 0.5** (sharpest — platform-effect DE-gene selection) + - `fc_cutoff_reg`: 0.1 -> **0.4, 0.75** (sharpest — regression DE-gene selection) + - `gene_cutoff`: 0.0 -> **0.0000625, 0.000125** (companion expression-mean filter) + - `gene_cutoff_reg`: 0.0 -> **0.0001, 0.0002** (companion expression-mean filter) + - `umi_min`: 20 -> **50, 100** (which spatial cells get deconvolved) + - `umi_min_sigma`: 20 -> **100, 300** (which cells fit the platform-effect variance) + Total variants = 1 default + 6 axes x 2 values = **13**. +- **Tier 3 — high-value SPLIT-specific knobs NOT surfaced by the component (future work; would + need a config arg + `viash ns build` + container rebuild — NOT in the submittable sweep):** + - `DO_purify_singlets` (SPLIT `purify`, hard-coded TRUE, script.R L113): whether to also + purify singlet cells, not just doublets. Arguably the #1 SPLIT-specific lever. Expose as a + boolean (`--do_purify_singlets`, default true) and sweep the flip to `false`. + - `DO_remove_residual_contamination` (SPLIT `purify`, unset, v0.3.0+): removes an extra 2-5% + of non-reference-supported counts for higher specificity. Expose as a boolean + (`--remove_residual_contamination`, default false) and sweep the flip to `true`. + - `doublet_mode` (RCTD `run.RCTD`, hard-coded "doublet", script.R L95): "doublet" is what + SPLIT's purification consumes; "full"/"multi" change behaviour materially — expose only if + the purification path is adapted to accept them. + + Exposing any of these needs: add to `config.vsh.yaml arguments:`, thread `par[...]` into the + matching `SPLIT::purify(...)` / `run.RCTD(...)` call in `script.R`, then `viash ns build` + a + container rebuild (see `check-component`) before a build/main Nebius run would honour them. + +## Wiring + +- Registered in the run_benchmark workflow: dep list + `src/workflows/run_benchmark/config.vsh.yaml` L186; present in the `--expression_correction_methods` + default string L107 (`no_correction:resolvi_correction:split`). +- Sweep: `scripts/run_benchmark/param_sweep/split_params.yaml` (committed params) + + `scripts/run_benchmark/param_sweep/run_test_split_nebius.sh` (Nebius launcher, standard CPU + compute env `5hfmdCBxMRd4nHZaJKYEQZ`, labels `task_ist_preprocessing,test,split` — no `gpu`). + The launch keeps `no_correction` enabled alongside `split` so the sweep is scored against the + uncorrected baseline; every other stage stays on its single default. + +## Risk points / gotchas + +- **Submittable-now constraint:** the Nebius run uses `--revision build/main`, so the sweep can + only vary args already in the build/main container (the six thresholds). Any newly-exposed + Tier-3 arg is silently ignored until rebuilt. +- **Zero-UMI reference-cell guard (script.R L56-60) is load-bearing** — without it small-panel + references abort RCTD with "fewer than 10 DE genes" (memory: `rctd-split-zero-umi-reference-nan`). +- **High-threshold sweep variants can legitimately fail** — pushing a threshold back to the + spacexr default on a very small panel can drop under the 10-DE-gene floor and abort inside + `create.RCTD`. That failure boundary is part of what the sweep measures, not a bug. +- **SPLIT installed from HEAD** (`remotes::install_github("bdsc-tds/SPLIT")`, unpinned) — the + `purify(counts=, rctd=, ...)` convenience signature the script relies on tracks whatever is on + the default branch; a breaking upstream API change would surface only at container rebuild. diff --git a/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml b/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml index 6f9ea9013..5be364d9c 100644 --- a/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml +++ b/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml @@ -31,7 +31,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0"] + pypi: ["anndata>=0.12.0,<0.13"] - type: native runners: diff --git a/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml b/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml index fca3c2bd6..624672965 100644 --- a/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml +++ b/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml @@ -33,7 +33,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0"] + pypi: ["anndata>=0.12.0,<0.13"] - type: native runners: diff --git a/src/methods_segmentation/binning/NOTES.md b/src/methods_segmentation/binning/NOTES.md new file mode 100644 index 000000000..646fd3538 --- /dev/null +++ b/src/methods_segmentation/binning/NOTES.md @@ -0,0 +1,201 @@ +# Binning segmentation — implementation notes + +Reference for anyone (human or Claude) working on this component. Binning is the +deliberate **"poor segmentation" baseline** of the iST preprocessing benchmark: it +does not look at the image at all — it lays a fixed square grid over the image +extent and calls each tile a "cell." The only knob is the tile (bin) edge length. + +Links: docs/repo https://github.com/openproblems-bio/task_ist_preprocessing. + +**Citation caveat:** the config's `references.doi` is currently +`10.1101/2023.02.13.528102` = Marco Salas et al., *"Optimizing Xenium In Situ data +utility by quality assessment and best practice analysis workflows"* (bioRxiv 2023, +Nilsson lab; theislab co-authors incl. Kuemmerle, Tiesmeyer, Luecken, Theis). That +is the **txsim-associated Xenium best-practices / benchmarking paper** that +benchmarks segmentation tools — an appropriate framework citation for a txsim +segmentation baseline (binning itself has no dedicated publication; it's a trivial +grid). This DOI is **fine as-is**. Note it was *corrected*: the introducing commit +`dd3dff390` had `10.1101/2024.11.07.622434`, which is an **unrelated** single-author +paper ("Integrative cell bin segmentation ... by Voronoi", Ming Lin, Nottingham) — +a wrong DOI, already fixed to the current one. + +## What this component is + +A method at the **segmentation** stage (`src/methods_segmentation/`), implementing +the standard API `src/api/comp_method_segmentation.yaml`: `raw_ist.zarr` in → +`segmentation.zarr` out (a 2D cell **label map** + a stub `table` carrying +`cell_id`/`region`). + +Unlike every other segmentation method here (cellpose, cellposev4, stardist, +watershed), binning **ignores pixel content**. It only uses the image's **shape** +(and, when tuning in microns, its coordinate transform) to tile the field of view +into `bin_size × bin_size` squares. Transcripts are then assigned to whichever tile +they fall in by the downstream transcript-assignment stage, so each square becomes a +"pseudo-cell." This is intentionally a lower bound on segmentation quality — it +captures no cell morphology, splits and merges real cells arbitrarily. + +The heavy lifting is a two-line txsim call, `txsim.preprocessing.segment_binning`. + +## script.py — step by step + +1. **Build hyperparameters** (`:61-65`) — `par.copy()`, map the literal string + `"None"` → `None`, drop `input`/`output`. Leaves `bin_size` (and `bin_size_um`). +2. **Read input** (`:67`) — `sd.read_zarr(par["input"])`. +3. **Pull image + transform** (`:70-71`) — `sdata['image']['scale0'].image` → + `.compute().to_numpy()` for the full-res `(c, y, x)` array, and `.transform.copy()` + for the `scale0` coordinate transform (reattached to the output label map so it + lands in the same coordinate space as the transcripts). +4. **Resolve bin size** (`:73-82`) — if `bin_size_um` is set, convert microns → + pixels via `pixel_size_um(...)` (see below) and `round`; else use `bin_size` + (raw pixels) directly. +5. **Tile** (`:84`) — `tx.preprocessing.segment_binning(image[0], bin_size)`. + `image[0]` is the first channel (a 2D `(y, x)` plane); only its **shape** is + used. The txsim function is literally: + ```python + x = np.floor(np.mgrid[0:n, 0:m][0] / bin_size) # row // bin_size + y = np.floor(np.mgrid[0:n, 0:m][1] / bin_size) # col // bin_size + bins = x*(np.ceil(m/bin_size)) + y + 1 # unique 1-based label per tile + ``` + → a full-coverage label map: every pixel belongs to exactly one square, labels + are `1..N`, no background (label 0). So there are **no gaps** — every transcript + lands in some pseudo-cell. +6. **Downcast + parse** (`:85-88`) — `convert_to_lower_dtype` shrinks the label + array to the smallest uint dtype that holds `max label`; wrap as an + `xarray.DataArray(dims=('y','x'))`, `Labels2DModel.parse` with the copied + transform, store as `labels['segmentation']`. +7. **Carry the metadata table** (`:89-108`) — copies a stub `table` from the input's + `tables["metadata"]`, resolving the downstream-required `cell_id` in priority + order: `obs["cell_id"]` → the table's `instance_key` column + (`uns["spatialdata_attrs"]["instance_key"]`) → the obs index; optional `region` + carried if present. **This block is duplicated in `cellpose/`, + `cellposev4/` and other segmentation scripts** (the instance_key fallback was + added for the Xenium WTA preview behind the Atera dataset). Keep the copies in + sync. Note the table describes the *reference* cells from metadata, not the + binning tiles — downstream counting/aggregation turns the label map into the + actual cell × gene matrix. +8. **Write** (`:110-113`) — `rmtree` any existing output, then `sd_output.write`. + +### `pixel_size_um` helper (`:28-49`) — added this tuning pass + +Reads microns-per-pixel from the image element's coordinate transform: +`element.transform` (dict of coord_system → transformation) → pick `"global"` (fall +back to the first) → `to_affine_matrix(('y','x'),('y','x'))` → average of the +`|diagonal|` scale magnitudes. On failure it prints a warning and returns **1.0** +(the standardized grid — see below), so a micron sweep still produces correct +results on current data even if the transform can't be read. + +## The bin_size unit gotcha (load-bearing) + +`bin_size` is in **PIXELS of the `image` element**, not microns — the txsim +function divides pixel indices by it. The physical size of a bin therefore depends +on the image's pixel size. + +The saving grace: this benchmark's `raw_ist.zarr` **images are rasterized onto a +standardized ~1 µm/pixel grid** (`target_unit_to_pixels: 1` in the loaders, e.g. +`src/datasets/loaders/vizgen_merscope/script.py:245`). Verified directly on the +mouse-brain Xenium test data: the `image` `scale0 → global` transform is +identity-scale + a translation (scale `[1,1,1]`), i.e. **1 px = 1 µm**, and s0 is +`[1, 1000, 1000]` (a 1 mm² crop). So on the standardized grid **`bin_size` in pixels +≈ bin size in microns**. (Contrast the *raw* Xenium morphology at ~0.2125 µm/px — +that resolution is gone by the time binning sees the rasterized `image`.) + +`--bin_size_um` (added this pass, opt-in) makes this explicit and robust: give a +physical edge length and the script converts to pixels using the actual transform, +so bins are comparable even if a future dataset is gridded at a different resolution. + +## config.vsh.yaml — setup + +- Base image `openproblems/base_python:1`; merges `setup_txsim_partial.yaml` + (`theislab/txsim@dev`, for `segment_binning`) + `setup_spatialdata_partial.yaml` + (`spatialdata>=0.7.3`, `anndata>=0.12.0`, `zarr>=3.0.0`). No GPU, no model + download — nothing exotic; it merges the base setups and works. +- Nextflow directives: `[ midtime, midcpu, midmem ]` (4 h / 15 CPU / 50 GB). This is + generous — the actual work is one `np.mgrid` over the image and a write; the cost + is dominated by reading the image and writing the label zarr. CPU-only. + +## Arguments + +| arg | type | default | maps to | +|---|---|---|---| +| `--bin_size` | integer | `30` | tile edge length in **pixels** → `segment_binning(img, bin_size)` | +| `--bin_size_um` | double | *(unset)* | tile edge length in **microns**; when set, overrides `--bin_size` after µm→px conversion | + +`--bin_size_um` is **newly exposed this pass** — it needs `viash ns build` + a +binning **container rebuild** before a run can use it (see check-component). Default +behaviour (bin_size_um unset → use bin_size pixels) is byte-identical to before. + +## Optimization / tuning + +The only lever is the **bin edge length**. There is genuinely nothing else — the +underlying `segment_binning` takes exactly `(img, bin_size)`. + +**Tier 0 — the input.** Binning uses only the image *extent* and pixel size, not +its content, so image channel/stain choice is irrelevant here (the opposite of the +real segmenters). What matters is the **pixel size** of the grid, which sets what a +given `bin_size` means physically — handled by `--bin_size_um` (see gotcha above). + +**Tier 1 — `bin_size` / `bin_size_um` (the whole knob).** Sets the pseudo-cell +size. Grounded in the ~1 µm/px grid and brain/Xenium cell scale (nuclei ~5–8 µm, +whole cells ~10–15 µm): +- too **small** (≈ 10 µm) → many bins per cell → over-segmentation, transcripts of + one cell split across tiles; +- ≈ **15 µm** → roughly one whole cell per bin → the best a blind grid can do; +- default **30 µm** → 2–3 cells per bin → coarse; +- too **large** (≈ 40–50 µm) → strong under-segmentation, several cells merged. +Because binning is a *baseline*, the point of the sweep is to bracket this from +sub-cellular to super-cellular, not to find a "great" value — it shows how much the +downstream benchmark score degrades as pseudo-cells drift from real cell size. + +**Tier 3 — worth adding for a serious study:** none beyond the micron knob already +added. (If ever needed: non-square / hexagonal binning, or an offset/anchor for the +grid origin, would require changes to `segment_binning` upstream in txsim, not just +this component.) + +**Recommended sweep** (see `scripts/run_benchmark/binning_params.yaml`): sweep +`bin_size_um` over `[10, 15, 20, 40, 50]` around the `30` default → 6 variants (a +"star" around the default, one value at a time). Light, as befits a one-knob +baseline. + +## Risk points / gotchas + +- **`bin_size` is pixels, not microns** — the whole "unit gotcha" section. Safe only + because the benchmark grid is standardized to ~1 µm/px; use `--bin_size_um` to be + explicit / dataset-robust. +- **`--bin_size_um` is not yet runtime-validated.** It's the one non-trivial new + code path (reads the transform, converts). The µm→px math and the + `pixel_size_um` transform read need a first real run to confirm (config parses + clean via `viash config view`; the default path is unchanged). The 1.0-µm/px + fallback means a bad transform read degrades to "treat µm as px," which is correct + on current data. +- **Full coverage, no background.** Every pixel gets a positive label — there is no + label 0 / background. Downstream aggregation therefore counts *all* transcripts + into some pseudo-cell; there are no "unassigned" transcripts from the grid itself. +- **Whole image loaded into RAM** (`:70`, `.compute().to_numpy()`), like the other + segmenters — the reason for `midmem` despite the trivial compute. +- **Shared metadata/cell_id block** (`:89-108`) is duplicated across segmentation + scripts; fix together. +- **Only channel 0's shape** is used (`image[0]`); irrelevant here since content is + ignored, but note it if the image is ever multi-channel with differing shapes. + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: dependency + `methods_segmentation/binning` (`:160`), in the default segmentation string + `custom_segmentation:cellpose:cellposev4:binning:stardist:watershed` (`:65`), and + it is the worked example in the `--method_parameters_yaml` doc (`:125`). +- `src/workflows/run_benchmark/main.nf:99`: imported. +- Run scripts: enabled in `run_test_local.sh`, `run_full_local.sh`, + `run_full_nebius.sh`, `run_mpii_nebius.sh`, `run_test_seqeracloud.sh`, + `run_full_seqeracloud.sh`; commented in `run_test_nebius.sh` (the canonical + params-file placement example). Dedicated sweep runners added this pass: + `run_test_binning_local.sh` / `run_test_binning_nebius.sh` + the committed + `binning_params.yaml`. + +## References + +- **Config DOI (appropriate):** Marco Salas S. et al., "Optimizing Xenium In Situ + data utility by quality assessment and best practice analysis workflows", bioRxiv + **10.1101/2023.02.13.528102** (2023) — the txsim-associated Xenium + best-practices/benchmarking paper. +- txsim (`theislab/txsim@dev`) — `segment_binning` in + `txsim/preprocessing/_segmentation.py`. diff --git a/src/methods_segmentation/cellpose/NOTES.md b/src/methods_segmentation/cellpose/NOTES.md new file mode 100644 index 000000000..504b98a7d --- /dev/null +++ b/src/methods_segmentation/cellpose/NOTES.md @@ -0,0 +1,316 @@ +# Cellpose (v3, cyto/nuclei CNN) segmentation — implementation notes + +Reference for anyone (human or Claude) working on this component. Captures how it +works, why it goes through the `txsim` wrapper, why it's a *separate* component +from `cellposev4`, what its defaults really are (they are **not** speed-tuned — +see below), and where the tuning levers are. + +Links: config `links.documentation`/`links.repository` both point at this task's +own repo (not upstream). Upstream Cellpose: docs +https://cellpose.readthedocs.io/en/latest/ · repo +https://github.com/MouseLand/cellpose. + +**Citation caveat (milder than cellposev4's).** The config's `references.doi` is +`10.1038/s41592-020-01018-x` = "Cellpose: a generalist algorithm for cellular +segmentation" (Stringer et al. 2021, Nat. Methods) — the original flow-dynamics +CNN. This component's **default model is `cyto`, which IS that 2021 model**, so the +DOI is appropriate for the shipped default (unlike `cellposev4`, whose DOI pointed +at the wrong paper). Caveat only: the pinned library is Cellpose **3.x**, and the +`--model_type` sweep can select `cyto2` (Cellpose 2.0, Pachitariu & Stringer 2022, +`10.1038/s41592-022-01663-4`) and `cyto3` (Cellpose3, Stringer & Pachitariu 2025, +`10.1038/s41592-025-02595-5`), which come from later papers. If the benchmark +settles on cyto2/cyto3, add those DOIs to the `references` list. + +## What this component is + +One of the methods at the **segmentation** stage of the iST preprocessing +benchmark (`src/methods_segmentation/`). It implements the standard API +`src/api/comp_method_segmentation.yaml`: `raw_ist.zarr` in → `segmentation.zarr` +out (a 2D cell **label map** + a stub `table`). + +It runs **Cellpose v3** (`cellpose<4.0.0`, currently resolves to 3.1.1.x) via the +**txsim wrapper** `txsim.preprocessing.segment_cellpose`, not by calling cellpose +directly. That indirection is the single most important thing to understand about +this component — see the "txsim wrapper" section. The default model is `cyto` +(the original CNN cytoplasm model). + +### Why it exists separately from `cellposev4` + +There are **two** Cellpose components in this repo, and the split is deliberate. +(This table is the mirror of the one in `cellposev4/NOTES.md`.) + +| | `cellpose/` (this one) | `cellposev4/` | +|---|---|---| +| cellpose version | `cellpose<4.0.0` (pinned) | `cellpose>=4.0.0` | +| model | `cyto` (CNN) via `models.Cellpose` | `cpsam` (ViT/SAM) via `CellposeModel` | +| call path | `txsim.preprocessing.segment_cellpose(...)` wrapper | direct `model.eval(...)` | +| extra deps | needs `txsim` | no `txsim` (leaner image) | +| arg surface | ~18 args (channels, model_type, invert, do_3D, …) | 6 args (diameter, flow_threshold, cellprob_threshold, niter, min_size, resample) | +| defaults | == Cellpose library defaults (quality) | deliberately **speed-tuned** | + +History (git `-- src/methods_segmentation/cellpose/`): +- `5487edcfa` "cellpose added" (2024-11) — original component; ~17 args, `base_python`, + txsim wrapper, `models.Cellpose` (v3 API). +- `35061648c` "Support cellpose v4" (2025-06) — reshuffled args to be v4-compatible + (moved `resample`, commented out `interp`/`net_avg`, note "diameter should be None + with v4"). This is why some args are commented out in the config today. +- `8015045ab` "Use cellpose v3 instead of v4 to stay in test time limits" — pinned + `cellpose<4.0.0`, because the v4 SAM/ViT model was **too slow** to finish inside + the CI/`viash test` budget. This component stayed v3; `cellposev4` (#166) was added + as the heavier GPU sibling. +- `7e850cc84` "Add gpu for cellpose" + `702ba3fd2` — added the `gpu` label and moved + to `base_pytorch_nvidia:1`. **But see the GPU gotcha below** — the txsim v3 path + does not actually use the GPU. +- `31d57ec1c` "Add method selection arguments" — merged `setup_spatialdata_partial`. + +## How Cellpose v3 works (background) + +Cellpose predicts, per pixel, a **flow vector** toward the center of the cell that +pixel belongs to, plus a **cell probability**. At inference each pixel "follows the +flow" for `niter` iterations to a fixed point; pixels converging to the same point +become one **mask**. This flow-dynamics core is shared across Cellpose 1–4. + +v3 exposes two model classes (both still present in `cellpose<4`): +- **`models.Cellpose`** — bundles a `SizeModel` (auto-diameter estimator) with a + `CellposeModel`. Constructor default `gpu=False`, `model_type="cyto3"`. This is + what the txsim wrapper uses. +- **`models.CellposeModel`** — the mask model alone (what `cellposev4` uses). + +Builtin models (`MODEL_NAMES`, cellpose 3.1.1.1): `cyto3`, `nuclei`, `cyto2`, +`cyto` (+ specialist `*_cp3`, tissuenet, livecell, bacterial, transformer models). +`cyto`/`nuclei`/`cyto2`/`cyto3` each have a paired size model, so auto-diameter +works for them. `diam_mean` is 30 px for cyto* models, 17 px for `nuclei`. + +## txsim `segment_cellpose` wrapper — the key indirection + +`script.py` calls `tx.preprocessing.segment_cellpose(image[0], hyperparameters)`. +The wrapper lives in `theislab/txsim@dev` +(`txsim/preprocessing/_segmentation.py`, function `segment_cellpose`). For +`cellpose<4` it does, in order: + +1. Reads `hyperparams["model_type"]` (defaults to `'nuclei'` if absent — but this + component always passes `cyto`). +2. `model = models.Cellpose(model_type=model_type)` — **note `gpu=` is NOT passed**, + so it defaults to `gpu=False`. See the GPU gotcha. +3. `del hyperparams["model_type"]` — the model_type key is consumed here and does + **not** reach `eval`. +4. `res, _, _, _ = model.eval(img, channels=[0, 0], **hyperparams)` — + **`channels=[0,0]` is hardcoded** (segment `img` as a single grayscale plane, no + separate nuclear channel). Every remaining config arg is forwarded as a kwarg. + +So the contract is: **all config args except `model_type` are forwarded verbatim +into `Cellpose.eval`** (which passes the non-explicit ones through `**kwargs` into +`CellposeModel.eval`). Adding a new config arg whose name matches an `eval` kwarg +is therefore enough to wire it — no script change needed (that is exactly how the +newly added `--niter` works). + +## script.py — step by step + +1. **Build hyperparameters** (`:34-38`) — `hyperparameters = par.copy()`, then + `{k:(v if v != "None" else None)}` converts the string `"None"` (how the + string-typed args `channel_axis`, `z_axis`, `rescale`, `anisotropy` carry "unset") + back to Python `None`, then `del`s `input`/`output`. Everything else — all the + segmentation knobs — stays in the dict. +2. **Read input** (`:40-43`) — `sd.read_zarr(par["input"])`; + `sdata['image']['scale0'].image.compute().to_numpy()` pulls the full-res image + into RAM; the `scale0` transform is copied to reattach to the output. +3. **Segment** (`:44`) — `tx.preprocessing.segment_cellpose(image[0], hyperparameters)`. + `image[0]` = first channel of a `(c, y, x)` array → a single 2D plane. Returns + the label array. +4. **Downcast + parse** (`:45-48`) — `convert_to_lower_dtype` shrinks the label array + to the smallest uint that holds `max label`; wrap as `xarray.DataArray(dims=('y','x'))`, + `Labels2DModel.parse` with the copied transform, store as + `sd_output.labels['segmentation']`. +5. **Carry the metadata table** (`:50-69`) — copies a stub `table` from the input's + `tables["metadata"]`. Downstream needs a `cell_id` column; resolved in priority + order: explicit `obs["cell_id"]` → the table's `instance_key` column + (`uns["spatialdata_attrs"]["instance_key"]`) → the obs index. Optional `region` + column carried if present. **This block is duplicated verbatim in + `cellposev4/script.py`** — the instance_key fallback was added for exports (e.g. + the Xenium WTA preview behind the Atera dataset) that lack an explicit `cell_id`. + Keep the two copies in sync. +6. **Write** (`:71-74`) — `rmtree` any existing output, then `sd_output.write`. + +Note: the output `table` describes the *reference* cells from metadata, not the +newly predicted masks — segmentation output here is the **label image**; the +counting/aggregation stage downstream turns masks into a cell × gene table. + +## config.vsh.yaml — the Docker setup + +- Base image `openproblems/base_pytorch_nvidia:1` (CUDA/torch), GPU-capable — but + the v3 code path does not use the GPU (gotcha below). +- `setup`: + - `pypi: cellpose<4.0.0` — the version pin that makes this "v3". + - `docker.run`: a **build-time model download** — instantiates + `models.Cellpose(gpu=False, model_type=m)` for `m in ['cyto','nuclei','cyto2','cyto3']`, + which caches each mask model **and its size model** into the image. This was + extended (this tuning session) from caching only `cyto` to caching all four, + because the `--model_type` sweep selects among them and several concurrent tasks + fetching weights caused intermittent `HTTP 504` from the weight server. + - Merges `/src/base/setup_txsim_partial.yaml` (txsim + squidpy + rasterio) **and** + `/src/base/setup_spatialdata_partial.yaml`. Unlike `cellposev4`, this one needs + txsim (for the wrapper). +- Nextflow directives: `[ midtime, midcpu, highmem, gpuhighmem ]` — 4 h / 15 CPU / + 100 GB / high-mem GPU. The `gpuhighmem` label schedules onto a GPU node, but see + the GPU gotcha — the GPU sits idle for the v3 path. + +## Arguments + +The forwarded `Cellpose.eval` / `CellposeModel.eval` args, their config defaults, +and how they compare to Cellpose 3.1.1.1's own defaults. **Almost every default +here already equals Cellpose's library default** — the only deviation is +`model_type` (`cyto` vs library `cyto3`). So, unlike `cellposev4`, this component +is **not** shipped with speed-tuned defaults; it is essentially "stock Cellpose." + +| arg | config default | Cellpose default | notes | +|---|---|---|---| +| `--model_type` | `cyto` | `cyto3` | selects the model; consumed by the wrapper, not passed to eval. `cyto`/`nuclei`/`cyto2`/`cyto3`. **deviation from library default** | +| `--diameter` | `30.0` | `30.` (class) | fixed size assumption. `0`/`None` → auto-estimate via SizeModel (slower). biggest lever | +| `--flow_threshold` | `0.4` | `0.4` | flow-error QC; higher = keep more (more permissive), lower = stricter shape filter | +| `--cellprob_threshold` | `0.0` | `0.0` | per-pixel logit gate (~−6…+6); lower = more/dimmer cells (recall), higher = fewer (precision) | +| `--min_size` | `15` | `15` | drop ROIs below N px; `0` keeps specks | +| `--resample` | `True` | `True` | run dynamics at full resolution (smoother, slower). already the quality setting | +| `--normalize` | `True` | `True` | 1/99-percentile normalization per channel | +| `--niter` | `0` (auto) | `None` (auto) | **newly exposed.** dynamics iterations; 0/None → ∝ diameter (~200 at resample). lower for speed | +| `--augment` | `False` | `False` | test-time augmentation (TTA); quality at ~4–8× cost | +| `--batch_size` | `8` | `8` | tiles per forward pass | +| `--tile_overlap` | `0.1` | `0.1` | tile overlap fraction | +| `--invert` | `False` | `False` | invert intensities before running | +| `--rescale` | `None` | `None` | manual resize factor (only used if diameter is None) | +| `--do_3D` | `False` | `False` | 3D segmentation (input here is 2D) | +| `--anisotropy` | `None` | `None` | 3D only | +| `--stitch_threshold` | `0.0` | `0.0` | 3D stitching of 2D masks | +| `--channel_axis` | `None` | `None` | auto-detected; input is a single 2D plane | +| `--z_axis` | `None` | `None` | 3D only | + +Commented-out in the config: `--net_avg` (removed in cellpose 2.2+), `--tile` +(deprecated), `--interp` (valid in v3 — default `True` — but commented out during +the abandoned v4-compat pass; low sweep value since flipping it off only degrades +2D dynamics). + +## Optimization / tuning + +Two things to understand before tuning (both differ from the `cellposev4` story): + +1. **Defaults are already Cellpose's own defaults, i.e. quality-oriented** (flow QC + on at 0.4, min-mask filtering on at 15, resample on, normalize on). The only + non-stock default is `model_type=cyto`. So "optimize for quality" here is **not** + "walk speed-tuned defaults back to library defaults" (that was the v4 job) — it's + *exploring the model choice, the object scale, and the recall/precision dials*. +2. Grounded in the cellpose 3.1.1.1 source (`models.py`, `dynamics.py`) and docs, + not memory. Ranges tied to Xenium morphology (~0.2125 µm/px). + +**Tier 0 — the input, not a parameter.** The wrapper hardcodes `channels=[0,0]` and +the script feeds only `image[0]` — a **single grayscale plane**. iST morphology is +effectively single-channel here, but if a dataset carried a membrane/boundary stain +alongside the nuclear one, feeding both (via `channels=[cyto_ch, nuc_ch]`) would +improve whole-cell boundaries more than any threshold tweak. Not wired up (would +require changing the hardcoded `channels` in txsim, which is out of this repo). + +**Tier 1 — highest impact on quality** + +- **`model_type`** *(default `cyto`; library default `cyto3`)*. The most + distinctive v3 lever (v4 has no model choice). On a single grayscale morphology + channel the model matters a lot: `nuclei` for a nuclear-dominant stain, `cyto2` + (Cellpose 2.0) and `cyto3` (Cellpose3 generalist super-model) as improved + generalists over the 2021 `cyto`. Sweep `nuclei`/`cyto2`/`cyto3`. +- **`diameter`** *(default 30.0; `0`/`None` = auto)*. Cellpose rescales so objects + land near the model's `diam_mean` (30 px cyto, 17 px nuclei); a fixed 30 assumes + ~30 px objects. Xenium (~0.2125 µm/px): an ~8 µm nucleus ≈ 38 px, a whole cell + ≈ 60–70 px, so 30 can be materially wrong. Sweep `0` (auto), ~40 (nucleus scale), + ~60 (whole-cell scale). +- **`cellprob_threshold`** *(default 0.0, range −6…+6)*. Pure recall↔precision dial, + no speed cost. Lower (−1…−2) recovers more/dimmer cells; raise (+1…+2) suppresses + dim detections. In iST you usually want all cells → this is the most direct dial. +- **`flow_threshold`** *(default 0.4)*. Already at the Cellpose default (QC on). + Raising it (0.6/0.8) keeps cells with higher flow error (more, possibly + ill-shaped); lowering it (0.2) is a stricter shape filter. + +**Tier 2 — quality/speed trade-offs** + +- **`niter`** *(newly exposed; default 0 = auto ≈ 200 at resample)*. Lower (50) is + faster — meaningful because the v3 path is CPU-bound (see gotcha); raise for + large/elongated cells that under-converge. +- **`resample`** *(default True = quality setting)*. Only non-default is `False` + (dynamics on the downsampled grid → faster, coarser boundaries). Include to + characterize the speed cost of the current default. +- **`min_size`** *(default 15)*. `0` keeps small specks (higher recall on tiny + cells); larger (50+) drops debris (higher precision). +- **`augment`** *(default False)*. TTA — the accuracy ceiling at ~4–8× cost; rarely + worth it in a benchmark, include one `True` to bound the gain. +- **`normalize`** *(default True)*. Flip to `False` to see intensity sensitivity; + usually worse for uneven-illumination iST images. + +**Tier 3 — not exposed, worth adding for a serious sweep** + +- **`normalize` as a dict** *(percentiles, `tile_norm`, `sharpen`)*. cellpose v3 + accepts a normalization dict (see `models.normalize_default`). For iST images with + uneven illumination, tile-wise normalization / sharpening can matter; the boolean + `--normalize` can't express it. Would need a JSON/dict-typed arg. +- **`bsize`** *(224)* / **`max_size_fraction`** *(0.4)* — throughput / large-mask + culling; low value for a first sweep. + +**Suggested quality-first order:** `model_type` (nuclei/cyto2/cyto3) → `diameter` +(auto or data scale) → `cellprob_threshold` (−2…+2) → `flow_threshold` → +`min_size`. The current sweep (`scripts/run_benchmark/cellpose_params.yaml`) +covers all of these; **20 variants** total (1 default + 19). + +## Risk points / gotchas + +- **GPU is scheduled but NOT used (v3 path).** The txsim wrapper calls + `models.Cellpose(model_type=...)` without `gpu=`, so it defaults to `gpu=False` and + runs on **CPU** — even though the config uses `base_pytorch_nvidia` and the + `gpuhighmem`/`gpu` label puts it on a GPU node. As of `txsim@dev` at time of + writing, this wastes the GPU allocation and makes runs slow. Fixing it properly + means either patching txsim to pass `gpu=core.use_gpu()` (as the v4 branch of the + same wrapper does) or dropping the GPU label. Left for the user — do not silently + flip infra labels. This is also why the `--niter` speed knob matters here. +- **`--model_type` beyond `cyto` triggers weight downloads unless pre-cached.** The + image now pre-caches `cyto`/`nuclei`/`cyto2`/`cyto3` (config `docker.run`), so the + sweep is safe; if you add another model to the sweep, add it to that list too or + concurrent tasks may hit `HTTP 504` from the weight server. +- **Whole image loaded into RAM** (`:42`) — `.compute().to_numpy()` on `scale0` is + the full-res plane; big panels are why the label is `highmem`. +- **Only channel 0 is segmented, as grayscale** (`image[0]` + hardcoded + `channels=[0,0]`). Any additional stain channels are ignored. +- **`model_type` is consumed by the wrapper, not passed to eval.** If you rename or + drop it, the wrapper falls back to `'nuclei'` — a silent behaviour change. +- **String "None" convention.** `channel_axis`/`z_axis`/`rescale`/`anisotropy` are + typed `string` with default `"None"`; the script converts the literal `"None"` → + Python `None`. A new numeric knob that needs an "auto/unset" sentinel should reuse + a real numeric sentinel the tool understands (as `--niter 0` does), not the string + hack, or it will be forwarded as a string and break eval. +- **Shared metadata block** (`:50-69`) is duplicated in `cellposev4/script.py`; fix + both together. +- **Don't confuse the classes.** This component uses `models.Cellpose` (with + SizeModel); `cellposev4` uses `models.CellposeModel`. `Cellpose` still exists in + v3 but was removed in v4. + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: listed as a dependency + (`methods_segmentation/cellpose`, L158) and included in the default segmentation + method string `custom_segmentation:cellpose:cellposev4:binning:stardist:watershed` + (L65). +- `src/workflows/run_benchmark/main.nf:97`: imported alongside `cellposev4`. +- Run scripts: enabled in `run_full_local.sh`, `run_full_nebius.sh`, + `run_full_seqeracloud.sh`, `run_gpu_nebius.sh`, `run_mpii_nebius.sh`, + `run_test_seqeracloud.sh`; commented out in `run_test_local.sh` / + `run_test_nebius.sh`. Dedicated sweep runners added this session: + `scripts/run_benchmark/run_test_cellpose_local.sh` and `_nebius.sh`, driven by the + committed `scripts/run_benchmark/cellpose_params.yaml`. + +## References + +- Docs: https://cellpose.readthedocs.io/en/latest/ (settings + api pages have the + `model.eval` parameter reference). +- Repo: https://github.com/MouseLand/cellpose (defaults verified against tag + `v3.1.1.1`, `cellpose/models.py` + `cellpose/dynamics.py`). +- txsim wrapper: `theislab/txsim@dev`, `txsim/preprocessing/_segmentation.py`, + `segment_cellpose`. +- **Original Cellpose (`cyto`, the config DOI):** Stringer, Wang, Michaelos, + Pachitariu (2021), Nat. Methods, **10.1038/s41592-020-01018-x**. +- **Cellpose 2.0 (`cyto2`):** Pachitariu & Stringer (2022), Nat. Methods, + **10.1038/s41592-022-01663-4**. +- **Cellpose3 (`cyto3`, image restoration):** Stringer & Pachitariu (2025), Nat. + Methods, **10.1038/s41592-025-02595-5**. diff --git a/src/methods_segmentation/cellposev4/NOTES.md b/src/methods_segmentation/cellposev4/NOTES.md new file mode 100644 index 000000000..145be5a16 --- /dev/null +++ b/src/methods_segmentation/cellposev4/NOTES.md @@ -0,0 +1,291 @@ +# Cellpose 4 (Cellpose-SAM) segmentation — implementation notes + +Reference for anyone (human or Claude) working on this component. Captures how it +works, why it's a *separate* component from `cellpose`, what Cellpose 4 actually +is, and where the settings come from. + +Links: docs https://cellpose.readthedocs.io/en/latest/ · repo +https://github.com/MouseLand/cellpose. + +**Citation caveat:** the config's `references.doi` is `10.1038/s41592-020-01018-x`, +which is the **original Cellpose paper** (Stringer et al. 2021, Nat. Methods) — the +flow-dynamics CNN, *not* the model this component runs. The model actually used +(Cellpose-SAM / `cpsam`) is described in a **different** paper (Pachitariu, +Rariden & Stringer, bioRxiv `10.1101/2025.04.28.651001`, 2025). The config DOI is +the general Cellpose reference, not the v4-specific one; treat it accordingly. + +## What this component is + +One of the methods at the **segmentation** stage of the iST preprocessing +benchmark (`src/methods_segmentation/`). It implements the standard API +`src/api/comp_method_segmentation.yaml`: `raw_ist.zarr` in → `segmentation.zarr` +out (a 2D cell **label map** + a stub `table`). + +It runs **Cellpose 4** (a.k.a. **Cellpose-SAM**, the `cpsam` model) directly via +`cellpose.models.CellposeModel` + `model.eval`. Adapted from +`openproblems-bio/task_spatial_segmentation` (the sibling OpenProblems task) — +the two scripts are nearly identical. + +### Why it exists separately from `cellpose` + +There are **two** Cellpose components in this repo, and this split is deliberate: + +| | `cellpose/` | `cellposev4/` (this one) | +|---|---|---| +| cellpose version | `cellpose<4.0.0` (pinned) | `cellpose>=4.0.0` | +| model | `cyto` (CNN) via `Cellpose` | `cpsam` (ViT/SAM) via `CellposeModel` | +| call path | `txsim.preprocessing.segment_cellpose(...)` wrapper | direct `model.eval(...)` | +| extra deps | needs `txsim` | no `txsim` (leaner image) | +| arg surface | ~17 args (channels, model_type, invert, do_3D, …) | 5 args (diameter, flow_threshold, niter, min_size, resample) | + +History (git): +- `35061648c` "Support cellpose v4" — first attempt to make the *original* + `cellpose` component v4-compatible (dropped `interp`, moved `resample`, noted + `diameter` "should be None with v4"). +- `8015045ab` "Use cellpose v3 instead of v4 to stay in test time limits" — + reverted: pinned `cellpose<4.0.0`, because the v4 SAM/ViT model was **too slow** + to finish inside the CI/`viash test` time budget on a lean/CPU runner. +- `041ab3e23` "Cellposev4 (#166)" — instead of forcing one component to straddle + both, added **this** dedicated `cellposev4` component with a speed-tuned default + param set (see Arguments) and registered it in the benchmark. + +So `cellpose` stays as the fast-enough v3 CNN baseline; `cellposev4` is the +heavier, more-generalizing transformer model run on GPU. + +## What Cellpose 4 / Cellpose-SAM is (background) + +Cellpose predicts, per pixel, a **flow vector** pointing toward the center of the +cell that pixel belongs to, plus a **cell probability**. At inference each pixel +"follows the flow" for a number of iterations to a fixed point; pixels converging +to the same point become one **mask**. (This flow-dynamics core is unchanged from +Cellpose 1–3.) + +**What's new in v4 (Cellpose-SAM, `cpsam`)** — from the paper itself (Pachitariu, +Rariden & Stringer, "Cellpose-SAM: superhuman generalization for cellular +segmentation", bioRxiv `10.1101/2025.04.28.651001`, v1 2025-05-01, corresponding +author C. Stringer, HHMI Janelia; CC-BY-NC). Verified against the bioRxiv abstract: +- They **adapt the pretrained transformer backbone of a foundation model (SAM) + into the Cellpose framework** — i.e. keep Cellpose's flow-dynamics + mask-reconstruction, swap the U-Net CNN backbone for SAM's ViT. The flow → + fixed-point → mask machinery is unchanged from Cellpose 1–3. +- **Headline claim:** the model "substantially outperforms inter-human agreement + and approaches the human-consensus bound" — the paper frames prior methods as + matching inter-human agreement, and argues a human *consensus* could roughly + halve error rates; Cellpose-SAM pushes toward that. Hence "superhuman." +- Generalization robustness was **explicitly trained in**, to: **channel + shuffling** (order-invariance), **cell size**, **shot noise**, **downsampling**, + and **isotropic + anisotropic blur**. This is why in practice you don't need to + set `diameter` or pick channels. +- Positioned as a **foundation model** that drops into the existing Cellpose + ecosystem: finetuning, human-in-the-loop training, image restoration, and 3D + segmentation. + +Consequences in the library API (from the docs / release notes, consistent with +the paper's channel- and size-invariance): +- **Trained around a mean object diameter of 30 px** (range ~7.5–120 px) and + largely **size-invariant** → `diameter` is optional. +- **`channels` is gone** — uses the first ≤3 channels; no channel/model-type + selection. +- **API consolidated:** the old `models.Cellpose` class (which bundled a + `SizeModel` for auto-diameter) and `models.SizeModel` were **removed**; only + `models.CellposeModel` remains, and `CellposeModel()` loads `cpsam` by default. + This is why the script uses `CellposeModel`, not `Cellpose`. +- Newer 4.2+ point releases add `cpsam_v2` / DINOv3 (`cpdino`) variants, but this + component just pins `cellpose>=4.0.0` and takes whatever default ships. + +## script.py — step by step + +1. **Device pick** (`:11-13`) — `torch.device('cuda' if available else 'cpu')`, + printed. Runs on CPU if no GPU, but see the speed caveat below. +2. **Read input** (`:46-49`) — `sd.read_zarr(par["input"])`, then + `sdata['image']['scale0'].image.compute().to_numpy()` pulls the full-res image + into memory, and the `scale0` transform is copied to reattach to the output. +3. **Init model** (`:51-52`) — `CellposeModel(gpu=torch.cuda.is_available())`. + No `model_type`/`pretrained_model` given → loads the default `cpsam` weights. +4. **Build eval params** (`:54`) — collects the 6 tunables + (`diameter, flow_threshold, cellprob_threshold, niter, min_size, resample`) from + `par`, dropping any that are `None`. Note `min_size: -1`, `flow_threshold: 0.0` + and `cellprob_threshold: 0.0` are **passed through** (they're not `None`) — see + Arguments for what those values mean. +5. **Segment** (`:56`) — `model.eval(image[0], progress=True, **eval_params)`. + `image[0]` = first channel of a `(c, y, x)` array → a single 2D plane. Returns + `(masks, flows, styles)`; only `masks` is kept. +6. **Post-process** (`:58-65`) — `convert_to_lower_dtype` downcasts the label + array to the smallest uint that holds `max label` (uint8/16/32/64) to shrink the + zarr; wrap as an `xarray.DataArray(dims=('y','x'))`, `Labels2DModel.parse` with + the copied transform, store as `sd_output.labels['segmentation']`. +7. **Carry the metadata table** (`:67-86`) — copies a stub `table` from the input's + `tables["metadata"]`. Downstream needs a `cell_id` column; it's resolved in + priority order: explicit `obs["cell_id"]` → the table's `instance_key` column + (`uns["spatialdata_attrs"]["instance_key"]`) → the obs index. Optional `region` + column is carried if present. **This exact block is shared verbatim with + `cellpose/script.py`** — the instance_key fallback was added for exports (e.g. + the Xenium WTA preview behind the Atera dataset) that lack an explicit + `cell_id`. Keep the two copies in sync if you touch it. +8. **Write** (`:88-91`) — `rmtree` any existing output, then `sd_output.write`. + +Note: the output `table` describes the *reference* cells from metadata, not the +newly predicted masks — segmentation output here is the **label image**; the +counting/aggregation stage downstream is what turns masks into a cell × gene table. + +## config.vsh.yaml — the Docker setup + +- Base image `openproblems/base_pytorch_nvidia:1` (CUDA/torch) — GPU-capable. +- `setup`: + - `pypi: cellpose>=4.0.0` — the version pin that makes this "v4". + - `script: from cellpose.models import CellposeModel; model = CellposeModel()` — + a **build-time model download**. Instantiating `CellposeModel()` fetches the + `cpsam` weights into the image so runtime tasks never hit the cellpose model + server. (The v3 `cellpose` component does the same thing via a `docker run` + line, for the same reason: several concurrent segmentation tasks all fetching + weights caused intermittent `HTTP 504` from the weight server.) + - Merges `/src/base/setup_spatialdata_partial.yaml` (spatialdata stack). Unlike + `cellpose`, it does **not** merge `setup_txsim_partial.yaml` — no txsim needed. +- Nextflow directives: `[ midtime, midcpu, highmem, gpuhighmem ]` — 4 h / 15 CPU / + 100 GB / high-mem GPU. GPU label is what schedules it onto GPU nodes; on CPU the + SAM model is very slow (the reason v4 was pulled from the CPU test path). + +## Arguments (defaults are speed-tuned) + +The defaults deliberately trade a little accuracy for speed so the transformer +model finishes in budget: + +- `--diameter` (default **30.0**) — expected cell diameter in px. 30 == the model's + training mean, so it's effectively "no rescaling." Left unset, v4 would run a + slower size estimate; fixing it at 30 skips that. Bump it for genuinely large + cells (downsamples → faster + better convergence). +- `--flow_threshold` (default **0.0**) — flow-error QC threshold. Cellpose's own + default is 0.4. **Setting it to 0 disables the flow-consistency check**, which + skips a recompute step → faster, at the cost of possibly keeping some ill-shaped + ROIs. Raise toward 0.4 if you see bad masks. +- `--cellprob_threshold` (default **0.0**) — cell-probability threshold; the + network emits a per-pixel logit in roughly **−6…+6** and only pixels **above** + this value seed masks. This is the Cellpose default (0.0). **Lower it** (e.g. −2) + to recover more / dimmer / low-contrast cells (higher recall); **raise it** + (e.g. +1…+2) to suppress detections in dim regions (higher precision). Unlike the + other knobs it's a pure recall↔precision dial, not a speed dial. Exposed here so + it can be swept per dataset. +- `--niter` (default **10**) — number of flow-dynamics iterations. Cellpose default + is `None` (≈ proportional to diameter, often ~200). **10 is aggressively low** → + fast, but pixels for large/elongated cells may not fully converge. Increase if + cells look under-segmented/split. +- `--min_size` (default **-1**) — minimum pixels per mask. `-1` **disables** + small-mask removal (Cellpose treats `min_size<0` as "off"); skips a filtering + pass. Set a positive value (v3 default was 15) to drop specks. +- `--resample` (default **false**) — whether to run dynamics at the original image + resolution. `false` runs them on the downsampled grid → faster, slightly coarser + boundaries. `true` gives smoother/more precise masks but is slower. + +All six are simply forwarded into `model.eval`. `flow_threshold=0`, +`cellprob_threshold=0` and `min_size=-1` are meaningful values (not "unset") and +*are* passed through — only `None` values get dropped in the eval-params +comprehension (`:54`). + +## Optimization / tuning + +Two things to understand before tuning: +1. **Every shipped default is biased toward speed**, not accuracy — the component + exists precisely because the v4 model is slow (see history). So "optimizing for + quality" here mostly means *walking those defaults back* toward Cellpose's own + defaults, spending the time budget you can afford on GPU. +2. Grounded in the Cellpose-SAM paper + docs (defaults quoted from the docs). The + levers below are ranked by expected impact on iST segmentation quality. + +**Tier 1 — highest impact on quality** + +- **`diameter`** *(exposed; default 30, Cellpose default = auto)*. The biggest + lever. Cellpose rescales the image so objects land near the ~30 px training mean; + `30` therefore means "assume objects are already ~30 px / no rescaling." For + Xenium morphology (~0.2125 µm/px) an ~8 µm nucleus ≈ 35–40 px and a whole cell + ≈ 60–70 px, so a fixed 30 can be materially wrong. v4 is size-*robust* but not + size-*invariant* — tune per dataset, or set `None`/`0` to auto-estimate (slower). +- **`cellprob_threshold`** *(now exposed; default 0.0, range −6…+6)*. Pure + recall↔precision dial. Lower (e.g. −2) → recover more / dimmer / low-contrast + cells; raise (e.g. +1…+2) → suppress dim-region detections. In iST you usually + want to capture *all* cells, so sweeping this (≈ −2 … +2) is the most direct way + to trade missed cells against false ones. No speed cost. +- **`flow_threshold`** *(exposed; default 0.0 here, Cellpose default 0.4)*. `0.0` + **disables** the flow-consistency QC (fast, keeps ill-shaped masks). Restore + ~0.4 to filter malformed ROIs; raise above 0.4 to recover more. + +**Tier 2 — quality/speed trade-offs (already exposed)** + +- **`niter`** *(default 10; Cellpose default None ≈ ∝ diameter, often ~200)*. `10` + is aggressively low → large/elongated cells may not converge → fragmentation. + Raise (or set `None`) for non-round cells. +- **`resample`** *(default False)*. `True` runs dynamics at full resolution → + smoother, more precise boundaries, slower. Turn on when boundary accuracy matters. +- **`min_size`** *(default −1 = off; Cellpose default 15)*. Set positive (≈ 15–50) + to drop debris/specks → higher precision. + +**Tier 3 — not exposed, worth adding for a serious sweep** + +- **`normalize`** *(Cellpose default True, 1/99-percentile)*. Accepts a dict + (percentile bounds, tile-wise normalization, sharpening). For iST images with + uneven illumination / low contrast, the normalization scheme can matter; expose + at least the percentiles if results look intensity-sensitive. +- **`batch_size` / `bsize`** *(tile size, 256 for cpsam)*. Throughput knobs — raise + `batch_size` to trade GPU memory for speed. +- **`augment` / TTA** — test-time augmentation is the accuracy ceiling at ~4–8× + cost; rarely worth it inside a benchmark. + +**Tier 0 — the biggest lever is the input, not a parameter** + +The script segments **`image[0]` — a single channel only** (`:56`). Cellpose-SAM +dropped the `channels` arg *because* it's trained to use up to 3 channels +(cytoplasm + nuclear, any order). If a dataset carries a membrane/boundary stain +alongside the nuclear channel, feeding **both** as a 2-channel stack would improve +*whole-cell* boundaries more than any threshold tweak — the current code discards +that signal. This is the highest-leverage change when input images are +multi-channel (but see the "only channel 0" gotcha below — it's not wired up). + +**Suggested quality-first sweep order:** `diameter` (per dataset or auto) → +`cellprob_threshold` (−2…+2) → `flow_threshold` → 0.4 → `niter` ↑ / `resample` → +True → `min_size` → ~15–50. + +## Risk points / gotchas + +- **Speed / GPU.** cpsam is a ViT — on CPU it can blow the `viash test`/CI time + limit (that's literally why the v3 pin exists on the sibling component). Run this + one on GPU (`gpuhighmem`). If someone flips the test path to include it on a + CPU-only runner, expect timeouts. +- **Whole image loaded into RAM** (`:48`) — `.compute().to_numpy()` on `scale0` is + the full-res plane; big panels are why the label is `highmem`. No tiling/chunking + is done here beyond what cellpose does internally. +- **Only channel 0 is segmented** (`image[0]`). iST morphology images here are + effectively single-channel; if a multi-channel image ever arrives, the other + channels are ignored (v4 could use up to 3 — not wired up). +- **Aggressive defaults can under-segment.** `niter=10`, `flow_threshold=0`, + `min_size=-1`, `resample=false` are all tuned for throughput. If a dataset gives + poor masks, the first knobs to relax are `niter` ↑, `flow_threshold` → 0.4, + `resample` → true. +- **Shared metadata block** (`:67-86`) is duplicated in `cellpose/script.py`; fix + both together. +- **Don't downgrade to the `Cellpose` class.** It no longer exists in v4; only + `CellposeModel` is valid. Likewise there is no `channels`/`model_type` in v4. + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: listed as a dependency + (`methods_segmentation/cellposev4`) and included in the default method string + `custom_segmentation:cellpose:cellposev4:binning:stardist:watershed`. +- `src/workflows/run_benchmark/main.nf:98`: imported alongside `cellpose`. +- Run scripts: enabled in the GPU/full/MPII/seqera run configs + (`run_gpu_nebius.sh`, `run_full_nebius.sh`, `run_mpii_nebius.sh`, + `run_full_seqeracloud.sh`); commented out in `run_test_local.sh` / + `run_test_nebius.sh` (GPU + time cost too high for the quick test path). + +## References + +- Docs: https://cellpose.readthedocs.io/en/latest/ (settings, api pages have the + `model.eval` parameter reference). +- Repo: https://github.com/MouseLand/cellpose +- **Cellpose-SAM (the model this component runs):** Pachitariu M., Rariden M., + Stringer C., "Cellpose-SAM: superhuman generalization for cellular + segmentation", bioRxiv **10.1101/2025.04.28.651001**, v1, 2025-05-01 + (https://www.biorxiv.org/content/10.1101/2025.04.28.651001v1). Not yet a + peer-reviewed journal version at time of writing. +- **Original Cellpose (the DOI actually in `config.vsh.yaml`):** Stringer et al. + (2021), Nat. Methods, DOI 10.1038/s41592-020-01018-x — the flow-dynamics CNN, + i.e. the *framework*, not the v4 model. See the "Citation caveat" up top. +- Adapted from `openproblems-bio/task_spatial_segmentation` (cellpose method). diff --git a/src/methods_segmentation/stardist/NOTES.md b/src/methods_segmentation/stardist/NOTES.md new file mode 100644 index 000000000..0121111b1 --- /dev/null +++ b/src/methods_segmentation/stardist/NOTES.md @@ -0,0 +1,228 @@ +# StarDist2D segmentation — implementation notes + +Reference for anyone (human or Claude) working on this component. Captures how it +works, why the Docker/version pins are the way they are, what the tunable knobs are, +and where it can break. + +Links: docs/repo https://github.com/stardist/stardist · +DOI 10.48550/arXiv.1806.03535. + +**Citation caveat (minor — the DOI is correct):** `references.doi` is +`10.48550/arXiv.1806.03535` = Schmidt, Weigert, Broaddus & Myers, *"Cell Detection +with Star-convex Polygons"*, MICCAI 2018. That **is** the StarDist**2D** method paper, +and this component runs `StarDist2D`, so — unlike `cellposev4` — the DOI matches the +model actually used. Only nuance: the pretrained weights (`2D_versatile_fluo`) were +released later with the library, and there is a separate 3D follow-up (Weigert et al., +*"Star-convex Polyhedra…"*, WACV 2020) that does **not** apply here. No fix needed. + +## What this component is + +One of the methods at the **segmentation** stage of the iST preprocessing benchmark +(`src/methods_segmentation/`). It implements the standard API +`src/api/comp_method_segmentation.yaml`: `raw_ist.zarr` in → `segmentation.zarr` out +(a 2D cell **label map** + a stub `table`). + +It runs **StarDist2D** — a CNN that, for every pixel, predicts a **star-convex +polygon** (a set of radial distances to the object boundary) plus an object +probability, then reconstructs instances by non-maximum suppression over those +polygons. Because objects are represented as star-convex polygons centered on each +nucleus, it excels at **round/blob-like nuclei in crowded fluorescence images** — +the DAPI-like nuclear morphology channel of iST data. It loads a **pretrained** model +(`StarDist2D.from_pretrained`), so there is no training step here. + +Unlike Cellpose-SAM (up to 3 channels), StarDist2D is a **single-channel nuclear +detector**: the script feeds it `image[0]` only (see Tier 0 below). + +## script.py — step by step + +1. **numpy-alias shim** (`:5-15`) — restores the deprecated `np.bool`/`np.int`/… + aliases that numpy≥1.24 removed but stardist/csbdeep still reference. Load-bearing + (see Setup); must run **before** `import stardist`. +2. **Read input** (`:52-54`) — `sd.read_zarr(par["input"])`, then + `sdata['image']['scale0'].image.compute().to_numpy()` pulls the full-res image + into RAM as a `(c, y, x)` array; the `scale0` transform is copied to reattach to + the output. (NB: the *original* 2025-07 version read `morphology_mip`; the current + raw_ist contract exposes the morphology plane as `image` — do not revert.) +3. **Load model** (`:59`) — `StarDist2D.from_pretrained(par['model'])`. Default + `2D_versatile_fluo`. This also loads the model's optimized `thresholds.json` + (prob=0.479071, nms=0.3 for that model) as the fallback thresholds. +4. **Percentile normalizer** (`:64-79`) — a csbdeep `Normalizer` subclass that min-max + scales the image to its **1st / 99.8th percentiles** (`normalize_mi_ma`), the + recommended StarDist preprocessing but with fixed percentile bounds. `block_size` + and `offset` (block overlap/context) are derived from the image width so a large + panel is processed in tiles. +5. **Build eval-params** (`:85-89`) — collects the newly exposed tunables + (`prob_thresh, nms_thresh, scale`) from `par`, **dropping any that are `None`**. + A dropped key ⇒ `predict_instances` uses the model's own optimized value (so the + no-args call is byte-for-byte the pre-tuning behaviour). Mirrors the cellposev4 + eval-params pattern. +6. **Segment** (`:92-95`) — `model.predict_instances_big(image[0], axes='YX', + block_size=…, min_overlap=…, context=…, normalizer=…, **eval_params)`. + `predict_instances_big` splits the image into `block_size` blocks, calls + `predict_instances` on each (forwarding `**eval_params` unchanged — it only + overrides `axes/overlap_label/return_labels/return_predict`), and reassembles the + labels into global coordinates. `image[0]` = first channel → a single 2D plane. +7. **Post-process** (`:100-104`) — `convert_to_lower_dtype` downcasts the label array + to the smallest uint that holds `max label`; wrap as an `xarray.DataArray`, + `Labels2DModel.parse` with the copied transform, store as + `sd_output.labels['segmentation']`. +8. **Carry the metadata table** (`:106-125`) — copies a stub `table` from the input's + `tables["metadata"]`; resolves a required `cell_id` in priority order (explicit + `obs["cell_id"]` → the table's `instance_key` column → the obs index) and carries + an optional `region`. **This block is duplicated verbatim in `cellpose/`, + `cellposev4/` — keep them in sync.** +9. **Write** (`:127-131`) — `rmtree` any existing output, then `sd_output.write`. + +Note: the output `table` describes the *reference* cells from metadata, not the newly +predicted masks — segmentation output here is the **label image**; downstream +aggregation turns masks into a cell × gene table. + +## config.vsh.yaml — the Docker setup (the highest-value section) + +Base image `openproblems/base_tensorflow_nvidia:1` (GPU-capable TF). The pip pins are +**fought-over and load-bearing** — see the history before changing them: + +- **`tensorflow==2.18.0`** + **`tf_keras==2.18.0`** + **`numpy>=2.0.0,<2.2.0`** + + **`scipy<1.15.0`**. The knot: `anndata`/`spatialdata` reference `np.bool` (re-added + in numpy **2.0**), so `numpy<2.0` makes `import anndata` crash; but TF **2.17** caps + `numpy<2.0`. Resolution = bump TF to **2.18** (allows numpy≥2.0), and because + csbdeep/stardist use the **Keras-2 API via `tf_keras`**, `tf_keras` must match the TF + version. `scipy<1.15.0` is pinned for a stardist/csbdeep compatibility break. +- Even with numpy≥2.0, stardist/csbdeep still reference the **removed** scalar aliases + (`np.bool`, `np.long`, …) at import → hence the runtime shim at the top of the script + (`:5-15`). Both halves (the TF/numpy pin **and** the shim) are required together. +- History of this pin (git): `06d998545` add (numpy<2.0, tf loose) → `5058466a7` + add `scipy<1.15.0` → `aa12d73fc` drop the numpy/scipy pins → `15a214eb3` (#167) the + current tf 2.18 / tf_keras 2.18 / numpy 2.0-2.2 / scipy pin + the alias shim. + Don't "simplify" these back. +- Dev note (in-config): on **macOS** the TF install fails; develop in a conda env + (`conda install -c conda-forge tensorflow`), test the Docker via gh-actions. +- Nextflow directives: `[ hightime, midcpu, highmem, gpu ]` — 8 h / 15 CPU / 100 GB / + GPU. TF uses the GPU when present; on CPU it still runs (much lighter than + cellposev4's SAM ViT) but slower. + +## Arguments + +| arg | type | default | maps to | +|-----|------|---------|---------| +| `--model` | string | `2D_versatile_fluo` | `StarDist2D.from_pretrained(model)` | +| `--prob_thresh` | double | *(unset → model 0.479071)* | `predict_instances(prob_thresh=)` | +| `--nms_thresh` | double | *(unset → model 0.3)* | `predict_instances(nms_thresh=)` | +| `--scale` | double | *(unset → no rescale)* | `predict_instances(scale=)` | + +`prob_thresh`, `nms_thresh`, `scale` were **added during tuning** (see below). They are +**optional with no static default on purpose**: leaving them unset makes StarDist use +each model's own optimized `thresholds.json` — hard-coding e.g. `0.479071` would be +wrong for a different `--model` whose optimized thresholds differ. **They need +`viash ns build` + a container rebuild to take effect** (see `check-component`). + +Not exposed: +- `--n_tiles` — a pure GPU-memory tiling knob. `predict_instances_big` **already** tiles + the image via `block_size` (≤4096 px), so `n_tiles` would only sub-tile each block for + the forward pass; it has **no effect on segmentation quality**, only on peak memory. + Deliberately left out of the sweep (would just pad it). +- The normalization percentiles (`1, 99.8`) are hardcoded in the script (`:76`); could + be exposed if results look intensity-sensitive, but not a priority (see Tier 2). + +## Optimization / tuning + +Two things to understand before tuning: +1. Unlike `cellposev4`, the shipped defaults here are **not** aggressively speed-tuned — + `predict_instances` runs with each model's paper-optimized thresholds. So "optimize + for quality" here is mostly **fitting StarDist to iST nuclei** (object size + the + recall/precision thresholds), not walking back throttled knobs. +2. Ranges are grounded in the StarDist paper/docs and the pixel geometry of the data. + The pretrained thresholds quoted (prob 0.479071 / nms 0.3) are those of + `2D_versatile_fluo`. + +**Tier 0 — the input, not a parameter.** The script segments **`image[0]` — a single +channel only** (`:93`). StarDist2D is a single-channel **nuclear** detector (no +membrane/multi-channel mode like Cellpose), so the lever here is **which stain is fed** +and its **normalization** — feeding the crisp nuclear (DAPI-like) plane and getting the +1/99.8-percentile normalization right matters more than any threshold tweak. Its output +is **nuclei**, not whole cells; whole-cell boundaries come from downstream expansion. + +**Tier 1 — highest impact on quality (all newly exposed)** + +- **`prob_thresh`** *(now exposed; default = model 0.479071, range ~0..1)*. Pure + recall↔precision dial and the most direct knob. **Lower** (e.g. 0.3–0.4) → recover + more / dimmer nuclei (higher recall); **raise** (e.g. 0.6–0.7) → keep only confident + detections (higher precision). In iST you usually want to capture all nuclei, so this + is the first axis to sweep. No speed cost. +- **`scale`** *(now exposed; default = none/1.0)*. StarDist's analogue of Cellpose's + `diameter`: rescales the image so nuclei land near the model's trained size. + `2D_versatile_fluo` was trained on a DSB2018 subset (varied, roughly ~20–40 px nuclei); + Xenium morphology (~0.2125 µm/px) makes an ~8 µm nucleus ≈ **38 px**, so a mild + down/upscale can matter. `>1` upscales (small nuclei appear larger), `<1` downscales. + Sweep both sides (≈ 0.5 … 2.0). +- **`model`** *(already exposed; default `2D_versatile_fluo`)*. The other **fluorescence** + pretrained model is `2D_paper_dsb2018`. `2D_versatile_he` is a **brightfield H&E RGB** + model — wrong for single-channel fluorescence — so it is NOT swept. + +**Tier 1/2 — merge control** + +- **`nms_thresh`** *(now exposed; default = model 0.3, range ~0..1)*. Non-max-suppression + IoU for overlapping polygon candidates. **Lower** → more aggressive suppression + (touching nuclei more likely merged/dropped); **raise** → keep more overlapping + detections in crowded tissue. Secondary to `prob_thresh`/`scale`; sweep a few values + (≈ 0.2 … 0.5). + +**Tier 3 — not exposed, worth adding only for a deeper sweep** + +- **Normalization percentiles** (`pmin, pmax`, currently `1, 99.8`). For panels with + uneven illumination / low contrast the percentile bounds change what counts as + background; expose `pmin`/`pmax` if results look intensity-sensitive. +- **`n_tiles`** — memory only, not quality (see Arguments). Not worth a sweep axis. + +**Suggested quality-first sweep order:** `prob_thresh` (0.3…0.7) → `scale` +(0.5…2.0, from pixel size) → `nms_thresh` (0.2…0.5) → `model` → `2D_paper_dsb2018`. +That is exactly the sweep encoded in `scripts/run_benchmark/stardist_params.yaml` +(**13 variants** = 1 default + 4 + 4 + 3 + 1). + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: listed as a dependency + (`methods_segmentation/stardist`, `:161`) and included in the default method string + `custom_segmentation:cellpose:cellposev4:binning:stardist:watershed` (`:65`). +- `src/workflows/run_benchmark/main.nf:100`: imported alongside the other seg methods. +- Sweep run scripts (this tuning): `scripts/run_benchmark/run_test_stardist_local.sh` + and `run_test_stardist_nebius.sh`, sharing the committed + `scripts/run_benchmark/stardist_params.yaml` (local reads it directly; Nebius reads + it from the GitHub raw URL — so it must be **pushed** before launching). +- Already enabled in the standing run configs: `run_gpu_nebius.sh`, `run_full_nebius.sh`, + `run_mpii_nebius.sh`, `run_test_seqeracloud.sh`, `run_full_seqeracloud.sh`; commented + out in `run_test_local.sh` / `run_test_nebius.sh` (GPU/time cost for the quick path). + +## Risk points / gotchas + +- **The TF/tf_keras/numpy/scipy pin quartet + the alias shim are a matched set.** Change + one and imports break (see Setup). This is the single most fragile thing here. +- **New args need a rebuild.** `prob_thresh`/`nms_thresh`/`scale` only exist after + `viash ns build` + a container rebuild; a stale image silently ignores them. The sweep + has **not yet been run end-to-end** with the new args — validated only by + `viash config view` + a script `ast.parse`. +- **`scale` through `predict_instances_big`.** It is forwarded per-block via `**kwargs`; + block geometry (`block_size`/`context`) stays in original-pixel units, so extreme + `scale` values interact with the block overlap. Sanity-check masks at the block seams + if using large `scale`. +- **Only channel 0 is segmented** (`image[0]`). Fine for single-channel iST morphology; + StarDist2D has no multi-channel mode anyway. +- **Whole image loaded into RAM** (`:53`) — full-res plane; big panels are why the label + is `highmem`. Tiling beyond that is `predict_instances_big`'s `block_size`. +- **Shared metadata block** (`:106-125`) is duplicated across the cellpose components; + fix all together. +- **StarDist detects nuclei, not whole cells.** Expect nuclear-sized masks; downstream + volume/expansion stages handle cell extent. + +## References + +- **StarDist2D (the model this component runs, = the config DOI):** Schmidt U., Weigert + M., Broaddus C., Myers G., *"Cell Detection with Star-convex Polygons"*, MICCAI 2018, + arXiv **1806.03535** (DOI 10.48550/arXiv.1806.03535). +- 3D follow-up (not used here): Weigert M., Schmidt U., et al., *"Star-convex Polyhedra + for 3D Object Detection and Segmentation in Microscopy"*, WACV 2020. +- Docs/repo: https://github.com/stardist/stardist (see `examples/other2D/ + predict_big_data.ipynb`, the source of the `predict_instances_big` + custom-normalizer + pattern used here). +- `2D_versatile_fluo` optimized thresholds (`thresholds.json`): prob 0.479071, nms 0.3; + trained on a subset of the DSB 2018 nuclei-segmentation dataset. diff --git a/src/methods_segmentation/watershed/NOTES.md b/src/methods_segmentation/watershed/NOTES.md new file mode 100644 index 000000000..56ceeb697 --- /dev/null +++ b/src/methods_segmentation/watershed/NOTES.md @@ -0,0 +1,238 @@ +# Watershed segmentation — implementation notes + +Reference for anyone (human or Claude) working on this component. Captures how the +classic watershed pipeline works, which of its ~40 exposed knobs actually matter, +and where the sharp edges are. + +Links: docs (this repo) https://github.com/openproblems-bio/task_ist_preprocessing · +implementation https://github.com/theislab/txsim (`txsim.preprocessing.segment_watershed`). + +**Citation caveat:** `references.doi` is `10.1109/34.87344` — Vincent & Soille (1991), +"Watersheds in digital spaces: an efficient algorithm based on immersion simulations", +IEEE TPAMI. That is a legitimate *foundational* watershed reference, so this is **not** a +wrong-paper mismatch like cellposev4. Minor nuance: the code uses +`skimage.segmentation.watershed`, which is the **marker-controlled / priority-flood** +variant (Meyer's flooding, plus Neubert & Protzel for the `compactness` option), not the +Vincent-Soille immersion algorithm per se. Also of historical note: the `repository` link +once (commits `d7c5ae0a0`..`5885bdac8`) pointed at `BennyStrobes/Watershed` — a completely +**unrelated genomics** method — before being corrected to `theislab/txsim` in `5885bdac8`. +The DOI is fine; leave it. + +## What this component is + +One of the methods at the **segmentation** stage of the iST preprocessing benchmark +(`src/methods_segmentation/`). It implements the standard API +`src/api/comp_method_segmentation.yaml`: `raw_ist.zarr` in → `segmentation.zarr` out (a 2D +cell **label map** + a stub `table`). + +Unlike the deep-learning segmenters (cellpose/cellposev4/stardist), this is a **classic, +CPU-only image-processing pipeline** with **no learned model** — a chain of skimage +operations wrapped by `txsim.preprocessing.segment_watershed`. Its distinguishing feature +is an unusually **large, fully function-selectable parameter surface** (~40 args): each +pipeline stage's algorithm is chosen by a `*_func` string, and each stage's parameters are +individually exposed. Most of that surface is low-value plumbing; see "Optimization / +tuning" for the handful of knobs that move the result. + +## How the watershed pipeline works (`segment_watershed`) + +`img[0]` (first channel, a single 2D plane) flows through this fixed chain in txsim +(`txsim/preprocessing/_segmentation.py:382`); each stage is toggled by whether its +`*_func` string maps to a known function: + +1. **Normalize** (`normalize_func` → `adjust_gamma`|`adjust_log`|`adjust_sigmoid`). Default + `gamma` with gamma=1, gain=1 → effectively identity. +2. **Contrast** (`contrast_adjustment_func` → `equalize_adapthist`|`equalize_hist`| + `rescale_intensity`). Default `equalize_adapthist` (CLAHE, clip_limit 0.01). +3. **Blur** (`blur_func` → `gaussian`|`median`). Default `gaussian`, `blur_sigma=1`. +4. **Threshold → binary nuclei mask** (`threshold_func` → `otsu`|`triangle`|`local_otsu`). + Image is first cast to ubyte (`skimage.util.img_as_ubyte`; the module-level `import + skimage` at `_segmentation.py:8` is what makes that name resolve). `local_otsu` + (`rank.otsu`, default) is **adaptive** over a footprint (`threshold_footprint`=square of + size `threshold_footprint_size`=50) → `nuclei = img >= local_threshold`. Global `otsu`/ + `triangle` pick one scalar threshold → `nuclei = img > threshold`. +5. **Mask post-processing** (loop over `post_processing_func_{i}`): default + `remove_small_objects` (min_size 64), then `remove_small_holes` (area_threshold 64), + applied to the **binary mask** before the distance transform. +6. **Distance transform** (`distance_transform_edt`) of the binary mask. +7. **Local-maxima markers** (`find_local_maxima` → `peak_local_max` with + `local_maxima_min_distance`, then `label`). Each maximum becomes one watershed seed → + one cell. +8. **Watershed** — `watershed(-distance, markers, mask=nuclei, watershed_line=True, + **watershed_params)`. Floods basins from the markers within the nuclei mask. +9. **Background-intensity filter** (`filter_cells_based_on_local_background_intensity`, + `_segmentation.py:330`) — if `bg_intensity_filter_bg_factor > 0` (default 0.3), drops + cells whose mean intensity is below `bg_factor ×` local background and reindexes labels. + +## script.py — step by step + +1. **Build hyperparams** (`:34-38`) — `par.copy()`, convert the string `"None"` → Python + `None` (viash passes unset optionals as the literal string `"None"`), drop `input`/ + `output`. Every remaining key is forwarded by name; `segment_watershed` reads what it + recognises via `hyperparams.get(...)`. +2. **Lift compactness into `watershed_params`** (`:40-46`) — see "Optimization"; pops the + exposed `--watershed_compactness` and nests it as `watershed_params={"compactness": …}`, + the only dict `segment_watershed` forwards to `skimage.watershed`. Default 0.0 == + standard watershed (unchanged behaviour). +3. **Read input** (`:48-52`) — `sd.read_zarr`, `sdata['image']['scale0']` pulled full-res + into memory (`.compute().to_numpy()`), and the `scale0` transform copied to reattach. +4. **Segment** (`:53`) — `tx.preprocessing.segment_watershed(image[0], hyperparameters)`. + `image[0]` = first channel only. +5. **Post-process** (`:54-57`) — `convert_to_lower_dtype` downcasts the label array to the + smallest uint that holds `max label`; wrap as `xarray.DataArray(dims=('y','x'))`, + `Labels2DModel.parse` with the copied transform, store as `labels['segmentation']`. +6. **Carry the metadata table** (`:59-78`) — copies a stub `table` from + `tables["metadata"]`; resolves `cell_id` in priority order explicit `obs["cell_id"]` → + `instance_key` column → obs index; carries optional `region`. **This block is shared + near-verbatim with the other segmentation components** (cellpose/cellposev4) — keep in + sync if you touch it. +7. **Write** (`:80-83`) — `rmtree` any existing output, then `sd_output.write`. + +## config.vsh.yaml — Docker setup + +- Base image `openproblems/base_python:1` (plain CPU Python — no CUDA/torch; this method + never uses a GPU). +- `setup` merges `/src/base/setup_txsim_partial.yaml` (needed — the whole method is + `txsim.preprocessing.segment_watershed`) and `/src/base/setup_spatialdata_partial.yaml`. + Commit `050fa15e1` reordered the installs ("Change order of installs in cellpose and + watershed") — the txsim/spatialdata install order is load-bearing for a clean env. +- Nextflow directives: `[ hightime, midcpu, midmem ]` — 8 h / 15 CPU / 50 GB. Note the + **hightime** budget: the classic pipeline (esp. `local_otsu` rank filter with a 50 px + footprint and the O(image) background-intensity filter over large windows) is slow on + big panels even though it is CPU-light per pixel. + +## Arguments (the surface is large; most of it is inert) + +~40 args across 8 stages. Two facts to internalise before reading the list: + +- **`*_func` args are function selectors with tiny, closed mappings.** A value **not** in + the mapping silently sets that stage's function to `None` (skips the stage) — and for + `threshold_func`/`local_maxima_func` that means downstream names (`nuclei`, + `distance_transform`, `local_maxima`) are never defined → the script **crashes** with a + `NameError`. Valid values: `normalize_func` ∈ {gamma, log, sigmoid}; `contrast_ + adjustment_func` ∈ {equalize_adapthist, equalize_hist, rescale_intensity}; `blur_func` ∈ + {gaussian, median}; `threshold_func` ∈ {otsu, triangle, local_otsu}; `local_maxima_func` + ∈ {find_local_maxima}; `post_processing_func_{i}` ∈ {remove_small_objects, + remove_small_holes}. +- **`filter_params` drops any param a chosen function doesn't accept.** So many exposed + args are no-ops for the default function (e.g. the `distance_transform_*` return/indices + args, the `threshold_shift_*`/`hist`/`out`/`mask` args). Don't expect them to do + anything unless you also switch the corresponding `*_func`. + +The knobs that actually move the output (defaults in **bold**): + +| Arg | Default | Effect | +|---|---|---| +| `--threshold_func` | **local_otsu** | foreground/background mask algorithm; see Tier 1 | +| `--threshold_footprint_size` | **50** | neighbourhood (px) for `local_otsu` only | +| `--blur_sigma` | **1** | gaussian pre-smoothing; higher = fewer spurious maxima | +| `--local_maxima_min_distance` | **5** | min px between watershed seeds = merge↔split dial | +| `--post_processing_min_size_1` | **64** | `remove_small_objects` min mask area (px) | +| `--post_processing_area_threshold_2` | **64** | `remove_small_holes` fill area (px) | +| `--bg_intensity_filter_bg_factor` | **0.3** | drop cells dimmer than 0.3× local bg; 0 disables | +| `--watershed_compactness` | **0.0** | NEW; compact-watershed shape regularity (Tier 3) | + +## Optimization / tuning + +Grounded in the txsim implementation and Xenium morphology geometry: at ~0.2125 µm/px an +~8 µm nucleus ≈ **37 px diameter, ~18 px radius, ~1000 px² area** — these anchor the +ranges below. Unlike cellposev4 the defaults here are not obviously *speed*-tuned; they +are a specific (and somewhat arbitrary) pipeline config, and the default `min_distance=5` +in particular tends to **over-split** real nuclei. + +**Tier 0 — the input, not a parameter.** Only `image[0]` (a single channel) is segmented +(`script.py:53`). watershed here targets **nuclei**; whole-cell boundaries are not +recoverable from a nuclear channel alone. If a membrane/boundary stain exists, no threshold +tweak substitutes for it. + +**Tier 1 — highest impact on quality (these are the swept knobs):** + +- **`local_maxima_min_distance`** *(default 5)*. The dominant lever — it sets how close two + watershed seeds may be, i.e. the merge↔split dial. 5 px is far below the nucleus radius + (~18 px), so single nuclei get multiple seeds → over-segmentation. Sweep **[10, 15, 20]** + (walking toward the radius) to merge fragments back into one cell. +- **`threshold_func`** *(default local_otsu)*. Decides *what is foreground* at all. + `local_otsu` is adaptive (good under uneven illumination but slow, footprint 50 px); + global **otsu**/**triangle** use one threshold (fast; `triangle` suits a dominant + background peak, typical of fluorescence). Categorical — sweep the two alternatives. +- **`blur_sigma`** *(default 1)*. Heavier gaussian smoothing removes spurious distance- + transform maxima (fewer over-split cells) but blurs small/dim nuclei. Sweep **[2, 3, 4]** + px (stays sub-nucleus). +- **`post_processing_min_size_1`** *(default 64)*. `remove_small_objects` on the binary mask + — the precision dial against debris/partial nuclei. Sweep **[128, 256, 512]**, all well + below a full nucleus area (~1000 px²) so real nuclei survive. + +**Tier 2 — quality knobs, exposed but not swept:** + +- **`threshold_footprint_size`** *(default 50)* — only active when `threshold_func= + local_otsu`; sets the local-Otsu neighbourhood. Couple it with any `threshold_func` + change. +- **`bg_intensity_filter_bg_factor`** *(default 0.3)* — a real recall↔precision post-filter + (drops dim cells); set 0 to disable, raise to prune harder. Its `window_size`/`bg_size` + (1000/2000) are large and drive much of the runtime. +- **`post_processing_area_threshold_2`** *(default 64)* — `remove_small_holes`; minor. +- **`contrast_adjustment_clip_limit`** *(default 0.01)* — CLAHE aggressiveness; can help + low-contrast panels but interacts unpredictably with thresholding. + +**Tier 3 — high-value knob newly exposed by this component:** + +- **`watershed_compactness`** *(NEW; default 0.0)*. `skimage.watershed`'s compactness was + previously unreachable (only settable via the nested `watershed_params` dict that txsim + forwards but the component never populated). Now exposed and lifted into that dict + (`script.py:40-46`). 0.0 = standard watershed; small positive values (try 1e-3…1e-1) + enable **compact watershed** (Neubert & Protzel) for rounder, more regular cell shapes — + useful when watershed produces jagged/leaky basins. Left out of the default sweep to keep + it at 12 variants and honour the Tier-1 focus. **Requires `viash ns build` + a container + rebuild to take effect** (see "Risk points"). + +**Not worth touching for a benchmark:** the normalize/contrast function *choices*, all the +`distance_transform_*` return/indices/sampling args, the `threshold_shift_*`/`hist`/`out`/ +`mask` args, `blur_mode`/`cval`/`preserve_range`/`truncate`, `*_connectivity_*`, `*_out_*`. +They are either inert under the default functions (dropped by `filter_params`) or pure +plumbing. + +**Suggested quality-first order:** `local_maxima_min_distance` ↑ (fix over-splitting) → +`threshold_func` / `threshold_footprint_size` → `blur_sigma` ↑ → `post_processing_min_size_1` +↑ → `bg_intensity_filter_bg_factor` → `watershed_compactness` for shape. + +## Risk points / gotchas + +- **Newly exposed arg needs a rebuild.** `--watershed_compactness` was added to + `config.vsh.yaml` and threaded in `script.py`; it has **no effect until** `viash ns build` + regenerates the target and the Docker container is rebuilt (see the `check-component` + skill). The rest of the sweep uses pre-existing args and works against the current image. +- **Unmapped `*_func` string = silent skip or crash.** Only sweep `threshold_func` among + its three valid values; never point a selector at a name outside its mapping (Arguments + §). `threshold_func`/`local_maxima_func` set to an unmapped value → `NameError` because + `nuclei`/`local_maxima` never get defined. +- **`min_distance` too small over-splits; too large merges cells.** The default 5 errs + toward over-segmentation — this is the first knob to raise if masks look fragmented. +- **Whole image loaded into RAM** (`script.py:51`) — full-res plane; part of why the label + is `midmem` and the stage is `hightime`. +- **Only channel 0 is segmented** (`image[0]`) — other channels ignored. +- **Shared metadata block** (`script.py:59-78`) is duplicated across the segmentation + components; fix together. +- **Runtime is threshold/filter-bound**, not model-bound: `local_otsu` (rank filter, + footprint 50) and the background-intensity filter (windows 1000/2000) dominate. + +## Wiring + +- `src/workflows/run_benchmark/config.vsh.yaml`: listed as a dependency + (`methods_segmentation/watershed`, `:162`) and included in the default segmentation + method string `custom_segmentation:cellpose:cellposev4:binning:stardist:watershed` + (`:65`). +- `src/workflows/run_benchmark/main.nf:101`: imported. +- Run scripts (this tuning task): `scripts/run_benchmark/run_test_watershed_local.sh` and + `run_test_watershed_nebius.sh`, driven by the committed sweep + `scripts/run_benchmark/watershed_params.yaml` (single source of truth; local reads it via + `$REPO_ROOT`, Nebius via its raw GitHub URL — so it must be pushed before a Nebius + launch). Otherwise watershed is commented out in the standard `run_test_*` scripts. + +## References + +- **Watershed transform (the DOI in config):** Vincent L. & Soille P. (1991), "Watersheds + in digital spaces: an efficient algorithm based on immersion simulations", IEEE TPAMI + 13(6):583-598, DOI **10.1109/34.87344**. +- Implementation: `theislab/txsim` (`dev`), `txsim/preprocessing/_segmentation.py:382` + `segment_watershed`. +- Underlying ops: `skimage.segmentation.watershed` (marker-controlled; `compactness` per + Neubert & Protzel), `skimage.feature.peak_local_max`, `skimage.filters.rank.otsu`. diff --git a/src/methods_transcript_assignment/fastreseg/NOTES.md b/src/methods_transcript_assignment/fastreseg/NOTES.md new file mode 100644 index 000000000..e57e325a8 --- /dev/null +++ b/src/methods_transcript_assignment/fastreseg/NOTES.md @@ -0,0 +1,139 @@ +# fastReseg — transcript assignment + +## What this component is + +fastReseg is a **transcript-assignment** method (API: `src/api/comp_method_transcript_assignment.yaml`; +inputs `raw_ist.zarr` + `segmentation.zarr` + `scrnaseq_reference.h5ad`, output +`transcript_assignments.zarr`). It corrects an initial image-based segmentation using the +spatial profile of transcripts. Wraps the R package +[`Nanostring-Biostats/FastReseg`](https://github.com/Nanostring-Biostats/FastReseg) +(DOI 10.1038/s41598-025-08733-5). + +**What's unusual: it is the only multi-script, bash-orchestrated component in the repo.** +It has four resources and the *first* (`orchestrator.sh`, a `bash_script`) is the Viash +entrypoint. It chains three sub-scripts through TSV/CSV files in a temp dir: + +``` +orchestrator.sh + → input.py (SpatialData zarr → counts.tsv, transcripts.tsv, cell_types.tsv) + → script.R (FastReseg::fastReseg_full_pipeline → cell_ids.csv, gene_names.csv, transcripts_out.csv) + → output.py (CSVs → transcript_assignments.zarr) +``` + +## orchestrator.sh — control flow + +- `set -eo pipefail` (line 6) — **load-bearing.** Without it, a crash in any sub-step is + swallowed (the script's last line is `echo $(date)` → exit 0) and Viash reports only the + generic *"Required output file is missing. Expression: par.required"*, hiding the real + Python/R traceback. This misdirection cost a whole debugging round; keep it. +- Sub-scripts are invoked as `"$meta_resources_dir/