diff --git a/merxen_stub/.gitignore b/merxen_stub/.gitignore new file mode 100644 index 0000000..9e7764a --- /dev/null +++ b/merxen_stub/.gitignore @@ -0,0 +1,10 @@ +work/ +results/ +results_*/ +.nextflow* +test/work/ +test/results_*/ +test/.nextflow* +test/.samplesheet_*.csv +.DS_Store +**/.DS_Store diff --git a/merxen_stub/README.md b/merxen_stub/README.md new file mode 100644 index 0000000..e894abe --- /dev/null +++ b/merxen_stub/README.md @@ -0,0 +1,177 @@ +# merxen_stub + +A structural stub of [bourdenxlab/MerXen](https://github.com/bourdenxlab/MerXen) +for validating Flow compatibility (FLOW-700). + +The aim is to reproduce MerXen's **shape** — the same processes, the same +directory-typed inputs and outputs, the same `${outdir}////` +tree — while doing no real computation, so we can iterate quickly on Flow +(especially directory inputs/outputs) without GPUs, conda, or reference data. + +## What is faithful vs. stubbed + +- **Copied verbatim from MerXen:** `workflows/main.nf`, `workflows/nextflow.config`, + `workflows/conf/dwight.config`, `workflows/samplesheet.example.csv`. All the + Groovy orchestration, parameter parsing, channel wiring, the `SEGMENT` + sub-workflow, `publishDir` targets, and process `input:`/`output:` tuple + shapes are unchanged. +- **Stubbed:** only each process's `script:` body. The real `merxen …` / + `pip` / `cargo` / `python -m …` / GPU-monitor calls are replaced with + `mkdir`/`echo` that create exactly the process's declared outputs with + placeholder content. `.zarr` outputs are directories containing a dummy + `.zattrs`; `*_out/` outputs are directories with a placeholder file. +- **Not reproduced:** MerXen's side-channel writes to hardcoded absolute paths + under `outdir` (e.g. `//latest/latest_spatialdata.zarr`, + the durable `segmentation/*.npy`). These are written outside Nextflow's + `publishDir`, so Flow discovers outputs from process executions cannot target + them anyway. Every `publishDir` directory output *is* reproduced. + +- **Input access is verified, not faked.** MerXen receives the raw + `merscope_dir`/`xenium_dir` as a *value* path (embedded in a config JSON, not + a Nextflow-staged `path()`) and its Python opens that absolute path directly. + So `BUILD_SPATIALDATA` reproduces the real input resolution + (`io/builders/pipeline.py`): it reads `input_path` from the config and **fails + the task** if that path is not an accessible directory, or — for a raw export + folder — if it lacks the files the builder reads (MERSCOPE: `images/`, + `micron_to_mosaic_pixel_transform.csv`, `cell_boundaries.parquet`, + `detected_transcripts.parquet|csv`; XENIUM: `experiment.xenium`/`specs.json` + and `transcripts.parquet`). A reusable `.zarr` cache path only needs to be a + directory. This is what makes the stub a real probe of Flow directory inputs: + if Flow doesn't mount the selected data path into the process container, or + the uploaded directory has the wrong structure, the stub fails exactly where a + real run would — it does not report a false success. + +Per FLOW-700, the pipeline code is not modified to accommodate Flow — only the +schema (`flow/schema/merxen.json`) and the execution config a deployment +supplies (see `test/local.config`). + +The equivalent access check is *not* yet applied to the file inputs consumed by +later stages (MECR/clustering/MapMyCells reference paths, cortical-depth / +distance GeoJSONs). Those are value paths too, so they carry the same +container-mount risk; `main.nf`'s preflight checks their existence only on the +launch node, not inside the process container. + +## Layout + +``` +merxen_stub/ +├── flow/schema/merxen.json # Flow schema (inputs + directory outputs) +├── workflows/ # mirrors the MerXen repo's workflows/ dir +│ ├── main.nf # verbatim +│ ├── nextflow.config # verbatim +│ ├── conf/dwight.config # verbatim +│ └── modules/*.nf # verbatim except stubbed script bodies +└── test/ # local smoke test (not part of the pipeline) + ├── run_smoke.sh + ├── local.config # disables conda/containers/GPU-lock for local runs + ├── inputs/ # fake raw MERSCOPE/Xenium directories + └── refs/ # dummy reference files + annotation GeoJSONs +``` + +## Flow registration + +- **Path:** `merxen_stub/workflows/main.nf` +- **Schema Path:** `merxen_stub/flow/schema/merxen.json` + +The raw MERSCOPE/Xenium inputs are **directories**, exposed in the schema as +`data` columns on the samplesheet. MerXen's hardcoded reference paths (MECR, +clustering annotation, MapMyCells) are exposed as `data` parameters. Every +curated output is a directory (`"filetype": ""`). + +## Running on Flow + +Register (admin panel → Pipelines): + +- **Repo:** `goodwright/nf-test` +- **Path:** `merxen_stub/workflows/main.nf` +- **Schema Path:** `merxen_stub/flow/schema/merxen.json` +- **Config Paths:** `merxen_stub/flow/flow.config` + +`flow/flow.config` adapts MerXen's workstation defaults (GPU queues, GPU `flock` +locks, apptainer/conda directives, 200–300 GB memory requests from +`conf/dwight.config`) to a Flow docker + slurm environment: it runs the stub's +pure-shell bodies in `ubuntu:22.04`, drops the GPU queues/locks, and requests +trivial resources. It changes execution only, never pipeline logic. + +Two Flow-backend behaviours (from `flow-api`) must be handled first — these are +the actual FLOW-700 outcomes: + +### Nextflow version (blocker) + +Flow's default is `DEFAULT_NEXTFLOW_VERSION = 23.04.3`, which **cannot compile** +`main.nf` (`String too long`; see below). Set the pipeline version's +`allowed_nextflow_versions` (or the global `NextflowVersion` table) to a **26.x** +release, or the run fails before any process starts. + +### Output flattening / collisions (the core finding) + +Flow appends a generated `publish.config` **last** (highest precedence) that +forces `publishDir` for every process to a flat +`flow-results//`, and a `withName` publishDir in a later config +**replaces** a module's body `publishDir` (verified). So on a normal run: + +- MerXen's nested `${outdir}////` tree is **not** + produced — everything is captured flat, keyed only by process name. +- The `pair_id` / `platform` / `segmentation` context is lost from the paths, + so parallel tasks of the same process **collide**: both platforms' + `BUILD_SPATIALDATA` publish `spatialdata_out/source_spatialdata.zarr`, every + `QC` task publishes `qc_out/`, etc., all into the same `flow-results//`. + +MerXen relies on its nested `publishDir` paths to disambiguate outputs; Flow's +forced flat publishing collapses them. This is the compatibility gap FLOW-700 +set out to find. The backend does have an escape hatch: pipelines registered +with **"imports samples" = true** get an *empty* `publish.config` (native +publishing), which lets MerXen's own nested `publishDir` directives take effect. +That is the most promising next thing to try — register the pipeline with +"imports samples" checked and confirm the nested tree and per-pair/platform +disambiguation survive. + +Directory outputs themselves are supported: the backend records +`is_directory` from the filesystem snapshot and derives `filetype` from the +file extension (so a `.zarr` directory registers as a directory with filetype +`zarr`). `flow/flow.config` has the lowest precedence (config_paths are applied +before deployment configs and `publish.config`), so it only sets keys the +deployment leaves unset (container, GPU params, queue, resources) and cannot +influence `publishDir`. + +### Minimal first run + +Isolate directory handling from reference data with `stop_stage=enrich` and +`mecr_enabled=false`: build → segment → enrich needs only the raw +MERSCOPE/Xenium directory inputs and emits `.zarr`/`*_out` directories. Then add +the reference `data` inputs, push `stop_stage` to `clustering_squidpy`, and +finally enable the opt-in stages. + +The uploaded directory must have the structure `BUILD_SPATIALDATA` checks (see +"Input access is verified" above), or the run fails there **by design** — that +failure is the signal that Flow either didn't mount the path or delivered the +wrong structure. Use `test/inputs/merscope/EX01` and `test/inputs/xenium/EX01` +as correctly-structured dummy directories to upload. A single-platform +(`analysis_mode=xenium`) run needs just the one Xenium directory. + +## Running the smoke test + +```bash +merxen_stub/test/run_smoke.sh +``` + +Requires Nextflow (auto-pinned to `NXF_VER=26.04.6`; see caveat below) and +Java 17+. It runs two configurations and asserts the expected directory tree: + +- **Test A** — paired defaults, `build_spatialdata → clustering_squidpy`. +- **Test B** — paired with all opt-in stages on (`enable_alignment`, + `cortical_depth_enabled`, `distance_from_object_enabled`) through `mapmycells`. + +## Known behaviours surfaced while stubbing + +- **Nextflow version:** MerXen's `main.nf` is ~114 KB; its `workflow` body + compiles a >64 KB constant that the legacy Groovy parser rejects + (`String too long`) on Nextflow 24.x/25.x. It compiles and runs on 26.x. Flow + must run this pipeline on Nextflow ≥ 26.x. +- **Only the first segmentation branch runs.** With `analysis_segmentation=both` + (default), `main.nf` joins the enriched-zarr channel (one item per + `pair_id|platform`) against a per-segmentation gate (two items per key) using + `.join()`, which is 1:1 — so only `reseg` proceeds through the downstream + analysis stages and `original_seg` is silently dropped, despite + `docs/outputs.md` documenting both. This is upstream behaviour reproduced + faithfully; worth confirming against however MerXen is run in production. diff --git a/merxen_stub/flow/flow.config b/merxen_stub/flow/flow.config new file mode 100644 index 0000000..a7df65e --- /dev/null +++ b/merxen_stub/flow/flow.config @@ -0,0 +1,47 @@ +/* + * Flow execution config for the stubbed MerXen pipeline. + * + * Pass this via the pipeline version's "Config Paths". It adapts MerXen's + * workstation defaults (conf/dwight.config, pulled in by the default `standard` + * profile) to a Flow docker + slurm environment: the stub bodies are pure + * shell, so they run in a stock image rather than MerXen's GPU containers / + * conda envs, with the workstation GPU queues, GPU locks, and oversized + * resource requests removed. + * + * Precedence note: Flow applies config_paths configs FIRST (lowest priority), + * then its deployment configs, then a generated publish.config LAST. So this + * file cannot set publishDir (publish.config overrides it) and is overridden by + * any deployment config that sets the same keys. It only covers keys the + * deployment leaves unset (container, GPU params, queue, resources). See the + * README "Running on Flow" section for the publishDir/output-flattening + * consequences. + */ + +conda.enabled = false + +params { + cellpose_gpu = false + gpu_process_lock_enabled = false + clustering_squidpy_use_gpu = false + clustering_squidpy_gpu_vram_monitor = false + alignment_device = "cpu" +} + +process { + container = "ubuntu:22.04" + + withName: ".*" { + conda = null + beforeScript = "" + queue = null + clusterOptions = null + cpus = 1 + memory = 1.GB + } + + // MerXen pins this process to a GPU apptainer .sif; replace it with the + // stock image so the stub body runs under docker. + withName: "CLUSTERING_SQUIDPY_COMPUTE" { + container = "ubuntu:22.04" + } +} diff --git a/merxen_stub/flow/schema/merxen-nosamples.json b/merxen_stub/flow/schema/merxen-nosamples.json new file mode 100644 index 0000000..10e4f99 --- /dev/null +++ b/merxen_stub/flow/schema/merxen-nosamples.json @@ -0,0 +1,92 @@ +{ + "inputs": [ + { + "name": "Samples", + "description": "One row per pair_id. No Flow samples required: type a pair_id and pick a directory for each platform. Raw inputs are directories; the stub treats them as opaque paths, so any directory works for a dummy run.", + "params": { + "samplesheet": { + "name": "Samplesheet", + "description": "CSV with one row per pair_id.", + "required": true, + "type": "csv", + "takes_samples": false, + "allow_custom_columns": true, + "columns": [ + { + "name": "pair_id", + "type": "string", + "required": true, + "render": true + }, + { + "name": "analysis_mode", + "type": "string", + "valid": ["paired", "merscope", "xenium"], + "required": false, + "render": true + }, + { + "name": "merscope_dir", + "description": "Directory selected as data (required for paired/merscope modes).", + "type": "data", + "required": false, + "render": true + }, + { + "name": "xenium_dir", + "description": "Directory selected as data (required for paired/xenium modes).", + "type": "data", + "required": false, + "render": true + } + ] + } + } + }, + { + "name": "Analysis options", + "description": "For a dummy run, set Stop stage to 'enrich' and turn MECR off so no reference data is needed.", + "params": { + "stop_stage": { + "name": "Stop stage", + "type": "string", + "default": "enrich", + "required": false + }, + "mecr_enabled": { + "name": "Enable MECR", + "type": "boolean", + "default": false, + "required": false + }, + "analysis_segmentation": { + "name": "Analysis segmentation branches", + "type": "string", + "valid": ["both", "all", "reseg", "original_seg", "proseg_hybrid"], + "default": "both", + "required": false + } + } + } + ], + "outputs": [ + { + "name": "Source SpatialData", + "description": "Per-platform source SpatialData .zarr store.", + "process": "BUILD_SPATIALDATA", + "filetype": "" + }, + { + "name": "Segmentation", + "description": "Cellpose masks and ProSeg base SpatialData .zarr.", + "process": "PROSEG_SEGMENT", + "filetype": "" + }, + { + "name": "Enrichment", + "description": "Per-shape assignment summaries and updated latest SpatialData .zarr.", + "process": "ENRICH", + "filetype": "" + } + ] +} diff --git a/merxen_stub/flow/schema/merxen.json b/merxen_stub/flow/schema/merxen.json new file mode 100644 index 0000000..19814cd --- /dev/null +++ b/merxen_stub/flow/schema/merxen.json @@ -0,0 +1,307 @@ +{ + "inputs": [ + { + "name": "Samples", + "description": "One row per tissue block / adjacent-section pair. Raw platform inputs are directories: on Flow they are selected as data and materialised to a path the pipeline reads.", + "params": { + "samplesheet": { + "name": "Samplesheet", + "description": "CSV with one row per pair_id. Directory columns point at raw MERSCOPE/Xenium export folders or reusable SpatialData .zarr stores.", + "required": true, + "type": "csv", + "takes_samples": true, + "allow_custom_columns": true, + "columns": [ + { + "name": "pair_id", + "type": "string", + "from_sample": "name", + "required": true, + "render": true + }, + { + "name": "analysis_mode", + "type": "string", + "valid": ["paired", "merscope", "xenium"], + "required": false, + "render": true + }, + { + "name": "merscope_dir", + "description": "Raw MERSCOPE region export directory (transcripts.parquet, cell_boundaries/, images/, ...).", + "type": "data", + "required": false, + "render": true + }, + { + "name": "xenium_dir", + "description": "Raw Xenium export directory.", + "type": "data", + "required": false, + "render": true + }, + { + "name": "merscope_spatialdata_path", + "description": "Reusable MERSCOPE SpatialData .zarr directory (skips build when present).", + "type": "data", + "pattern": "zarr$", + "required": false, + "render": true + }, + { + "name": "xenium_spatialdata_path", + "description": "Reusable Xenium SpatialData .zarr directory.", + "type": "data", + "pattern": "zarr$", + "required": false, + "render": true + } + ] + } + } + }, + { + "name": "Reference data", + "description": "Assets MerXen hardcodes as absolute paths on the workstation. On Flow they are provided as data inputs. Required only for the stages that consume them (MECR, hierarchical clustering annotation, MapMyCells).", + "advanced": true, + "params": { + "mecr_neurons_h5ad_path": { + "name": "WHB neurons H5AD", + "description": "Whole-brain 10Xv3 neuron raw expression matrix (MECR reference).", + "type": "data", + "required": false + }, + "mecr_nonneurons_h5ad_path": { + "name": "WHB non-neurons H5AD", + "description": "Whole-brain 10Xv3 non-neuron raw expression matrix (MECR reference).", + "type": "data", + "required": false + }, + "mecr_cell_metadata_path": { + "name": "WHB cell metadata", + "description": "Whole-brain reference cell metadata CSV (MECR).", + "type": "data", + "required": false + }, + "mecr_taxonomy_metadata_path": { + "name": "WHB taxonomy metadata", + "description": "Cluster annotation term metadata CSV (MECR + clustering annotation).", + "type": "data", + "required": false + }, + "mecr_cluster_membership_path": { + "name": "WHB cluster membership", + "description": "Cluster-to-cluster-annotation membership metadata CSV.", + "type": "data", + "required": false + }, + "clustering_squidpy_broad_marker_lookup_path": { + "name": "Broad marker lookup", + "description": "Query-marker JSON used for hierarchical clustering broad annotation.", + "type": "data", + "required": false + }, + "clustering_squidpy_broad_taxonomy_metadata_path": { + "name": "Broad taxonomy metadata", + "description": "Taxonomy metadata CSV for clustering broad annotation.", + "type": "data", + "required": false + }, + "mapmycells_marker_lookup_path": { + "name": "MapMyCells marker lookup", + "description": "Whole-brain query-marker JSON for MapMyCells.", + "type": "data", + "required": false + }, + "mapmycells_precomputed_stats_path": { + "name": "MapMyCells precomputed stats", + "description": "Precomputed reference statistics for MapMyCells.", + "type": "data", + "required": false + } + } + }, + { + "name": "Analysis options", + "description": "Which downstream segmentation branches and optional stages run. These control which output directories are produced.", + "params": { + "analysis_segmentation": { + "name": "Analysis segmentation branches", + "description": "Downstream branch set. both=reseg,original_seg; all adds proseg_hybrid.", + "type": "string", + "valid": ["both", "all", "reseg", "original_seg", "proseg_hybrid"], + "default": "both", + "required": false + }, + "enable_alignment": { + "name": "Enable alignment", + "description": "Register paired sections with Spateo (paired mode only).", + "type": "boolean", + "default": false, + "required": false + }, + "mecr_enabled": { + "name": "Enable MECR", + "description": "Run mutually-exclusive co-expression rate analysis.", + "type": "boolean", + "default": true, + "required": false + }, + "mask_image_quantification_enabled": { + "name": "Enable mask image quantification", + "type": "boolean", + "default": true, + "required": false + }, + "viewer_cache_enabled": { + "name": "Enable viewer caches", + "type": "boolean", + "default": true, + "required": false + }, + "spatial_gene_analysis_enabled": { + "name": "Enable spatial gene analysis", + "type": "boolean", + "default": true, + "required": false + }, + "cortical_depth_enabled": { + "name": "Enable cortical depth", + "type": "boolean", + "default": false, + "required": false + }, + "distance_from_object_enabled": { + "name": "Enable distance from object", + "type": "boolean", + "default": false, + "required": false + }, + "start_stage": { + "name": "Start stage", + "description": "First stage to run.", + "type": "string", + "default": "build_spatialdata", + "required": false + }, + "stop_stage": { + "name": "Stop stage", + "description": "Final stage to run.", + "type": "string", + "default": "clustering_squidpy", + "required": false + } + } + } + ], + "outputs": [ + { + "name": "Source SpatialData", + "description": "Per-platform source SpatialData .zarr store built from raw data or symlinked from cache.", + "process": "BUILD_SPATIALDATA", + "filetype": "" + }, + { + "name": "Segmentation", + "description": "Cellpose masks, ProSeg base SpatialData .zarr, and stitching diagnostics.", + "process": "PROSEG_SEGMENT", + "filetype": "" + }, + { + "name": "Enrichment", + "description": "Per-shape assignment summaries and the updated latest SpatialData .zarr.", + "process": "ENRICH", + "filetype": "" + }, + { + "name": "Viewer cache", + "description": "Napari viewer derived caches written into the latest SpatialData .zarr.", + "process": "VIEWER_CACHE", + "filetype": "" + }, + { + "name": "Mask image quantification", + "description": "Per-channel Cellpose cell x image-stat matrices.", + "process": "MASK_IMAGE_QUANTIFICATION", + "filetype": "" + }, + { + "name": "QC", + "description": "Per-dataset, per-segmentation QC metrics directory.", + "process": "QC", + "filetype": "" + }, + { + "name": "MECR reference", + "description": "Shared whole-brain reference markers and statistics.", + "process": "MECR_REFERENCE", + "filetype": "" + }, + { + "name": "MECR", + "description": "Per-pair, per-segmentation mutually-exclusive co-expression results.", + "process": "MECR", + "filetype": "" + }, + { + "name": "Alignment", + "description": "Spateo alignment transform and centroid coordinate tables.", + "process": "ALIGN", + "filetype": "" + }, + { + "name": "Alignment QC", + "description": "Alignment QC metrics and overlays.", + "process": "ALIGN_QC", + "filetype": "" + }, + { + "name": "Comparison", + "description": "Paired MERSCOPE vs Xenium gene-count comparison tables.", + "process": "COMPARE", + "filetype": "" + }, + { + "name": "Visualization", + "description": "Paired or single-platform visualisation plots.", + "process": "VISUALIZE", + "filetype": "" + }, + { + "name": "Spatial gene analysis", + "description": "Spatial autocorrelation and transcript pattern outputs.", + "process": "SPATIAL_GENE_ANALYSIS", + "filetype": "" + }, + { + "name": "Squidpy clustering", + "description": "Clustered AnnData, UMAP/spatial plots, and hierarchical annotation.", + "process": "CLUSTERING_SQUIDPY_FINALIZE", + "filetype": "" + }, + { + "name": "Cortical depth", + "description": "Laplace/equivolumetric depth outputs (opt-in).", + "process": "COMPUTE_CORTICAL_DEPTH", + "filetype": "" + }, + { + "name": "Distance from object", + "description": "Per-block nearest-object distance outputs (opt-in).", + "process": "DISTANCE_FROM_OBJECT_ANNOTATE", + "filetype": "" + }, + { + "name": "Distance from object cohort", + "description": "Cohort paired near-vs-far differential expression (opt-in).", + "process": "DISTANCE_FROM_OBJECT_COHORT", + "filetype": "" + }, + { + "name": "MapMyCells", + "description": "Per-cell cell-type assignments and QC (opt-in).", + "process": "MAPMYCELLS", + "filetype": "" + } + ] +} diff --git a/merxen_stub/test/.gitignore b/merxen_stub/test/.gitignore new file mode 100644 index 0000000..5931300 --- /dev/null +++ b/merxen_stub/test/.gitignore @@ -0,0 +1,4 @@ +work/ +results_*/ +.nextflow* +.samplesheet_*.csv diff --git a/merxen_stub/test/inputs/merscope/EX01/cell_boundaries.parquet b/merxen_stub/test/inputs/merscope/EX01/cell_boundaries.parquet new file mode 100644 index 0000000..a9aee81 --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/cell_boundaries.parquet @@ -0,0 +1 @@ +stub boundaries diff --git a/merxen_stub/test/inputs/merscope/EX01/cell_by_gene.csv b/merxen_stub/test/inputs/merscope/EX01/cell_by_gene.csv new file mode 100644 index 0000000..0ffddf6 --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/cell_by_gene.csv @@ -0,0 +1 @@ +gene,count diff --git a/merxen_stub/test/inputs/merscope/EX01/cell_metadata.csv b/merxen_stub/test/inputs/merscope/EX01/cell_metadata.csv new file mode 100644 index 0000000..896990d --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/cell_metadata.csv @@ -0,0 +1 @@ +cell,x,y diff --git a/merxen_stub/test/inputs/merscope/EX01/detected_transcripts.parquet b/merxen_stub/test/inputs/merscope/EX01/detected_transcripts.parquet new file mode 100644 index 0000000..fb446d3 --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/detected_transcripts.parquet @@ -0,0 +1 @@ +stub detected transcripts diff --git a/merxen_stub/test/inputs/merscope/EX01/images/mosaic_DAPI_z0.tif b/merxen_stub/test/inputs/merscope/EX01/images/mosaic_DAPI_z0.tif new file mode 100644 index 0000000..13b6ab1 --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/images/mosaic_DAPI_z0.tif @@ -0,0 +1 @@ +stub image diff --git a/merxen_stub/test/inputs/merscope/EX01/micron_to_mosaic_pixel_transform.csv b/merxen_stub/test/inputs/merscope/EX01/micron_to_mosaic_pixel_transform.csv new file mode 100644 index 0000000..cc8404f --- /dev/null +++ b/merxen_stub/test/inputs/merscope/EX01/micron_to_mosaic_pixel_transform.csv @@ -0,0 +1 @@ +col,row diff --git a/merxen_stub/test/inputs/xenium/EX01/cell_boundaries.parquet b/merxen_stub/test/inputs/xenium/EX01/cell_boundaries.parquet new file mode 100644 index 0000000..3a248d5 --- /dev/null +++ b/merxen_stub/test/inputs/xenium/EX01/cell_boundaries.parquet @@ -0,0 +1 @@ +stub cell boundaries diff --git a/merxen_stub/test/inputs/xenium/EX01/cells.parquet b/merxen_stub/test/inputs/xenium/EX01/cells.parquet new file mode 100644 index 0000000..ffc20bc --- /dev/null +++ b/merxen_stub/test/inputs/xenium/EX01/cells.parquet @@ -0,0 +1 @@ +stub cells diff --git a/merxen_stub/test/inputs/xenium/EX01/experiment.xenium b/merxen_stub/test/inputs/xenium/EX01/experiment.xenium new file mode 100644 index 0000000..8fb8eeb --- /dev/null +++ b/merxen_stub/test/inputs/xenium/EX01/experiment.xenium @@ -0,0 +1 @@ +stub xenium experiment diff --git a/merxen_stub/test/inputs/xenium/EX01/transcripts.parquet b/merxen_stub/test/inputs/xenium/EX01/transcripts.parquet new file mode 100644 index 0000000..53865bf --- /dev/null +++ b/merxen_stub/test/inputs/xenium/EX01/transcripts.parquet @@ -0,0 +1 @@ +stub xenium transcripts diff --git a/merxen_stub/test/local.config b/merxen_stub/test/local.config new file mode 100644 index 0000000..79976be --- /dev/null +++ b/merxen_stub/test/local.config @@ -0,0 +1,30 @@ +/* + * Local execution overrides for smoke-testing the stubbed pipeline. + * + * This file is NOT part of the pipeline. It supplies the execution-environment + * settings a real deployment (Flow, an HPC profile) would provide: it disables + * conda/containers and the GPU process locks so the stubbed bash bodies run + * natively on a workstation with no MerXen dependencies installed. + */ + +conda.enabled = false + +params { + cellpose_gpu = false + gpu_process_lock_enabled = false + clustering_squidpy_use_gpu = false + clustering_squidpy_gpu_vram_monitor = false + alignment_device = "cpu" +} + +process.executor = "local" + +// dwight.config attaches a GPU-lock beforeScript (flock on a shared lock file) +// to the GPU processes. flock is Linux-only and the lock is a workstation +// concern, so clear it for local stub runs. +process { + withName: "CELLPOSE_SEGMENT" { beforeScript = "" } + withName: "CELLPOSE_NUCLEI_SEGMENT" { beforeScript = "" } + withName: "ALIGN" { beforeScript = "" } + withName: "CLUSTERING_SQUIDPY_COMPUTE" { beforeScript = "" } +} diff --git a/merxen_stub/test/refs/EX01_merscope_objects.geojson b/merxen_stub/test/refs/EX01_merscope_objects.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_merscope_objects.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/EX01_merscope_pia.geojson b/merxen_stub/test/refs/EX01_merscope_pia.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_merscope_pia.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/EX01_merscope_side.geojson b/merxen_stub/test/refs/EX01_merscope_side.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_merscope_side.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/EX01_xenium_objects.geojson b/merxen_stub/test/refs/EX01_xenium_objects.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_xenium_objects.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/EX01_xenium_pia.geojson b/merxen_stub/test/refs/EX01_xenium_pia.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_xenium_pia.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/EX01_xenium_side.geojson b/merxen_stub/test/refs/EX01_xenium_side.geojson new file mode 100644 index 0000000..66b7952 --- /dev/null +++ b/merxen_stub/test/refs/EX01_xenium_side.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[]} diff --git a/merxen_stub/test/refs/precomputed_stats.h5 b/merxen_stub/test/refs/precomputed_stats.h5 new file mode 100644 index 0000000..af8eafb --- /dev/null +++ b/merxen_stub/test/refs/precomputed_stats.h5 @@ -0,0 +1 @@ +stub reference precomputed_stats.h5 diff --git a/merxen_stub/test/refs/query_markers.json b/merxen_stub/test/refs/query_markers.json new file mode 100644 index 0000000..58f39eb --- /dev/null +++ b/merxen_stub/test/refs/query_markers.json @@ -0,0 +1 @@ +stub reference query_markers.json diff --git a/merxen_stub/test/refs/whb_cell_metadata.csv b/merxen_stub/test/refs/whb_cell_metadata.csv new file mode 100644 index 0000000..28714bc --- /dev/null +++ b/merxen_stub/test/refs/whb_cell_metadata.csv @@ -0,0 +1 @@ +stub reference whb_cell_metadata.csv diff --git a/merxen_stub/test/refs/whb_cluster_membership.csv b/merxen_stub/test/refs/whb_cluster_membership.csv new file mode 100644 index 0000000..1b17a45 --- /dev/null +++ b/merxen_stub/test/refs/whb_cluster_membership.csv @@ -0,0 +1 @@ +stub reference whb_cluster_membership.csv diff --git a/merxen_stub/test/refs/whb_neurons.h5ad b/merxen_stub/test/refs/whb_neurons.h5ad new file mode 100644 index 0000000..1d5f95e --- /dev/null +++ b/merxen_stub/test/refs/whb_neurons.h5ad @@ -0,0 +1 @@ +stub reference whb_neurons.h5ad diff --git a/merxen_stub/test/refs/whb_nonneurons.h5ad b/merxen_stub/test/refs/whb_nonneurons.h5ad new file mode 100644 index 0000000..fdc05ce --- /dev/null +++ b/merxen_stub/test/refs/whb_nonneurons.h5ad @@ -0,0 +1 @@ +stub reference whb_nonneurons.h5ad diff --git a/merxen_stub/test/refs/whb_taxonomy.csv b/merxen_stub/test/refs/whb_taxonomy.csv new file mode 100644 index 0000000..d6b73be --- /dev/null +++ b/merxen_stub/test/refs/whb_taxonomy.csv @@ -0,0 +1 @@ +stub reference whb_taxonomy.csv diff --git a/merxen_stub/test/run_smoke.sh b/merxen_stub/test/run_smoke.sh new file mode 100755 index 0000000..25a63fc --- /dev/null +++ b/merxen_stub/test/run_smoke.sh @@ -0,0 +1,116 @@ +#!/usr/bin/env bash +# +# Local smoke test for the stubbed MerXen pipeline. +# +# Runs the pipeline in two configurations and asserts that the expected +# directory-typed outputs appear. Verifies the DAG wiring and the output tree, +# not any scientific result (every process body is stubbed). +set -euo pipefail + +# MerXen's main.nf is large; it only compiles on Nextflow's newer language +# parser (>= 26.x). Older versions fail with "String too long". +export NXF_VER="${NXF_VER:-26.04.6}" + +TEST_DIR="$(cd "$(dirname "$0")" && pwd)" +STUB_DIR="$(cd "$TEST_DIR/.." && pwd)" +MAIN_NF="$STUB_DIR/workflows/main.nf" +LOCAL_CONFIG="$TEST_DIR/local.config" +REFS="$TEST_DIR/refs" +INPUTS="$TEST_DIR/inputs" +WORK="$TEST_DIR/work" + +PASS=0 +FAIL=0 + +report() { + if [ "$1" = "PASS" ]; then echo " PASS: $2"; PASS=$((PASS + 1)) + else echo " FAIL: $2"; FAIL=$((FAIL + 1)); fi +} + +assert_dir() { + [ -d "$1" ] && report "PASS" "$2" || report "FAIL" "$2 (missing $1)" +} + +# Reference paths that MerXen hardcodes; on Flow these become data inputs. +ref_args=( + --mecr_neurons_h5ad_path="$REFS/whb_neurons.h5ad" + --mecr_nonneurons_h5ad_path="$REFS/whb_nonneurons.h5ad" + --mecr_cell_metadata_path="$REFS/whb_cell_metadata.csv" + --mecr_taxonomy_metadata_path="$REFS/whb_taxonomy.csv" + --mecr_cluster_membership_path="$REFS/whb_cluster_membership.csv" + --clustering_squidpy_broad_marker_lookup_path="$REFS/query_markers.json" + --clustering_squidpy_broad_taxonomy_metadata_path="$REFS/whb_taxonomy.csv" + --clustering_squidpy_broad_cluster_membership_path="$REFS/whb_cluster_membership.csv" + --mapmycells_marker_lookup_path="$REFS/query_markers.json" + --mapmycells_precomputed_stats_path="$REFS/precomputed_stats.h5" +) + +run_pipeline() { + local outdir="$1"; shift + local samplesheet="$1"; shift + rm -rf "$outdir" "$WORK" + nextflow -log "$TEST_DIR/.nextflow.log" run "$MAIN_NF" \ + -c "$LOCAL_CONFIG" \ + -work-dir "$WORK" \ + --samplesheet "$samplesheet" \ + --outdir "$outdir" \ + "${ref_args[@]}" \ + "$@" +} + +# --------------------------------------------------------------------------- +echo "=== Test A: paired defaults (build -> clustering_squidpy) ===" +SS_A="$TEST_DIR/.samplesheet_a.csv" +{ + echo "pair_id,merscope_dir,xenium_dir" + echo "EX01,$INPUTS/merscope/EX01,$INPUTS/xenium/EX01" +} > "$SS_A" +OUT_A="$TEST_DIR/results_a" +if run_pipeline "$OUT_A" "$SS_A" --spatial_gene_analysis_transcript_analysis_enabled false; then + assert_dir "$OUT_A/EX01/merscope/spatialdata/spatialdata_out/source_spatialdata.zarr" "MERSCOPE source spatialdata zarr" + assert_dir "$OUT_A/EX01/xenium/spatialdata/spatialdata_out/source_spatialdata.zarr" "XENIUM source spatialdata zarr" + assert_dir "$OUT_A/EX01/merscope/enrichment/enrich_out" "MERSCOPE enrich_out dir" + assert_dir "$OUT_A/EX01/merscope/viewer_cache/viewer_cache_out" "MERSCOPE viewer_cache_out dir" + assert_dir "$OUT_A/EX01/merscope/mask_image_quantification/mask_image_quantification_out" "MERSCOPE mask_image_quantification_out dir" + assert_dir "$OUT_A/EX01/merscope/segmentation/segment_out/proseg_base_latest.zarr" "MERSCOPE proseg base latest zarr" + assert_dir "$OUT_A/EX01/merscope/reseg/qc/qc_out" "MERSCOPE reseg qc_out dir" + assert_dir "$OUT_A/EX01/xenium/reseg/qc/qc_out" "XENIUM reseg qc_out dir" + assert_dir "$OUT_A/mecr_reference/mecr_reference_out" "shared MECR reference dir" + assert_dir "$OUT_A/EX01/reseg/mecr/mecr_out" "reseg mecr_out dir" + assert_dir "$OUT_A/EX01/reseg/comparison/compare_out" "reseg compare_out dir (paired)" + assert_dir "$OUT_A/EX01/reseg/visualization/visualize_out" "reseg visualize_out dir" + assert_dir "$OUT_A/EX01/reseg/spatial_gene_analysis/spatial_gene_analysis_out" "reseg spatial_gene_analysis_out dir" + assert_dir "$OUT_A/EX01/reseg/clustering_squidpy/clustering_squidpy_out" "reseg clustering_squidpy_out dir" + assert_dir "$OUT_A/EX01/xenium/enrichment/latest_input.zarr" "XENIUM enrichment latest_input.zarr (symlinked dir)" +else + report "FAIL" "pipeline runs with paired defaults" +fi + +# --------------------------------------------------------------------------- +echo "" +echo "=== Test B: paired, all opt-in stages (-> mapmycells) ===" +SS_B="$TEST_DIR/.samplesheet_b.csv" +{ + echo "pair_id,merscope_dir,xenium_dir,merscope_pial_boundary_geojson,merscope_side_boundaries_geojson,xenium_pial_boundary_geojson,xenium_side_boundaries_geojson,merscope_distance_object_annotation_geojson,xenium_distance_object_annotation_geojson" + echo "EX01,$INPUTS/merscope/EX01,$INPUTS/xenium/EX01,$REFS/EX01_merscope_pia.geojson,$REFS/EX01_merscope_side.geojson,$REFS/EX01_xenium_pia.geojson,$REFS/EX01_xenium_side.geojson,$REFS/EX01_merscope_objects.geojson,$REFS/EX01_xenium_objects.geojson" +} > "$SS_B" +OUT_B="$TEST_DIR/results_b" +if run_pipeline "$OUT_B" "$SS_B" \ + --enable_alignment true \ + --cortical_depth_enabled true \ + --distance_from_object_enabled true \ + --stop_stage mapmycells; then + assert_dir "$OUT_B/EX01/alignment/align_out/alignment_coords" "alignment coords dir (paired)" + assert_dir "$OUT_B/EX01/alignment_qc/alignment_qc_out" "alignment_qc_out dir" + assert_dir "$OUT_B/EX01/merscope/compute_cortical_depth/compute_cortical_depth_out" "MERSCOPE compute_cortical_depth_out dir" + assert_dir "$OUT_B/EX01/merscope/distance_from_object/distance_from_object_out" "MERSCOPE distance_from_object_out dir" + assert_dir "$OUT_B/distance_from_object/cohort/merscope/distance_from_object_cohort_out" "MERSCOPE distance cohort dir" + assert_dir "$OUT_B/EX01/reseg/mapmycells/mapmycells_out" "reseg mapmycells_out dir" +else + report "FAIL" "pipeline runs with all opt-in stages" +fi + +rm -rf "$WORK" +echo "" +echo "=== Results: $PASS passed, $FAIL failed ===" +[ "$FAIL" -eq 0 ] diff --git a/merxen_stub/workflows/conf/dwight.config b/merxen_stub/workflows/conf/dwight.config new file mode 100644 index 0000000..fe10301 --- /dev/null +++ b/merxen_stub/workflows/conf/dwight.config @@ -0,0 +1,152 @@ +/* + * Execution settings for the-dwight workstation. + * + * Scientific and algorithmic defaults remain in ../nextflow.config. This file + * contains host paths, hardware choices, concurrency limits, and local + * executor capacity that a future HPC profile will replace. + */ + +params { + // Local software and GPU availability + cellpose_gpu = true + proseg_binary = "/usr/local/bin/proseg" + proseg_search_paths = [ + "/usr/bin/proseg", + "/usr/local/bin/proseg", + ] + proseg_install_path = "/usr/local/bin/proseg" + proseg_num_threads = 32 + alignment_device = "auto" + alignment_pytorch_cuda_alloc_conf = "expandable_segments:True,max_split_size_mb:256" + clustering_squidpy_use_gpu = true + clustering_squidpy_gpu_vram_monitor = true + clustering_squidpy_gpu_vram_monitor_interval_seconds = 2 + + // Host-level concurrency + mask_image_quantification_max_forks = 3 + // Nine 60 GB / 8 CPU reservations use at most 540 GB and 72 CPUs, + // leaving a conservative margin on the 775 GB / 80-thread workstation. + viewer_cache_max_forks = 9 + cortical_depth_max_forks = 3 + distance_from_object_n_cpus = 8 + distance_from_object_max_forks = 3 + alignment_max_forks = 1 + mecr_max_forks = 4 + clustering_squidpy_max_forks = 4 + spatial_gene_analysis_max_forks = 4 + + // Local reference assets + mecr_neurons_h5ad_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/expression_matrices/WHB-10Xv3/20240330/WHB-10Xv3-Neurons-raw.h5ad" + mecr_nonneurons_h5ad_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/expression_matrices/WHB-10Xv3/20240330/WHB-10Xv3-Nonneurons-raw.h5ad" + mecr_cell_metadata_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/metadata/WHB-10Xv3/20241115/cell_metadata.csv" + mecr_taxonomy_metadata_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/metadata/WHB-taxonomy/20240330/cluster_annotation_term.csv" + mecr_cluster_membership_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/metadata/WHB-taxonomy/20240330/cluster_to_cluster_annotation_membership.csv" + clustering_squidpy_broad_marker_lookup_path = "/media/mathieubo/SSD1/MerXen/mapmycells/query_markers.n10.20240221800.json" + clustering_squidpy_broad_taxonomy_metadata_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/metadata/WHB-taxonomy/20240330/cluster_annotation_term.csv" + clustering_squidpy_broad_cluster_membership_path = "/media/mathieubo/SSD1/MerXen/mapmycells/abc_whb/metadata/WHB-taxonomy/20240330/cluster_to_cluster_annotation_membership.csv" + clustering_squidpy_broad_reference_cache_dir = "/media/mathieubo/SSD1/MerXen/mapmycells" + mapmycells_marker_lookup_path = "/media/mathieubo/SSD1/MerXen/mapmycells/query_markers.n10.20240221800.json" + mapmycells_precomputed_stats_path = "/media/mathieubo/SSD1/MerXen/mapmycells/precomputed_stats.siletti.training.h5" + mapmycells_region_cache_dir = "/media/mathieubo/SSD1/MerXen/mapmycells" + + // Workstation memory and scheduling limits + max_ram_gb = 640 + warn_ram_gb = 600 + transcript_chunk_rows = 1000000 + build_spatialdata_max_forks = 3 + cellpose_segment_max_forks = 1 + proseg_segment_max_forks = 2 + gpu_process_lock_enabled = true + // All tasks and concurrent launches on Dwight share one GPU lock. + gpu_process_lock_file = "/tmp/merxen-dwight-gpu.lock" + mapmycells_n_processors = 8 +} + +process.executor = "local" + +// Advertise host capacity to the local executor without turning it into the +// default resource request for every process. +executor { + cpus = 72 + memory = "640 GB" +} + +process { + withName: "BUILD_SPATIALDATA" { + maxForks = params.build_spatialdata_max_forks + } + withName: "CELLPOSE_SEGMENT" { + maxForks = params.cellpose_segment_max_forks + beforeScript = """ + if [[ "${params.gpu_process_lock_enabled}" == "true" && "${params.cellpose_gpu}" == "true" ]]; then + MERXEN_GPU_LOCK_FILE="${params.gpu_process_lock_file}" + mkdir -p "\$(dirname "\${MERXEN_GPU_LOCK_FILE}")" + exec 9>"\${MERXEN_GPU_LOCK_FILE}" + echo "Waiting for MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + flock 9 + echo "Acquired MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + fi + """.stripIndent().trim() + } + withName: "CELLPOSE_NUCLEI_SEGMENT" { + maxForks = params.cellpose_segment_max_forks + beforeScript = """ + if [[ "${params.gpu_process_lock_enabled}" == "true" && "${params.cellpose_gpu}" == "true" ]]; then + MERXEN_GPU_LOCK_FILE="${params.gpu_process_lock_file}" + mkdir -p "\$(dirname "\${MERXEN_GPU_LOCK_FILE}")" + exec 9>"\${MERXEN_GPU_LOCK_FILE}" + echo "Waiting for MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + flock 9 + echo "Acquired MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + fi + """.stripIndent().trim() + } + withName: "PROSEG_SEGMENT" { + maxForks = params.proseg_segment_max_forks + } + withName: "VIEWER_CACHE" { + memory = "60 GB" + maxForks = params.viewer_cache_max_forks + } + withName: "MASK_IMAGE_QUANTIFICATION" { + maxForks = params.mask_image_quantification_max_forks + } + withName: "COMPUTE_CORTICAL_DEPTH" { + maxForks = params.cortical_depth_max_forks + } + withName: "DISTANCE_FROM_OBJECT_ANNOTATE" { + maxForks = params.distance_from_object_max_forks + } + withName: "ALIGN" { + maxForks = params.alignment_max_forks + beforeScript = """ + if [[ "${params.gpu_process_lock_enabled}" == "true" && "${params.alignment_device}" != "cpu" ]]; then + MERXEN_GPU_LOCK_FILE="${params.gpu_process_lock_file}" + mkdir -p "\$(dirname "\${MERXEN_GPU_LOCK_FILE}")" + exec 9>"\${MERXEN_GPU_LOCK_FILE}" + echo "Waiting for MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + flock 9 + echo "Acquired MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + fi + """.stripIndent().trim() + } + withName: "CLUSTERING_SQUIDPY_COMPUTE" { + maxForks = params.clustering_squidpy_max_forks + beforeScript = """ + if [[ "${params.gpu_process_lock_enabled}" == "true" && "${params.clustering_squidpy_use_gpu}" == "true" ]]; then + MERXEN_GPU_LOCK_FILE="${params.gpu_process_lock_file}" + mkdir -p "\$(dirname "\${MERXEN_GPU_LOCK_FILE}")" + exec 9>"\${MERXEN_GPU_LOCK_FILE}" + echo "Waiting for MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + flock 9 + echo "Acquired MerXen GPU process lock: \${MERXEN_GPU_LOCK_FILE}" >&2 + fi + """.stripIndent().trim() + } + withName: "SPATIAL_GENE_ANALYSIS" { + maxForks = params.spatial_gene_analysis_max_forks + } + withName: "MECR" { + maxForks = params.mecr_max_forks + } +} diff --git a/merxen_stub/workflows/main.nf b/merxen_stub/workflows/main.nf new file mode 100644 index 0000000..4df84b1 --- /dev/null +++ b/merxen_stub/workflows/main.nf @@ -0,0 +1,3057 @@ +nextflow.enable.dsl = 2 + +include { BUILD_SPATIALDATA } from "./modules/spatialdata_build" +include { ENSURE_PROSEG } from "./modules/proseg_bootstrap" +include { CELLPOSE_NUCLEI_SEGMENT; SEGMENT } from "./modules/segmentation" +include { ENRICH } from "./modules/enrichment" +include { VIEWER_CACHE } from "./modules/viewer_cache" +include { MASK_IMAGE_QUANTIFICATION } from "./modules/mask_image_quantification" +include { COMPUTE_CORTICAL_DEPTH } from "./modules/compute_cortical_depth" +include { + DISTANCE_FROM_OBJECT_ANNOTATE; + DISTANCE_FROM_OBJECT_COHORT +} from "./modules/distance_from_object" +include { VALIDATE_ANALYSIS_LAYER; QC } from "./modules/qc" +include { ALIGN; ALIGN_QC } from "./modules/alignment" +include { COMPARE } from "./modules/comparison" +include { VISUALIZE } from "./modules/visualization" +include { SPATIAL_GENE_ANALYSIS } from "./modules/spatial_gene_analysis" +include { MECR_REFERENCE; MECR } from "./modules/mecr" +include { + CLUSTERING_SQUIDPY_PREPARE; + CLUSTERING_SQUIDPY_COMPUTE; + CLUSTERING_SQUIDPY_FINALIZE +} from "./modules/clustering_squidpy" +include { MAPMYCELLS } from "./modules/mapmycells" + +def parseChannels(rawValue, defaults) { + if (rawValue == null) { + return defaults + } + def values = rawValue + .toString() + .split(",") + .collect { value -> value.trim() } + .findAll { value -> value.length() > 0 } + return values ? values : defaults +} + +def parseRange(rawValue, fallbackStart = 0, fallbackEnd = 6) { + if (rawValue == null || rawValue.toString().trim().isEmpty()) { + return [fallbackStart as int, fallbackEnd as int] + } + def parts = rawValue.toString().split("-").collect { value -> value.trim() } + if (parts.size() != 2) { + return [fallbackStart as int, fallbackEnd as int] + } + return [parts[0] as int, parts[1] as int] +} + +def intOrDefault(rawValue, defaultValue) { + if (rawValue == null || rawValue.toString().trim().isEmpty()) { + return defaultValue as int + } + return rawValue as int +} + +def floatOrDefault(rawValue, defaultValue) { + if (rawValue == null || rawValue.toString().trim().isEmpty()) { + return defaultValue as float + } + return rawValue as float +} + +def chooseField(row, names) { + def result = names.find { name -> + def value = row[name] + value != null && value.toString().trim().length() > 0 + } + return result ? row[result].toString().trim() : null +} + +def normalizeAnalysisMode(rawValue) { + def raw = rawValue == null ? "paired" : rawValue.toString().trim() + if (!raw) { + raw = "paired" + } + def key = raw + .toLowerCase() + .replaceAll(/[^a-z0-9]+/, "_") + .replaceAll(/^_+|_+$/, "") + def aliases = [ + "paired": "paired", + "pair": "paired", + "both": "paired", + "merscope": "merscope", + "merfish": "merscope", + "m": "merscope", + "xenium": "xenium", + "x": "xenium", + ] + if (!aliases.containsKey(key)) { + throw new IllegalArgumentException( + "Unknown analysis_mode '${raw}'. Valid values: paired, merscope, xenium" + ) + } + return aliases[key] +} + +def activePlatformsForMode(analysisMode) { + if (analysisMode == "paired") { + return ["MERSCOPE", "XENIUM"] + } + if (analysisMode == "merscope") { + return ["MERSCOPE"] + } + if (analysisMode == "xenium") { + return ["XENIUM"] + } + throw new IllegalArgumentException("Unknown analysis mode: ${analysisMode}") +} + +def normalizeAnalysisSegmentation(rawValue) { + def raw = rawValue == null ? "both" : rawValue.toString().trim() + if (!raw) { + raw = "both" + } + def aliases = [ + "both": ["reseg", "original_seg"], + "all": ["reseg", "original_seg", "proseg_hybrid"], + "reseg": ["reseg"], + "resegmented": ["reseg"], + "proseg": ["reseg"], + "mosaik": ["reseg"], + "hybrid": ["proseg_hybrid"], + "proseg_hybrid": ["proseg_hybrid"], + "hybrid_seg": ["proseg_hybrid"], + "original": ["original_seg"], + "original_seg": ["original_seg"], + "original_segmentation": ["original_seg"], + "instrument": ["original_seg"], + "instrument_seg": ["original_seg"], + "instrument_segmentation": ["original_seg"], + ] + def selected = [] + raw + .split(",") + .collect { value -> value.trim() } + .findAll { value -> value.length() > 0 } + .each { value -> + def key = value + .toLowerCase() + .replaceAll(/[^a-z0-9]+/, "_") + .replaceAll(/^_+|_+$/, "") + if (!aliases.containsKey(key)) { + throw new IllegalArgumentException( + "Unknown analysis_segmentation '${value}'. Valid values: " + + "both, all, reseg, original_seg, proseg_hybrid" + ) + } + aliases[key].each { segmentation -> + if (!selected.contains(segmentation)) { + selected << segmentation + } + } + } + return selected ? selected : ["reseg", "original_seg"] +} + +def requirePlatformInput(row, pairId, platform) { + if (platform == "MERSCOPE") { + def merscopeDir = chooseField(row, ["merscope_dir"]) + def merscopeSpatialdataPath = chooseField( + row, + ["merscope_spatialdata_path", "merscope_zarr_path"] + ) + if (!merscopeDir && !merscopeSpatialdataPath) { + error( + "Samplesheet row for ${pairId} must provide " + + "merscope_dir or merscope_spatialdata_path" + ) + } + return [inputDir: merscopeDir, spatialdataPath: merscopeSpatialdataPath] + } + + if (platform == "XENIUM") { + def xeniumDir = chooseField(row, ["xenium_dir"]) + def xeniumSpatialdataPath = chooseField(row, ["xenium_spatialdata_path"]) + if (!xeniumDir && !xeniumSpatialdataPath) { + error( + "Samplesheet row for ${pairId} must provide " + + "xenium_dir or xenium_spatialdata_path" + ) + } + return [inputDir: xeniumDir, spatialdataPath: xeniumSpatialdataPath] + } + + throw new IllegalArgumentException("Unknown platform: ${platform}") +} + +def buildConfigForPlatform(row, pairId, platform) { + def input = requirePlatformInput(row, pairId, platform) + if (platform == "MERSCOPE") { + def zRange = parseRange(row.merscope_z_range, 0, 6) + return [ + dataset_name: "${pairId}_MERSCOPE", + platform: "MERSCOPE", + input_path: input.inputDir ?: input.spatialdataPath, + output_path: "spatialdata_out/source_spatialdata.zarr", + persistent_output_path: input.spatialdataPath ?: null, + merscope_transform_path: chooseField(row, ["merscope_transform_path"]) ?: null, + merscope: [ + z_layers: (zRange[0]..zRange[1]).collect { layer -> layer as int }, + ], + xenium: [:], + ] + } + + if (platform == "XENIUM") { + return [ + dataset_name: "${pairId}_XENIUM", + platform: "XENIUM", + input_path: input.inputDir ?: input.spatialdataPath, + output_path: "spatialdata_out/source_spatialdata.zarr", + persistent_output_path: input.spatialdataPath ?: null, + xenium_spec_path: chooseField(row, ["xenium_spec_path"]) ?: null, + merscope: [:], + xenium: [:], + ] + } + + throw new IllegalArgumentException("Unknown platform: ${platform}") +} + +def segmentMetaForPlatform(row, platform, params) { + if (platform == "MERSCOPE") { + return [ + channels: parseChannels(row.merscope_channels, ["DAPI", "PolyT"]), + image_prefix: chooseField(row, ["merscope_image_prefix"]) ?: null, + z_range: parseRange(row.merscope_z_range, 0, 6), + transform_path: chooseField(row, ["merscope_transform_path"]) ?: null, + xenium_spec_path: null, + min_qv: null, + voxel_layers: intOrDefault( + row.merscope_voxel_layers, + params.default_merscope_voxel_layers + ), + ] + } + + if (platform == "XENIUM") { + return [ + channels: parseChannels(row.xenium_channels, ["DAPI", "18S"]), + image_prefix: null, + z_range: null, + transform_path: null, + xenium_spec_path: chooseField(row, ["xenium_spec_path"]) ?: null, + min_qv: floatOrDefault(row.xenium_min_qv, params.xenium_min_qv), + voxel_layers: intOrDefault( + row.xenium_voxel_layers, + params.default_xenium_voxel_layers + ), + ] + } + + throw new IllegalArgumentException("Unknown platform: ${platform}") +} + +def samplesJsonForPlatforms(pairId, platforms, platformPaths = [:]) { + def samples = platforms.collect { platform -> + def sample = [ + sample_id: "${pairId}_${platform}", + platform: platform, + ] + if (platformPaths.containsKey(platform)) { + sample.zarr_path = platformPaths[platform].toString() + } + return sample + } + return groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(samples)) +} + +def analysisLayerKeys(platform, segmentation) { + if (segmentation in ["proseg", "reseg"]) { + return [ + table_key: "table_MOSAIK_proseg", + shape_key: "MOSAIK_proseg", + ] + } + if (segmentation in ["original", "original_seg"]) { + return [ + table_key: "table_original", + shape_key: platform == "MERSCOPE" + ? "merscope_cell_boundaries" + : "xenium_cell_boundaries", + ] + } + if (segmentation == "proseg_hybrid") { + return [ + table_key: "table_MOSAIK_proseg_hybrid", + shape_key: "MOSAIK_proseg_hybrid", + ] + } + if (segmentation in ["cellpose", "proseg_mask"]) { + return [ + table_key: "table_MOSAIK_cellpose", + shape_key: "MOSAIK_cellpose", + ] + } + if (segmentation == "proseg_geometry_assignment") { + return [ + table_key: "table_MOSAIK_proseg_geometry_assignment", + shape_key: "MOSAIK_proseg", + ] + } + throw new IllegalArgumentException("Unknown analysis segmentation: ${segmentation}") +} + +def normalizeDistanceFromObjectSegmentations(rawValue) { + def rawValues = rawValue instanceof Collection + ? rawValue + : rawValue.toString().split(",") + def aliases = [ + "reseg": "proseg", + "resegmented": "proseg", + "proseg": "proseg", + "original": "original", + "original_seg": "original", + "instrument": "original", + "proseg_mask": "cellpose", + "mask": "cellpose", + "cellpose": "cellpose", + "cellpose_mask": "cellpose", + "proseg_geometry": "proseg_geometry_assignment", + "proseg_geometry_assignment": "proseg_geometry_assignment", + "hybrid": "proseg_hybrid", + "proseg_hybrid": "proseg_hybrid", + ] + def selected = [] + rawValues.each { value -> + def key = value + .toString() + .trim() + .toLowerCase() + .replaceAll(/[^a-z0-9]+/, "_") + .replaceAll(/^_+|_+$/, "") + if (!key) { + return + } + if (!aliases.containsKey(key)) { + throw new IllegalArgumentException( + "Unknown distance_from_object segmentation '${value}'. Valid " + + "values: proseg, original, cellpose, proseg_geometry_assignment" + ) + } + if (!selected.contains(aliases[key])) { + selected << aliases[key] + } + } + return selected ?: ["proseg", "original", "cellpose"] +} + +def normalizeOptionalStringList(rawValue) { + if (rawValue == null) { + return null + } + def values = rawValue instanceof Collection + ? rawValue + : rawValue.toString().split(",") + def cleaned = values + .collect { value -> value.toString().trim() } + .findAll { value -> value.length() > 0 } + .unique() + return cleaned ?: null +} + +def samplesJsonForSegmentation(pairId, platforms, platformPaths, segmentation) { + def samples = platforms.collect { platform -> + def layerKeys = analysisLayerKeys(platform, segmentation) + def sample = [ + sample_id: "${pairId}_${platform}", + platform: platform, + segmentation: segmentation, + table_key: layerKeys.table_key, + shape_key: layerKeys.shape_key, + ] + if (platformPaths.containsKey(platform)) { + sample.zarr_path = platformPaths[platform].toString() + } + return sample + } + return groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(samples)) +} + +def mergeMecrSamplesJson(samplesJsonValues) { + def samplesByInput = [:] + samplesJsonValues.each { samplesJson -> + def samples = new groovy.json.JsonSlurper().parseText(samplesJson.toString()) + samples.each { sample -> + def inputKey = [sample.zarr_path, sample.table_key].join("|") + samplesByInput[inputKey] = sample + } + } + def mergedSamples = samplesByInput.values().toList().sort { sample -> + [sample.sample_id, sample.table_key, sample.zarr_path].join("|") + } + return groovy.json.JsonOutput.prettyPrint( + groovy.json.JsonOutput.toJson(mergedSamples) + ) +} + +def spatialGeneSamplesJson(samplesJson, row) { + def samples = new groovy.json.JsonSlurper().parseText(samplesJson.toString()) + samples.each { sample -> + def platform = sample.platform.toString() + sample.nuclei_shape_key = "cellpose_nuclei" + sample.pial_boundary_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "pial") + ) + sample.wm_boundary_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "wm") + ) + sample.side_boundary_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "side") + ) + sample.exclusion_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "exclusion") + ) + sample.ribbon_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "ribbon") + ) + sample.annotation_path = optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "annotation") + ) + } + return groovy.json.JsonOutput.prettyPrint( + groovy.json.JsonOutput.toJson(samples) + ) +} + +def corticalDepthConfigForPlatform( + row, + pairId, + platform, + analysisSegmentations, + params +) { + def tables = analysisSegmentations.collect { segmentation -> + def layerKeys = analysisLayerKeys(platform, segmentation) + [ + segmentation: segmentation, + table_key: layerKeys.table_key, + shape_key: layerKeys.shape_key, + ] + } + return [ + dataset_name: "${pairId}_${platform}", + platform: platform, + latest_zarr_path: "latest_input.zarr", + output_dir: "compute_cortical_depth_out", + tables: tables, + pial_boundary_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "pial") + ), + wm_boundary_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "wm") + ), + side_boundary_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "side") + ), + exclusion_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "exclusion") + ), + ribbon_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "ribbon") + ), + annotation_path: optionalNormalizedPathString( + corticalDepthAnnotationPath(row, platform, "annotation") + ), + coordinate_unit_um: params.cortical_depth_coordinate_unit_um, + raster_resolution_um: params.cortical_depth_raster_resolution_um, + raster_padding_um: params.cortical_depth_raster_padding_um, + boundary_band_um: params.cortical_depth_boundary_band_um, + boundary_smoothing_window: params.cortical_depth_boundary_smoothing_window, + streamline_spacing_um: params.cortical_depth_streamline_spacing_um, + streamline_step_um: params.cortical_depth_streamline_step_um, + streamline_max_steps: params.cortical_depth_streamline_max_steps, + streamline_resample_points: params.cortical_depth_streamline_resample_points, + side_boundary_distance_um: params.cortical_depth_side_boundary_distance_um, + contour_levels: params.cortical_depth_contour_levels, + write_spatialdata_table: params.cortical_depth_write_spatialdata_table, + ] +} + +def distanceFromObjectConfigForPlatform( + pairId, + platform, + segmentations, + params +) { + def tables = segmentations.collect { segmentation -> + def layerKeys = analysisLayerKeys(platform, segmentation) + [ + segmentation: segmentation, + table_key: layerKeys.table_key, + shape_key: layerKeys.shape_key, + ] + } + return [ + pair_id: pairId, + dataset_name: "${pairId}_${platform}", + platform: platform, + latest_zarr_path: "latest_input.zarr", + output_dir: "distance_from_object_out", + object_annotation_path: "object_annotations.geojson", + tables: tables, + object_types: normalizeOptionalStringList( + params.distance_from_object_object_types + ), + coordinate_unit_um: params.distance_from_object_coordinate_unit_um, + near_distance_um: params.distance_from_object_near_distance_um, + far_distance_um: params.distance_from_object_far_distance_um, + max_distance_um: params.distance_from_object_max_distance_um, + tissue_annotation_column: "cortical_depth_annotation", + included_tissue_annotations: ["grey_matter"], + min_cells_per_pseudobulk: params.distance_from_object_min_cells_per_pseudobulk, + write_spatialdata_table: params.distance_from_object_write_spatialdata_table, + ] +} + +def samplesJsonFromGroupedZarrs(pairId, activePlatforms, platformNames, zarrPaths) { + def names = platformNames.collect { platformName -> platformName.toString() } + def pathsByPlatform = [:] + activePlatforms.each { platform -> + def idx = names.indexOf(platform) + if (idx < 0) { + error "Missing ${platform} latest zarr for ${pairId}" + } + pathsByPlatform[platform] = zarrPaths[idx].toString() + } + return samplesJsonForPlatforms(pairId, activePlatforms, pathsByPlatform) +} + +def normalizeStage(rawValue, paramName) { + def raw = rawValue == null ? "" : rawValue.toString().trim() + if (!raw) { + throw new IllegalArgumentException("Missing required stage value: ${paramName}") + } + def key = raw + .toLowerCase() + .replaceAll(/[^a-z0-9]+/, "_") + .replaceAll(/^_+|_+$/, "") + def aliases = [ + "build": "build_spatialdata", + "build_spatialdata": "build_spatialdata", + "spatialdata": "build_spatialdata", + "spatialdata_build": "build_spatialdata", + "segment": "segment", + "segmentation": "segment", + "segment_nuclei": "segment_nuclei", + "nuclei_segment": "segment_nuclei", + "cellpose_nuclei": "segment_nuclei", + "nuclei_cellpose": "segment_nuclei", + "enrich": "enrich", + "enrichment": "enrich", + "build_viewer_caches": "build_viewer_caches", + "viewer_caches": "build_viewer_caches", + "viewer_cache": "build_viewer_caches", + "pyramids": "build_viewer_caches", + "mask_image_quantification": "mask_image_quantification", + "image_quantification": "mask_image_quantification", + "quantify_images": "mask_image_quantification", + "mask_quantification": "mask_image_quantification", + "compute_cortical_depth": "compute_cortical_depth", + "cortical_depth": "compute_cortical_depth", + "laplace_depth": "compute_cortical_depth", + "distance_from_object": "distance_from_object", + "object_distance": "distance_from_object", + "distance_to_object": "distance_from_object", + "plaque_distance": "distance_from_object", + "qc": "qc", + "align": "align", + "alignment": "align", + "align_qc": "align_qc", + "alignment_qc": "align_qc", + "compare": "compare", + "comparison": "compare", + "visualize": "visualize", + "visualise": "visualize", + "visualization": "visualize", + "visualisation": "visualize", + "spatial_gene_analysis": "spatial_gene_analysis", + "spatial_genes": "spatial_gene_analysis", + "gene_spatial_analysis": "spatial_gene_analysis", + "spatial_autocorrelation": "spatial_gene_analysis", + "mecr": "mecr", + "mutually_exclusive_coexpression": "mecr", + "mutually_exclusive_co_expression": "mecr", + "cluster": "clustering_squidpy", + "clustering": "clustering_squidpy", + "clustering_squidpy": "clustering_squidpy", + "squidpy": "clustering_squidpy", + "cell_type_mapping": "mapmycells", + "celltype_mapping": "mapmycells", + "map_my_cells": "mapmycells", + "mapmycells": "mapmycells", + ] + if (!aliases.containsKey(key)) { + throw new IllegalArgumentException( + "Unknown ${paramName} '${raw}'. Valid stages: " + + "build_spatialdata, segment_nuclei, segment, enrich, " + + "build_viewer_caches, mask_image_quantification, " + + "compute_cortical_depth, distance_from_object, qc, align, align_qc, " + + "compare, visualize, spatial_gene_analysis, mecr, " + + "clustering_squidpy, mapmycells" + ) + } + return aliases[key] +} + +def activeStageOrder( + alignmentEnabled, + pairedMode, + maskImageQuantificationEnabled, + spatialGeneAnalysisEnabled, + corticalDepthEnabled, + distanceFromObjectEnabled, + viewerCacheEnabled, + mecrEnabled +) { + def stages = ["build_spatialdata", "segment_nuclei", "segment", "enrich"] + if (viewerCacheEnabled) { + stages += ["build_viewer_caches"] + } + if (maskImageQuantificationEnabled) { + stages += ["mask_image_quantification"] + } + stages += ["qc"] + if (mecrEnabled) { + stages += ["mecr"] + } + if (pairedMode && alignmentEnabled) { + stages += ["align", "align_qc"] + } + if (pairedMode) { + stages += ["compare"] + } + stages += ["visualize"] + if (spatialGeneAnalysisEnabled) { + stages += ["spatial_gene_analysis"] + } + stages += ["clustering_squidpy"] + // compute_cortical_depth is sequenced AFTER clustering_squidpy so the per-cell + // cluster annotations exist for its depth violin plots. It is a terminal + // opt-in stage; see the default stop_stage handling in rowSampleSettings. + if (corticalDepthEnabled) { + stages += ["compute_cortical_depth"] + } + if (distanceFromObjectEnabled) { + stages += ["distance_from_object"] + } + stages += ["mapmycells"] + return stages +} + +def validateStage( + stage, + stages, + paramName, + alignmentEnabled, + spatialGeneAnalysisEnabled +) { + if (!stages.contains(stage)) { + def hint = "" + if (stage in ["align", "align_qc"] && !alignmentEnabled) { + hint = " Pass --enable_alignment true to use alignment stages." + } + if (stage == "compute_cortical_depth") { + hint = " Pass --cortical_depth_enabled true to use cortical depth." + } + if (stage == "distance_from_object") { + hint = " Pass --distance_from_object_enabled true to use object distance." + } + if (stage == "mecr") { + hint = " Pass --mecr_enabled true to use MECR." + } + if (stage == "spatial_gene_analysis" && !spatialGeneAnalysisEnabled) { + hint = " Pass --spatial_gene_analysis_enabled true to use spatial gene analysis." + } + throw new IllegalArgumentException( + "${paramName} '${stage}' is not active for this run.${hint} " + + "Active stages: ${stages.join(', ')}" + ) + } +} + +def stageInRange(stage, startStage, stopStage, stages) { + def stageIdx = stages.indexOf(stage) + if (stageIdx < 0) { + return false + } + return stageIdx >= stages.indexOf(startStage) && stageIdx <= stages.indexOf(stopStage) +} + +def requireExistingPath(rawPath, label) { + def p = java.nio.file.Paths.get(rawPath.toString()).toAbsolutePath().normalize() + if (!p.toFile().exists()) { + throw new IllegalArgumentException( + "Missing expected ${label}: ${p}\n" + + "This path is required because an upstream stage was skipped. " + + "Either run from an earlier stage or restore the published output." + ) + } + return p +} + +def isBlankPath(rawPath) { + if (rawPath == null) { + return true + } + def value = rawPath.toString().trim() + return !value || value.toLowerCase() == "null" +} + +def normalizedPath(rawPath) { + return java.nio.file.Paths.get(rawPath.toString()).toAbsolutePath().normalize() +} + +def appendPreflightFileCheck(errors, rawPath, label) { + if (isBlankPath(rawPath)) { + errors << "Missing required path parameter for ${label}" + return + } + def path = normalizedPath(rawPath) + if (!path.toFile().exists()) { + errors << "Missing ${label}: ${path}" + } +} + +def findCachedTaxonomyMetadataPath(rawCacheDir) { + if (isBlankPath(rawCacheDir)) { + return null + } + def root = normalizedPath(rawCacheDir) + def taxonomyRoot = root.resolve("abc_whb/metadata/WHB-taxonomy").toFile() + if (!taxonomyRoot.isDirectory()) { + return null + } + def matches = [] + taxonomyRoot.eachDir { versionDir -> + def candidate = versionDir.toPath().resolve("cluster_annotation_term.csv") + if (candidate.toFile().exists()) { + matches << candidate + } + } + return matches ? matches.sort { path -> path.toString() }[-1] : null +} + +def appendClusteringSquidpyPreflightChecks(errors, settings, params) { + def hierarchicalEnabled = boolOrDefault( + params.clustering_squidpy_hierarchical_enabled, + true, + "clustering_squidpy_hierarchical_enabled", + ) + if (!(settings.run_clustering_squidpy && hierarchicalEnabled)) { + return + } + appendPreflightFileCheck( + errors, + params.clustering_squidpy_broad_marker_lookup_path, + "CLUSTERING_SQUIDPY broad marker lookup", + ) + if (isBlankPath(params.clustering_squidpy_broad_taxonomy_metadata_path)) { + def cachedTaxonomy = findCachedTaxonomyMetadataPath( + params.clustering_squidpy_broad_reference_cache_dir, + ) + if (cachedTaxonomy == null) { + errors << ( + "Missing CLUSTERING_SQUIDPY broad taxonomy metadata: " + + "set --clustering_squidpy_broad_taxonomy_metadata_path " + + "or provide a cache containing " + + "abc_whb/metadata/WHB-taxonomy/*/cluster_annotation_term.csv" + ) + } + } else { + appendPreflightFileCheck( + errors, + params.clustering_squidpy_broad_taxonomy_metadata_path, + "CLUSTERING_SQUIDPY broad taxonomy metadata", + ) + } + if (!isBlankPath(params.clustering_squidpy_broad_cluster_membership_path)) { + appendPreflightFileCheck( + errors, + params.clustering_squidpy_broad_cluster_membership_path, + "CLUSTERING_SQUIDPY broad cluster membership metadata", + ) + } +} + +def appendMapMyCellsPreflightChecks(errors, settings, params) { + if (!settings.run_mapmycells) { + return + } + validateMapMyCellsParams(params) + def plotsOnly = params.mapmycells_plots_only == null + ? false + : params.mapmycells_plots_only.toString().trim().toLowerCase() == "true" + if (plotsOnly) { + return + } + def referenceMode = params.mapmycells_reference_mode == null + ? "both" + : params.mapmycells_reference_mode.toString().trim().toLowerCase() + if (referenceMode in ["whole_brain", "both"]) { + appendPreflightFileCheck( + errors, + params.mapmycells_marker_lookup_path, + "MAPMYCELLS whole-brain marker lookup", + ) + appendPreflightFileCheck( + errors, + params.mapmycells_precomputed_stats_path, + "MAPMYCELLS whole-brain precomputed stats", + ) + } +} + +def optionalNormalizedPathString(rawPath) { + return isBlankPath(rawPath) ? null : normalizedPath(rawPath).toString() +} + +def corticalDepthAnnotationPath(row, platform, role) { + def prefix = platform.toString().toLowerCase() + def candidatesByRole = [ + annotation: [ + "${prefix}_cortical_depth_annotation_geojson", + "${prefix}_cortical_depth_annotations_geojson", + "${prefix}_cortical_depth_annotation_path", + "cortical_depth_annotation_geojson", + "cortical_depth_annotations_geojson", + "cortical_depth_annotation_path", + ], + pial: [ + "${prefix}_pial_boundary_geojson", + "${prefix}_pia_boundary_geojson", + "${prefix}_pial_boundary_path", + "pial_boundary_geojson", + "pia_boundary_geojson", + "pial_boundary_path", + ], + wm: [ + "${prefix}_wm_boundary_geojson", + "${prefix}_grey_white_boundary_geojson", + "${prefix}_gray_white_boundary_geojson", + "${prefix}_gm_wm_boundary_geojson", + "${prefix}_wm_boundary_path", + "wm_boundary_geojson", + "grey_white_boundary_geojson", + "gray_white_boundary_geojson", + "gm_wm_boundary_geojson", + "wm_boundary_path", + ], + side: [ + "${prefix}_side_boundary_geojson", + "${prefix}_side_boundaries_geojson", + "${prefix}_tissue_edge_geojson", + "side_boundary_geojson", + "side_boundaries_geojson", + "tissue_edge_geojson", + ], + exclusion: [ + "${prefix}_exclusion_mask_geojson", + "${prefix}_exclusion_masks_geojson", + "${prefix}_cortical_depth_exclusion_geojson", + "exclusion_mask_geojson", + "exclusion_masks_geojson", + "cortical_depth_exclusion_geojson", + ], + ribbon: [ + "${prefix}_cortical_ribbon_geojson", + "${prefix}_ribbon_geojson", + "${prefix}_cortical_ribbon_path", + "cortical_ribbon_geojson", + "ribbon_geojson", + "cortical_ribbon_path", + ], + ] + return chooseField(row, candidatesByRole[role] ?: []) +} + +def distanceFromObjectAnnotationPath(row, platform) { + def prefix = platform.toString().toLowerCase() + return chooseField(row, [ + "${prefix}_distance_object_annotation_geojson", + "${prefix}_distance_from_object_annotation_geojson", + "${prefix}_object_annotation_geojson", + "${prefix}_plaque_annotation_geojson", + "distance_object_annotation_geojson", + "distance_from_object_annotation_geojson", + "object_annotation_geojson", + "plaque_annotation_geojson", + ]) +} + +def appendCorticalDepthPreflightChecks(errors, row, settings, _params) { + if (!settings.run_compute_cortical_depth) { + return + } + settings.active_platforms.each { platform -> + def labelPrefix = "COMPUTE_CORTICAL_DEPTH ${settings.pair_id}:${platform}" + def annotationPath = corticalDepthAnnotationPath(row, platform, "annotation") + if (!isBlankPath(annotationPath)) { + appendPreflightFileCheck( + errors, + annotationPath, + "${labelPrefix} combined annotation GeoJSON", + ) + } else { + appendPreflightFileCheck( + errors, + corticalDepthAnnotationPath(row, platform, "pial"), + "${labelPrefix} pial boundary GeoJSON", + ) + } + ["side", "exclusion", "ribbon"].each { role -> + def optionalPath = corticalDepthAnnotationPath(row, platform, role) + if (!isBlankPath(optionalPath)) { + appendPreflightFileCheck( + errors, + optionalPath, + "${labelPrefix} ${role} annotation GeoJSON", + ) + } + } + } +} + +def appendSpatialGeneTranscriptPreflightChecks(errors, row, settings, _params) { + if ( + !settings.run_spatial_gene_analysis || + !settings.spatial_gene_analysis_transcript_analysis_enabled + ) { + return + } + settings.active_platforms.each { platform -> + def labelPrefix = "SPATIAL_GENE_ANALYSIS ${settings.pair_id}:${platform}" + def annotationPath = corticalDepthAnnotationPath(row, platform, "annotation") + if (!isBlankPath(annotationPath)) { + appendPreflightFileCheck( + errors, + annotationPath, + "${labelPrefix} combined pia/tissue-edge annotation GeoJSON", + ) + } else { + appendPreflightFileCheck( + errors, + corticalDepthAnnotationPath(row, platform, "pial"), + "${labelPrefix} pial boundary GeoJSON", + ) + appendPreflightFileCheck( + errors, + corticalDepthAnnotationPath(row, platform, "side"), + "${labelPrefix} tissue-edge GeoJSON", + ) + } + ["wm", "exclusion", "ribbon"].each { role -> + def optionalPath = corticalDepthAnnotationPath(row, platform, role) + if (!isBlankPath(optionalPath)) { + appendPreflightFileCheck( + errors, + optionalPath, + "${labelPrefix} ${role} annotation GeoJSON", + ) + } + } + } +} + +def appendDistanceFromObjectPreflightChecks(errors, row, settings, _params) { + if (!settings.run_distance_from_object) { + return + } + settings.active_platforms.each { platform -> + appendPreflightFileCheck( + errors, + distanceFromObjectAnnotationPath(row, platform), + "DISTANCE_FROM_OBJECT ${settings.pair_id}:${platform} object GeoJSON", + ) + } +} + +def appendMecrPreflightChecks(errors, settings, params) { + if (!settings.run_mecr) { + return + } + [ + [params.mecr_neurons_h5ad_path, "WHB neuron raw H5AD"], + [params.mecr_nonneurons_h5ad_path, "WHB non-neuron raw H5AD"], + [params.mecr_cell_metadata_path, "WHB cell metadata"], + [params.mecr_taxonomy_metadata_path, "WHB taxonomy metadata"], + [params.mecr_cluster_membership_path, "WHB cluster membership metadata"], + ].each { rawPath, label -> + appendPreflightFileCheck(errors, rawPath, "MECR ${label}") + } +} + +def runPreflightChecks(row, settings, params) { + def errors = [] + appendClusteringSquidpyPreflightChecks(errors, settings, params) + appendMapMyCellsPreflightChecks(errors, settings, params) + appendCorticalDepthPreflightChecks(errors, row, settings, params) + appendSpatialGeneTranscriptPreflightChecks(errors, row, settings, params) + appendDistanceFromObjectPreflightChecks(errors, row, settings, params) + appendMecrPreflightChecks(errors, settings, params) + if (errors) { + throw new IllegalArgumentException( + "Preflight checks failed for sample ${settings.pair_id} " + + "(selected stages: ${settings.selected_stages.join(' -> ')}):\n" + + errors.collect { errorMessage -> " - ${errorMessage}" }.join("\n") + ) + } +} + +def publishedDatasetPath(outdir, pairId, platform, suffix) { + return "${outdir}/${pairId}/${platform.toLowerCase()}/${suffix}" +} + +def publishedPairPath(outdir, pairId, suffix) { + return "${outdir}/${pairId}/${suffix}" +} + +def rowFieldOrDefault(row, fieldName, fallback) { + def value = chooseField(row, [fieldName]) + return value == null ? fallback : value +} + +def boolOrDefault(rawValue, defaultValue, label) { + if (rawValue == null || rawValue.toString().trim().isEmpty()) { + return defaultValue as boolean + } + def key = rawValue + .toString() + .trim() + .toLowerCase() + if (key in ["true", "t", "yes", "y", "1"]) { + return true + } + if (key in ["false", "f", "no", "n", "0"]) { + return false + } + throw new IllegalArgumentException( + "Unknown boolean value for ${label}: '${rawValue}'. " + + "Use true or false." + ) +} + +def rowSampleSettings(row, params) { + def pairId = row.pair_id?.toString()?.trim() + if (!pairId) { + error "Found samplesheet row with missing pair_id: ${row}" + } + + def analysisMode = normalizeAnalysisMode( + rowFieldOrDefault(row, "analysis_mode", params.analysis_mode) + ) + def pairedMode = analysisMode == "paired" + def activePlatforms = activePlatformsForMode(analysisMode) + def analysisSegmentations = normalizeAnalysisSegmentation( + rowFieldOrDefault(row, "analysis_segmentation", params.analysis_segmentation) + ) + def requestedAlignmentEnabled = boolOrDefault( + rowFieldOrDefault(row, "enable_alignment", params.enable_alignment), + false, + "enable_alignment for ${pairId}", + ) + def alignmentEnabled = pairedMode && requestedAlignmentEnabled + def maskImageQuantificationEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "mask_image_quantification_enabled", + params.mask_image_quantification_enabled, + ), + true, + "mask_image_quantification_enabled for ${pairId}", + ) + def viewerCacheEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "viewer_cache_enabled", + params.viewer_cache_enabled, + ), + true, + "viewer_cache_enabled for ${pairId}", + ) + def spatialGeneAnalysisEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "spatial_gene_analysis_enabled", + params.spatial_gene_analysis_enabled, + ), + true, + "spatial_gene_analysis_enabled for ${pairId}", + ) + def spatialGeneTranscriptAnalysisEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "spatial_gene_analysis_transcript_analysis_enabled", + params.spatial_gene_analysis_transcript_analysis_enabled, + ), + true, + "spatial_gene_analysis_transcript_analysis_enabled for ${pairId}", + ) + def corticalDepthEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "cortical_depth_enabled", + params.cortical_depth_enabled, + ), + false, + "cortical_depth_enabled for ${pairId}", + ) + def distanceFromObjectEnabled = boolOrDefault( + rowFieldOrDefault( + row, + "distance_from_object_enabled", + params.distance_from_object_enabled, + ), + false, + "distance_from_object_enabled for ${pairId}", + ) + def mecrEnabled = boolOrDefault( + rowFieldOrDefault(row, "mecr_enabled", params.mecr_enabled), + true, + "mecr_enabled for ${pairId}", + ) + def distanceFromObjectSegmentations = normalizeDistanceFromObjectSegmentations( + rowFieldOrDefault( + row, + "distance_from_object_segmentations", + params.distance_from_object_segmentations, + ) + ) + + def rowOnlyStageRaw = chooseField(row, ["only_stage"]) + def rowStartStageRaw = chooseField(row, ["start_stage"]) + def rowStopStageRaw = chooseField(row, ["stop_stage"]) + def hasRowStageRange = rowStartStageRaw != null || rowStopStageRaw != null + def globalOnlyStageRaw = params.only_stage == null + ? null + : params.only_stage.toString().trim() + if (!globalOnlyStageRaw) { + globalOnlyStageRaw = null + } + + def onlyStageRaw = rowOnlyStageRaw ?: (hasRowStageRange ? null : globalOnlyStageRaw) + def startStageRaw = onlyStageRaw ?: (rowStartStageRaw ?: params.start_stage) + def stopStageRaw = onlyStageRaw ?: (rowStopStageRaw ?: params.stop_stage) + def startParamName = onlyStageRaw ? "only_stage" : "start_stage" + def stopParamName = onlyStageRaw ? "only_stage" : "stop_stage" + if (rowOnlyStageRaw) { + startParamName = "samplesheet only_stage for ${pairId}" + stopParamName = "samplesheet only_stage for ${pairId}" + } else if (hasRowStageRange) { + startParamName = "samplesheet start_stage for ${pairId}" + stopParamName = "samplesheet stop_stage for ${pairId}" + } + + def stageOrder = activeStageOrder( + alignmentEnabled, + pairedMode, + maskImageQuantificationEnabled, + spatialGeneAnalysisEnabled, + corticalDepthEnabled, + distanceFromObjectEnabled, + viewerCacheEnabled, + mecrEnabled, + ) + def startStage = normalizeStage(startStageRaw, startParamName) + def stopStage = normalizeStage(stopStageRaw, stopParamName) + validateStage( + startStage, + stageOrder, + startParamName, + alignmentEnabled, + spatialGeneAnalysisEnabled, + ) + validateStage( + stopStage, + stageOrder, + stopParamName, + alignmentEnabled, + spatialGeneAnalysisEnabled, + ) + // Opt-in terminal analyses extend the historical clustering stop default. + // Explicit only_stage selections and non-default stop stages remain exact. + if (!onlyStageRaw && stopStage == "clustering_squidpy") { + if (distanceFromObjectEnabled && stageOrder.contains("distance_from_object")) { + stopStage = "distance_from_object" + } else if (corticalDepthEnabled && stageOrder.contains("compute_cortical_depth")) { + stopStage = "compute_cortical_depth" + } + } + if (stageOrder.indexOf(startStage) > stageOrder.indexOf(stopStage)) { + error( + "Samplesheet row ${pairId} has start_stage '${startStage}' " + + "after stop_stage '${stopStage}'." + ) + } + + def runBuild = stageInRange("build_spatialdata", startStage, stopStage, stageOrder) + def runSegmentNuclei = stageInRange( + "segment_nuclei", + startStage, + stopStage, + stageOrder, + ) + def runSegment = stageInRange("segment", startStage, stopStage, stageOrder) + def runEnrich = stageInRange("enrich", startStage, stopStage, stageOrder) + def runBuildViewerCaches = stageInRange( + "build_viewer_caches", + startStage, + stopStage, + stageOrder, + ) + def runMaskImageQuantification = stageInRange( + "mask_image_quantification", + startStage, + stopStage, + stageOrder, + ) + def runComputeCorticalDepth = stageInRange( + "compute_cortical_depth", + startStage, + stopStage, + stageOrder, + ) + def runDistanceFromObject = stageInRange( + "distance_from_object", + startStage, + stopStage, + stageOrder, + ) + def runQc = stageInRange("qc", startStage, stopStage, stageOrder) + def runMecr = stageInRange("mecr", startStage, stopStage, stageOrder) + def runAlign = stageInRange("align", startStage, stopStage, stageOrder) + def runAlignQc = stageInRange("align_qc", startStage, stopStage, stageOrder) + def runCompare = stageInRange("compare", startStage, stopStage, stageOrder) + def runVisualize = stageInRange("visualize", startStage, stopStage, stageOrder) + def runSpatialGeneAnalysis = stageInRange( + "spatial_gene_analysis", + startStage, + stopStage, + stageOrder, + ) + def runClusteringSquidpy = stageInRange( + "clustering_squidpy", + startStage, + stopStage, + stageOrder, + ) + def runMapMyCells = stageInRange("mapmycells", startStage, stopStage, stageOrder) + def needAnalysisZarrs = + runQc || + runAlign || + runAlignQc || + runMecr || + runCompare || + runVisualize || + runSpatialGeneAnalysis || + runClusteringSquidpy + def needAlignmentDownstream = + pairedMode && + alignmentEnabled && + ( + runCompare || + runVisualize || + runSpatialGeneAnalysis || + runClusteringSquidpy + ) + def needAlignmentResults = + pairedMode && + alignmentEnabled && + (runAlign || runAlignQc || needAlignmentDownstream) + + return [ + pair_id: pairId, + analysis_mode: analysisMode, + enable_alignment: alignmentEnabled, + active_platforms: activePlatforms, + paired_mode: pairedMode, + analysis_segmentations: analysisSegmentations, + distance_from_object_segmentations: distanceFromObjectSegmentations, + spatial_gene_analysis_enabled: spatialGeneAnalysisEnabled, + spatial_gene_analysis_transcript_analysis_enabled: ( + spatialGeneAnalysisEnabled && spatialGeneTranscriptAnalysisEnabled + ), + stage_order: stageOrder, + start_stage: startStage, + stop_stage: stopStage, + selected_stages: stageOrder.findAll { stage -> + stageInRange(stage, startStage, stopStage, stageOrder) + }, + run_build: runBuild, + run_segment_nuclei: runSegmentNuclei, + run_segment: runSegment, + run_enrich: runEnrich, + run_build_viewer_caches: runBuildViewerCaches, + run_mask_image_quantification: runMaskImageQuantification, + run_compute_cortical_depth: runComputeCorticalDepth, + run_distance_from_object: runDistanceFromObject, + run_qc: runQc, + run_mecr: runMecr, + run_align: runAlign, + run_align_qc: runAlignQc, + run_compare: runCompare, + run_visualize: runVisualize, + run_spatial_gene_analysis: runSpatialGeneAnalysis, + run_clustering_squidpy: runClusteringSquidpy, + run_mapmycells: runMapMyCells, + need_build_results: runSegmentNuclei || runSegment || runEnrich, + need_enriched_zarrs: ( + runBuildViewerCaches || + runMaskImageQuantification || + runComputeCorticalDepth || + runDistanceFromObject || + runQc || + needAnalysisZarrs + ), + need_quantified_zarrs: runMaskImageQuantification, + need_analysis_zarrs: needAnalysisZarrs, + need_alignment_results: needAlignmentResults, + need_alignment_downstream: needAlignmentDownstream, + ] +} + +def validateMapMyCellsParams(params) { + def mapMyCellsReferenceMode = params.mapmycells_reference_mode == null + ? "both" + : params.mapmycells_reference_mode.toString().trim().toLowerCase() + def mapMyCellsPlotsOnly = params.mapmycells_plots_only == null + ? false + : params.mapmycells_plots_only.toString().trim().toLowerCase() == "true" + if (!(mapMyCellsReferenceMode in ["whole_brain", "region", "both"])) { + throw new IllegalArgumentException( + "Invalid MAPMYCELLS --mapmycells_reference_mode " + + "'${params.mapmycells_reference_mode}'. Valid values: " + + "whole_brain, region, both" + ) + } + if (!mapMyCellsPlotsOnly && + mapMyCellsReferenceMode in ["whole_brain", "both"] && + !params.mapmycells_marker_lookup_path) { + throw new IllegalArgumentException( + "Missing required parameter for MAPMYCELLS: --mapmycells_marker_lookup_path" + ) + } + if (!mapMyCellsPlotsOnly && + mapMyCellsReferenceMode in ["whole_brain", "both"] && + !params.mapmycells_precomputed_stats_path) { + throw new IllegalArgumentException( + "Missing required parameter for MAPMYCELLS: --mapmycells_precomputed_stats_path" + ) + } + if (!mapMyCellsPlotsOnly && + mapMyCellsReferenceMode in ["region", "both"] && + !params.mapmycells_region_labels) { + throw new IllegalArgumentException( + "Missing required parameter for MAPMYCELLS region mode: " + + "--mapmycells_region_labels" + ) + } +} + +workflow { + if (!params.samplesheet) { + error "Missing required parameter: --samplesheet" + } + + samplesheet_ch = channel + .fromPath(params.samplesheet, checkIfExists: true) + .splitCsv(header: true, sep: ",", quote: '"', strip: true) + + sample_rows_raw_ch = samplesheet_ch.map { row -> + def settings = rowSampleSettings(row, params) + log.info( + "Sample ${settings.pair_id}: analysis_mode=${settings.analysis_mode}; " + + "enable_alignment=${settings.enable_alignment}; " + + "spatial_gene_analysis=${settings.run_spatial_gene_analysis}; " + + "spatial_gene_transcripts=" + + "${settings.spatial_gene_analysis_transcript_analysis_enabled}; " + + "cortical_depth=${settings.run_compute_cortical_depth}; " + + "distance_from_object=${settings.run_distance_from_object}; " + + "active platforms=${settings.active_platforms.join(', ')}; " + + "analysis segmentations=${settings.analysis_segmentations.join(', ')}; " + + "distance segmentations=" + + "${settings.distance_from_object_segmentations.join(', ')}; " + + "selected stages=${settings.selected_stages.join(' -> ')}" + ) + tuple(settings.pair_id, row, settings) + } + + preflight_done_ch = sample_rows_raw_ch + .map { _pairId, row, settings -> + runPreflightChecks(row, settings, params) + true + } + .collect() + .map { true } + + sample_rows_ch = sample_rows_raw_ch + .combine(preflight_done_ch) + .map { pairId, row, settings, _doneFlag -> + tuple(pairId, row, settings) + } + + build_inputs_ch = sample_rows_ch.flatMap { pairId, row, settings -> + if (!settings.run_build) { + [] + } else { + settings.active_platforms.collect { platform -> + def key = "${pairId}|${platform}" + def buildConfig = buildConfigForPlatform(row, pairId, platform) + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(buildConfig)), + ) + } + } + } + + build_task_results_ch = BUILD_SPATIALDATA(build_inputs_ch) + + build_published_results_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_build_results && !settings.run_build)) { + [] + } else { + settings.active_platforms.collect { platform -> + def key = "${pairId}|${platform}" + def sourceSpatialdata = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "spatialdata/spatialdata_out/source_spatialdata.zarr", + ), + "BUILD_SPATIALDATA output for ${pairId}:${platform}", + ) + tuple(key, pairId, platform, sourceSpatialdata) + } + } + } + + build_results_ch = build_task_results_ch.mix(build_published_results_ch) + + segment_meta_ch = sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_segment_nuclei || settings.run_segment)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple( + "${pairId}|${platform}", + segmentMetaForPlatform(row, platform, params), + settings.run_segment_nuclei, + settings.run_segment, + ) + } + } + } + + segmentation_config_inputs_ch = build_results_ch + .join(segment_meta_ch) + .map { + key, pairId, platform, sourceSpatialdata, meta, + runSegmentNuclei, runSegment -> + def persistentLatestZarrPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "latest/latest_spatialdata.zarr", + ) + ).toAbsolutePath().toString() + def persistentMaskPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_masks_tiled.npy", + ) + ).toAbsolutePath().toString() + def persistentCellprobPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_cellprobs_tiled.npy", + ) + ).toAbsolutePath().toString() + def persistentTranscriptsPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/transcripts_for_proseg.csv", + ) + ).toAbsolutePath().toString() + def persistentStitchingStatsPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_stitching_stats.json", + ) + ).toAbsolutePath().toString() + def persistentNucleiMaskPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_nuclei_masks_tiled.npy", + ) + ).toAbsolutePath().toString() + def persistentNucleiStitchingStatsPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_nuclei_stitching_stats.json", + ) + ).toAbsolutePath().toString() + def baseConfig = [ + cellpose: [ + model_type: params.cellpose_model_type, + gpu: params.cellpose_gpu, + diameter: params.cellpose_diameter, + flow_threshold: params.cellpose_flow_threshold, + cellprob_threshold: params.cellpose_cellprob, + tile_overlap: params.cellpose_tile_overlap, + bsize: params.cellpose_bsize, + factor_rescale: 1.0, + use_bfloat16: params.cellpose_use_bfloat16, + ], + nuclei_cellpose: [ + model_type: "nuclei", + gpu: params.cellpose_gpu, + diameter: params.cellpose_diameter, + flow_threshold: params.cellpose_flow_threshold, + cellprob_threshold: params.cellpose_cellprob, + tile_overlap: params.cellpose_tile_overlap, + bsize: params.cellpose_bsize, + factor_rescale: 1.0, + use_bfloat16: params.cellpose_use_bfloat16, + ], + tiling: [ + tile_size_candidates: params.cellpose_tile_size_candidates, + stitch_overlap_px: params.cellpose_stitch_overlap_px, + min_tile_size: params.cellpose_min_tile_size, + status_every_tiles: params.cellpose_stitch_status_every_tiles, + filter_per_tile: params.cellpose_filter_per_tile, + duplicate_iou_threshold: params.cellpose_duplicate_iou_threshold, + duplicate_overlap_fraction: params.cellpose_duplicate_overlap_fraction, + min_remaining_fraction: params.cellpose_min_remaining_fraction, + edge_touch_policy: params.cellpose_edge_touch_policy, + write_stitching_stats: params.cellpose_write_stitching_stats, + ], + mask_filter: [ + final_min_area_um2: params.cellpose_final_min_area_um2, + final_max_area_um2: params.cellpose_final_max_area_um2, + final_filter_chunk_mb: params.cellpose_final_filter_chunk_mb, + ], + nuclei_mask_filter: [ + final_min_area_um2: params.cellpose_final_min_area_um2, + final_max_area_um2: params.cellpose_final_max_area_um2, + final_filter_chunk_mb: params.cellpose_final_filter_chunk_mb, + ], + proseg: [ + binary_path: params.proseg_binary, + samples: params.proseg_samples, + voxel_size: params.proseg_voxel_size, + burnin_voxel_size: params.proseg_burnin_voxel_size, + nuclear_reassignment_prob: params.proseg_nuclear_reassignment_prob, + diffusion_probability: params.proseg_diffusion_probability, + cell_compactness: params.proseg_cell_compactness, + num_threads: params.proseg_num_threads, + voxel_layers: 2, + ], + proseg_hybrid: [ + enabled: params.proseg_hybrid_enabled, + min_transcripts: params.proseg_hybrid_min_transcripts, + outlier_neighbors: params.proseg_hybrid_outlier_neighbors, + outlier_mad_multiplier: params.proseg_hybrid_outlier_mad_multiplier, + minimum_external_group: params.proseg_hybrid_minimum_external_group, + chain_radius_scale: params.proseg_hybrid_chain_radius_scale, + near_surface_radius_fraction: params.proseg_hybrid_near_surface_radius_fraction, + maximum_expansion_radius_fraction: params.proseg_hybrid_maximum_expansion_radius_fraction, + attachment_arc_width_scale: params.proseg_hybrid_attachment_arc_width_scale, + rounding_radius_fraction: params.proseg_hybrid_rounding_radius_fraction, + smoothing_radius_um: params.proseg_hybrid_smoothing_radius_um, + outward_rounding_um: params.proseg_hybrid_outward_rounding_um, + smoothing_quad_segs: params.proseg_hybrid_smoothing_quad_segs, + containment_tolerance_um: params.proseg_hybrid_containment_tolerance_um, + ], + memory: [ + max_system_ram_gb: params.max_ram_gb, + memory_warn_gb: params.warn_ram_gb, + transcript_chunk_rows: params.transcript_chunk_rows, + ], + ] + + def segmentConfig = baseConfig + [ + dataset: [ + name: "${pairId}_${platform}", + platform: platform, + data_path: sourceSpatialdata.toString(), + channels: meta.channels, + output_dir: "segment_out", + persistent_latest_zarr_path: persistentLatestZarrPath, + persistent_mask_path: persistentMaskPath, + persistent_cellpose_cellprob_path: persistentCellprobPath, + persistent_transcripts_path: persistentTranscriptsPath, + persistent_cellpose_stitching_stats_path: persistentStitchingStatsPath, + persistent_nuclei_mask_path: persistentNucleiMaskPath, + persistent_nuclei_stitching_stats_path: persistentNucleiStitchingStatsPath, + image_prefix: meta.image_prefix, + z_range: meta.z_range, + transform_path: meta.transform_path, + xenium_spec_path: meta.xenium_spec_path, + min_qv: meta.min_qv, + proseg_overrides: [voxel_layers: meta.voxel_layers], + ], + ] + + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(segmentConfig)), + runSegmentNuclei, + runSegment, + ) + } + + nuclei_segment_inputs_ch = segmentation_config_inputs_ch + .filter { + _key, _pairId, _platform, _segmentConfigJson, + runSegmentNuclei, _runSegment -> + runSegmentNuclei + } + .map { + key, pairId, platform, segmentConfigJson, + _runSegmentNuclei, _runSegment -> + tuple(key, pairId, platform, segmentConfigJson) + } + + cellpose_segment_inputs_ch = segmentation_config_inputs_ch + .filter { + _key, _pairId, _platform, _segmentConfigJson, + _runSegmentNuclei, runSegment -> + runSegment + } + .map { + key, pairId, platform, segmentConfigJson, + _runSegmentNuclei, _runSegment -> + tuple(key, pairId, platform, segmentConfigJson) + } + + nuclei_task_results_ch = CELLPOSE_NUCLEI_SEGMENT(nuclei_segment_inputs_ch) + + nuclei_published_results_ch = segmentation_config_inputs_ch + .filter { + _key, _pairId, _platform, _segmentConfigJson, + runSegmentNuclei, runSegment -> + runSegment && !runSegmentNuclei + } + .map { + key, pairId, platform, segmentConfigJson, + _runSegmentNuclei, _runSegment -> + def nucleiMask = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_nuclei_masks_tiled.npy", + ), + "Cellpose nuclei mask for ${pairId}:${platform}", + ) + def nucleiStats = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_nuclei_stitching_stats.json", + ), + "Cellpose nuclei stitching stats for ${pairId}:${platform}", + ) + tuple( + key, + pairId, + platform, + segmentConfigJson, + nucleiMask, + nucleiStats, + ) + } + + nuclei_results_ch = nuclei_task_results_ch.mix(nuclei_published_results_ch) + proseg_trigger_ch = cellpose_segment_inputs_ch.map { true }.take(1) + proseg_path_ch = ENSURE_PROSEG(proseg_trigger_ch) + segment_task_results_ch = SEGMENT( + cellpose_segment_inputs_ch, + nuclei_results_ch, + proseg_path_ch, + ) + + segment_published_results_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.run_enrich && !settings.run_segment)) { + [] + } else { + settings.active_platforms.collect { platform -> + def key = "${pairId}|${platform}" + def latestZarr = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/segment_out/proseg_base_latest.zarr", + ), + "SEGMENT latest zarr for ${pairId}:${platform}", + ) + def maskPath = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_masks_tiled.npy", + ), + "SEGMENT Cellpose mask for ${pairId}:${platform}", + ) + def transcriptsCsv = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/transcripts_for_proseg.csv", + ), + "SEGMENT transcript CSV for ${pairId}:${platform}", + ) + def nucleiMaskPath = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_nuclei_masks_tiled.npy", + ), + "SEGMENT Cellpose nuclei mask for ${pairId}:${platform}", + ) + tuple( + key, + pairId, + platform, + latestZarr, + maskPath, + transcriptsCsv, + nucleiMaskPath, + ) + } + } + } + + segment_results_ch = segment_task_results_ch.mix(segment_published_results_ch) + + metadata_ch = build_results_ch.map { key, pairId, platform, sourceSpatialdata -> + tuple(key, pairId, platform, sourceSpatialdata.toString()) + } + + enrich_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!settings.run_enrich) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", true) + } + } + } + + enrich_inputs_ch = segment_results_ch + .join(metadata_ch) + .join(enrich_gate_ch) + .map { + key, pairId, platform, latestZarr, maskPath, _transcriptsCsv, + nucleiMaskPath, + pairMeta, platformMeta, originalDataPath, _runFlag -> + if (pairId != pairMeta || platform != platformMeta) { + error( + "Internal channel mismatch for key=${key}: " + + "${pairId}/${platform} vs ${pairMeta}/${platformMeta}" + ) + } + + def persistentLatestZarrPath = file( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "latest/latest_spatialdata.zarr", + ) + ).toAbsolutePath().toString() + + def enrichConfig = [ + dataset_name: "${pairId}_${platform}", + platform: platform, + latest_zarr_path: "latest_input.zarr", + mask_path: "enrich_input_mask.npy", + nuclei_mask_path: "enrich_input_nuclei_mask.npy", + original_data_path: originalDataPath, + output_dir: "enrich_out", + persistent_output_path: persistentLatestZarrPath, + transform_path: null, + ] + + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(enrichConfig)), + latestZarr, + maskPath, + nucleiMaskPath, + ) + } + + enrich_task_results_ch = ENRICH(enrich_inputs_ch) + + enrich_published_results_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_enriched_zarrs && !settings.run_enrich)) { + [] + } else { + settings.active_platforms.collect { platform -> + def key = "${pairId}|${platform}" + def latestZarr = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "latest/latest_spatialdata.zarr", + ), + "ENRICH latest zarr for ${pairId}:${platform}", + ) + // The published latest zarr is the only ENRICH artifact consumed + // downstream. Older completed runs may not retain enrich_out, so + // use the required latest zarr as the unused provenance placeholder + // instead of rejecting an otherwise valid downstream stage start. + def enrichOut = latestZarr + tuple(key, pairId, platform, latestZarr, enrichOut) + } + } + } + + enrich_results_ch = enrich_task_results_ch.mix(enrich_published_results_ch) + + enriched_zarrs_ch = enrich_results_ch.map { + key, pairId, platform, enrichedLatestZarr, _enrichOutDir -> + tuple( + key, + pairId, + platform, + java.nio.file.Paths.get(enrichedLatestZarr.toString()).toRealPath().toString(), + ) + } + + // Pre-build the napari viewer's derived caches (label masks + label/outline + // pyramids + image pyramid) into the enriched latest zarr. This is a separate + // writer of the shared store, so it is serialized into the post-enrich chain: + // the stages that consume the enriched zarr read `post_enrich_zarrs_ch` below, + // which is the viewer-cache output when the stage runs and the enriched zarr + // passthrough otherwise. + viewer_cache_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!settings.run_build_viewer_caches) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", true) + } + } + } + + viewer_cache_transform_ch = sample_rows_ch.flatMap { pairId, row, settings -> + if (!settings.run_build_viewer_caches) { + [] + } else { + settings.active_platforms.collect { platform -> + def transformPath = platform == "MERSCOPE" + ? chooseField(row, ["merscope_transform_path"]) + : chooseField(row, ["xenium_spec_path"]) + tuple("${pairId}|${platform}", transformPath) + } + } + } + + viewer_cache_inputs_ch = enriched_zarrs_ch + .join(viewer_cache_gate_ch) + .join(viewer_cache_transform_ch) + .map { key, pairId, platform, enrichedLatestZarr, _runFlag, transformPath -> + def resolvedTransform = (transformPath && transformPath.toString().trim()) + ? file(transformPath).toAbsolutePath().toString() + : null + def viewerCacheConfig = [ + dataset_name: "${pairId}_${platform}", + platform: platform, + latest_zarr_path: "latest_input.zarr", + original_data_path: "latest_input.zarr", + output_dir: "viewer_cache_out", + transform_path: resolvedTransform, + downsample: params.viewer_cache_downsample, + label_chunk_size: params.viewer_cache_label_chunk_size, + contour_width: params.viewer_cache_contour_width, + min_size: params.viewer_cache_min_size, + build_image_pyramid: params.viewer_cache_build_image_pyramid, + ] + + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(viewerCacheConfig)), + enrichedLatestZarr, + ) + } + + viewer_cache_results_ch = VIEWER_CACHE(viewer_cache_inputs_ch) + + viewer_cached_zarrs_ch = viewer_cache_results_ch.map { + key, pairId, platform, cachedLatestZarr, _viewerCacheOutDir -> + tuple( + key, + pairId, + platform, + java.nio.file.Paths.get(cachedLatestZarr.toString()).toRealPath().toString(), + ) + } + + viewer_cache_passthrough_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_enriched_zarrs && !settings.run_build_viewer_caches)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", true) + } + } + } + + viewer_cache_passthrough_ch = enriched_zarrs_ch + .join(viewer_cache_passthrough_gate_ch) + .map { key, pairId, platform, enrichedLatestZarr, _runFlag -> + tuple(key, pairId, platform, enrichedLatestZarr) + } + + post_enrich_zarrs_ch = viewer_cached_zarrs_ch.mix(viewer_cache_passthrough_ch) + + mask_image_quantification_gate_ch = sample_rows_ch.flatMap { + pairId, _row, settings -> + if (!settings.run_mask_image_quantification) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", true) + } + } + } + + mask_image_quantification_task_masks_ch = segment_task_results_ch.map { + key, _pairId, _platform, _latestZarr, maskPath, _transcriptsCsv, + _nucleiMaskPath -> + tuple(key, maskPath) + } + + mask_image_quantification_published_masks_ch = sample_rows_ch.flatMap { + pairId, _row, settings -> + if (!( + settings.run_mask_image_quantification && + !settings.run_segment + )) { + [] + } else { + settings.active_platforms.collect { platform -> + def key = "${pairId}|${platform}" + def maskPath = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + platform, + "segmentation/cellpose_masks_tiled.npy", + ), + "SEGMENT Cellpose mask for ${pairId}:${platform}", + ) + tuple(key, maskPath) + } + } + } + + mask_image_quantification_masks_ch = + mask_image_quantification_task_masks_ch.mix( + mask_image_quantification_published_masks_ch + ) + + mask_image_quantification_inputs_ch = post_enrich_zarrs_ch + .join(mask_image_quantification_masks_ch) + .join(mask_image_quantification_gate_ch) + .map { key, pairId, platform, enrichedLatestZarr, maskPath, _runFlag -> + def quantConfig = [ + dataset_name: "${pairId}_${platform}", + platform: platform, + latest_zarr_path: "latest_input.zarr", + mask_path: "mask_image_quantification_input_mask.npy", + output_dir: "mask_image_quantification_out", + ] + + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(quantConfig)), + enrichedLatestZarr, + maskPath, + ) + } + + mask_image_quantification_results_ch = MASK_IMAGE_QUANTIFICATION( + mask_image_quantification_inputs_ch + ) + + quantified_zarrs_ch = mask_image_quantification_results_ch.map { + key, pairId, platform, quantifiedLatestZarr, _quantOutDir -> + tuple( + key, + pairId, + platform, + java.nio.file.Paths.get(quantifiedLatestZarr.toString()).toRealPath().toString(), + ) + } + + enriched_downstream_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_enriched_zarrs && !settings.need_quantified_zarrs)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", true) + } + } + } + + enriched_downstream_zarrs_ch = post_enrich_zarrs_ch + .join(enriched_downstream_gate_ch) + .map { key, pairId, platform, enrichedLatestZarr, _runFlag -> + tuple(key, pairId, platform, enrichedLatestZarr) + } + + downstream_zarrs_ch = enriched_downstream_zarrs_ch.mix(quantified_zarrs_ch) + + // Cortical depth now runs after clustering_squidpy (wired below) so that the + // per-cell broad_class / subcluster_label annotations exist and the depth + // violin plots can be produced. The analysis branch therefore reads the + // enriched/quantified zarrs directly and does not wait on cortical depth, + // which also breaks what would otherwise be a cycle: + // clustering -> cortical depth -> clustering. + analysis_ready_zarrs_ch = downstream_zarrs_ch + + analysis_layer_validation_gate_ch = sample_rows_ch.flatMap { + pairId, _row, settings -> + if (!settings.need_analysis_zarrs) { + [] + } else { + settings.active_platforms.collectMany { platform -> + settings.analysis_segmentations.collect { segmentation -> + tuple("${pairId}|${platform}", segmentation, settings) + } + } + } + } + + analysis_layer_validation_inputs_ch = analysis_ready_zarrs_ch + .join(analysis_layer_validation_gate_ch) + .map { + key, pairId, platform, enrichedLatestZarr, segmentation, settings -> + def layerKeys = analysisLayerKeys(platform, segmentation) + tuple( + "${key}|${segmentation}", + pairId, + platform, + segmentation, + enrichedLatestZarr, + layerKeys.table_key, + layerKeys.shape_key, + settings, + ) + } + + analysis_layer_validation_results_ch = VALIDATE_ANALYSIS_LAYER( + analysis_layer_validation_inputs_ch + ) + + qc_inputs_ch = analysis_layer_validation_results_ch + .filter { + _key, _pairId, _platform, _segmentation, _latestZarr, + _tableKey, _shapeKey, settings, _validationJson -> + settings.run_qc + } + .map { + key, pairId, platform, segmentation, latestZarr, + tableKey, shapeKey, _settings, _validationJson -> + tuple( + key, + pairId, + platform, + segmentation, + latestZarr, + tableKey, + shapeKey, + ) + } + + qc_results_ch = QC(qc_inputs_ch) + + analysis_branch_settings_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!settings.need_analysis_zarrs) { + [] + } else { + settings.active_platforms.collectMany { platform -> + settings.analysis_segmentations.collect { segmentation -> + tuple("${pairId}|${platform}|${segmentation}", settings) + } + } + } + } + + analysis_from_qc_ch = qc_results_ch + .map { + _key, pairId, platform, segmentation, enrichedLatestZarr, + _qcOutDir, tableKey, shapeKey -> + tuple( + "${pairId}|${platform}|${segmentation}", + pairId, + segmentation, + platform, + java.nio.file.Paths.get(enrichedLatestZarr.toString()).toRealPath().toString(), + tableKey, + shapeKey, + ) + } + .join(analysis_branch_settings_ch) + .map { + _branchKey, pairId, segmentation, platform, zarrPath, + tableKey, shapeKey, settings -> + tuple( + pairId, + segmentation, + platform, + zarrPath, + tableKey, + shapeKey, + settings, + ) + } + + analysis_without_qc_ch = analysis_layer_validation_results_ch + .filter { + _key, _pairId, _platform, _segmentation, _latestZarr, + _tableKey, _shapeKey, settings, _validationJson -> + !settings.run_qc + } + .map { + _key, pairId, platform, segmentation, latestZarr, + tableKey, shapeKey, settings, _validationJson -> + tuple( + pairId, + segmentation, + platform, + latestZarr, + tableKey, + shapeKey, + settings, + ) + } + + analysis_dataset_zarrs_ch = analysis_from_qc_ch.mix(analysis_without_qc_ch) + + mecr_samples_ch = analysis_dataset_zarrs_ch + .filter { + _pairId, _segmentation, _platform, _zarrPath, _tableKey, _shapeKey, + settings -> + settings.run_mecr + } + .map { + pairId, segmentation, platform, zarrPath, _tableKey, _shapeKey, + settings -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + settings, + platform, + zarrPath, + ) + } + .groupTuple() + .map { + _branchKey, pairIds, segmentations, settingsList, platformNames, + zarrPaths -> + def pairId = pairIds[0].toString() + def segmentation = segmentations[0].toString() + def settings = settingsList[0] + def pathsByPlatform = [:] + platformNames.eachWithIndex { platform, idx -> + pathsByPlatform[platform.toString()] = zarrPaths[idx].toString() + } + tuple( + pairId, + segmentation, + samplesJsonForSegmentation( + pairId, + settings.active_platforms, + pathsByPlatform, + segmentation, + ), + ) + } + + mecr_reference_inputs_ch = mecr_samples_ch + .map { _pairId, _segmentation, samplesJson -> samplesJson } + .collect() + .map { samplesJsonValues -> mergeMecrSamplesJson(samplesJsonValues) } + mecr_reference_results_ch = MECR_REFERENCE(mecr_reference_inputs_ch) + mecr_inputs_ch = mecr_samples_ch + .combine(mecr_reference_results_ch) + .map { pairId, segmentation, samplesJson, referenceOut -> + tuple(pairId, segmentation, samplesJson, referenceOut) + } + mecr_results_ch = MECR(mecr_inputs_ch) + + merscope_zarr_ch = analysis_ready_zarrs_ch + .filter { _key, _pairId, platform, _zarrPath -> platform == "MERSCOPE" } + .map { _key, pairId, _platform, zarrPath -> tuple(pairId, zarrPath) } + + xenium_zarr_ch = analysis_ready_zarrs_ch + .filter { _key, _pairId, platform, _zarrPath -> platform == "XENIUM" } + .map { _key, pairId, _platform, zarrPath -> tuple(pairId, zarrPath) } + + paired_need_zarrs_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.paired_mode && settings.need_analysis_zarrs)) { + [] + } else { + [tuple(pairId, settings)] + } + } + + paired_zarrs_ch = merscope_zarr_ch + .join(xenium_zarr_ch) + .join(paired_need_zarrs_ch) + .map { pairId, merscopePath, xeniumPath, settings -> + tuple(pairId, merscopePath, xeniumPath, settings) + } + + merscope_analysis_ch = analysis_dataset_zarrs_ch + .filter { + _pairId, _segmentation, platform, _zarrPath, _tableKey, _shapeKey, settings -> + settings.paired_mode && platform == "MERSCOPE" + } + .map { + pairId, segmentation, _platform, zarrPath, tableKey, shapeKey, settings -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + zarrPath, + tableKey, + shapeKey, + settings, + ) + } + + xenium_analysis_ch = analysis_dataset_zarrs_ch + .filter { + _pairId, _segmentation, platform, _zarrPath, _tableKey, _shapeKey, settings -> + settings.paired_mode && platform == "XENIUM" + } + .map { + pairId, segmentation, _platform, zarrPath, tableKey, shapeKey, _settings -> + tuple( + "${pairId}|${segmentation}", + zarrPath, + tableKey, + shapeKey, + ) + } + + paired_analysis_zarrs_ch = merscope_analysis_ch + .join(xenium_analysis_ch) + .map { + _branchKey, pairId, segmentation, merscopePath, merscopeTableKey, + merscopeShapeKey, settings, xeniumPath, xeniumTableKey, xeniumShapeKey -> + tuple( + pairId, + segmentation, + merscopePath, + xeniumPath, + merscopeTableKey, + merscopeShapeKey, + xeniumTableKey, + xeniumShapeKey, + settings, + ) + } + + align_inputs_ch = paired_zarrs_ch + .filter { _pairId, _merscopePath, _xeniumPath, settings -> settings.run_align } + .map { pairId, merscopePath, xeniumPath, _settings -> + tuple(pairId, merscopePath, xeniumPath) + } + + alignment_task_results_ch = ALIGN(align_inputs_ch) + + alignment_published_results_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_alignment_results && !settings.run_align)) { + [] + } else { + def merscopeLatest = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + "MERSCOPE", + "latest/latest_spatialdata.zarr", + ), + "QC/enriched MERSCOPE latest zarr for ${pairId}", + ) + def xeniumLatest = requireExistingPath( + publishedDatasetPath( + params.outdir, + pairId, + "XENIUM", + "latest/latest_spatialdata.zarr", + ), + "QC/enriched XENIUM latest zarr for ${pairId}", + ) + def transformJson = requireExistingPath( + publishedPairPath( + params.outdir, + pairId, + "alignment/align_out/alignment_transform.json", + ), + "ALIGN transform JSON for ${pairId}", + ) + def coordsDir = requireExistingPath( + publishedPairPath( + params.outdir, + pairId, + "alignment/align_out/alignment_coords", + ), + "ALIGN coordinate directory for ${pairId}", + ) + [ + tuple( + pairId, + merscopeLatest.toRealPath().toString(), + xeniumLatest.toRealPath().toString(), + transformJson, + coordsDir, + ) + ] + } + } + + alignment_results_ch = alignment_task_results_ch.mix(alignment_published_results_ch) + + align_qc_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!settings.run_align_qc) { + [] + } else { + [tuple(pairId, true)] + } + } + + align_qc_inputs_ch = alignment_results_ch + .join(align_qc_gate_ch) + .map { pairId, merscopeLatest, xeniumLatest, transformJson, coordsDir, _runFlag -> + tuple(pairId, merscopeLatest, xeniumLatest, transformJson, coordsDir) + } + + alignment_qc_results_ch = ALIGN_QC(align_qc_inputs_ch) + + alignment_done_no_qc_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_alignment_downstream && !settings.run_align_qc)) { + [] + } else { + [tuple(pairId, settings.analysis_segmentations)] + } + } + + alignment_done_no_qc_ch = alignment_results_ch + .join(alignment_done_no_qc_gate_ch) + .flatMap { + pairId, _merscopeLatest, _xeniumLatest, _transformJson, _coordsDir, analysisSegmentations -> + analysisSegmentations.collect { segmentation -> + tuple("${pairId}|${segmentation}", true) + } + } + + alignment_done_after_qc_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.need_alignment_downstream && settings.run_align_qc)) { + [] + } else { + [tuple(pairId, settings.analysis_segmentations)] + } + } + + alignment_done_after_qc_ch = alignment_qc_results_ch + .join(alignment_done_after_qc_gate_ch) + .flatMap { pairId, _alignQcOut, analysisSegmentations -> + analysisSegmentations.collect { segmentation -> + tuple("${pairId}|${segmentation}", true) + } + } + + alignment_done_branch_ch = alignment_done_no_qc_ch.mix(alignment_done_after_qc_ch) + + paired_downstream_no_align_ch = paired_analysis_zarrs_ch + .filter { + _pairId, _segmentation, _merscopePath, _xeniumPath, _merscopeTableKey, + _merscopeShapeKey, _xeniumTableKey, _xeniumShapeKey, settings -> + !settings.enable_alignment && + (settings.run_compare || + settings.run_visualize || + settings.run_spatial_gene_analysis || + settings.run_clustering_squidpy) + } + + paired_downstream_align_ch = paired_analysis_zarrs_ch + .filter { + _pairId, _segmentation, _merscopePath, _xeniumPath, _merscopeTableKey, + _merscopeShapeKey, _xeniumTableKey, _xeniumShapeKey, settings -> + settings.enable_alignment && + (settings.run_compare || + settings.run_visualize || + settings.run_spatial_gene_analysis || + settings.run_clustering_squidpy) + } + .map { + pairId, segmentation, merscopePath, xeniumPath, merscopeTableKey, + merscopeShapeKey, xeniumTableKey, xeniumShapeKey, settings -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + merscopePath, + xeniumPath, + merscopeTableKey, + merscopeShapeKey, + xeniumTableKey, + xeniumShapeKey, + settings, + ) + } + .join(alignment_done_branch_ch) + .map { + _branchKey, pairId, segmentation, merscopePath, xeniumPath, + merscopeTableKey, merscopeShapeKey, xeniumTableKey, + xeniumShapeKey, settings, _doneFlag -> + tuple( + pairId, + segmentation, + merscopePath, + xeniumPath, + merscopeTableKey, + merscopeShapeKey, + xeniumTableKey, + xeniumShapeKey, + settings, + ) + } + + paired_downstream_zarrs_ch = paired_downstream_no_align_ch.mix(paired_downstream_align_ch) + + compare_inputs_ch = paired_downstream_zarrs_ch + .filter { + _pairId, _segmentation, _merscopePath, _xeniumPath, _merscopeTableKey, + _merscopeShapeKey, _xeniumTableKey, _xeniumShapeKey, settings -> + settings.run_compare + } + .map { + pairId, segmentation, merscopePath, xeniumPath, merscopeTableKey, + _merscopeShapeKey, xeniumTableKey, _xeniumShapeKey, _settings -> + tuple( + pairId, + segmentation, + merscopePath, + xeniumPath, + merscopeTableKey, + xeniumTableKey, + ) + } + + compare_results_ch = COMPARE(compare_inputs_ch) + + compare_done_ch = compare_results_ch.map { pairId, segmentation, _compareOutDir -> + tuple("${pairId}|${segmentation}", true) + } + + analysis_samples_from_aligned_pairs_ch = paired_downstream_zarrs_ch + .filter { + _pairId, _segmentation, _merscopePath, _xeniumPath, _merscopeTableKey, + _merscopeShapeKey, _xeniumTableKey, _xeniumShapeKey, settings -> + settings.enable_alignment && + ( + settings.run_visualize || + settings.run_spatial_gene_analysis || + settings.run_clustering_squidpy + ) + } + .map { + pairId, segmentation, merscopePath, xeniumPath, _merscopeTableKey, + _merscopeShapeKey, _xeniumTableKey, _xeniumShapeKey, settings -> + tuple( + pairId, + segmentation, + samplesJsonForSegmentation( + pairId, + ["MERSCOPE", "XENIUM"], + ["MERSCOPE": merscopePath, "XENIUM": xeniumPath], + segmentation, + ), + settings, + ) + } + + analysis_samples_from_dataset_ch = analysis_dataset_zarrs_ch + .filter { + _pairId, _segmentation, _platform, _zarrPath, _tableKey, _shapeKey, settings -> + (!settings.paired_mode || !settings.enable_alignment) && + ( + settings.run_visualize || + settings.run_spatial_gene_analysis || + settings.run_clustering_squidpy + ) + } + .map { + pairId, segmentation, platform, zarrPath, _tableKey, _shapeKey, settings -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + settings, + platform, + zarrPath, + ) + } + .groupTuple() + .map { _branchKey, pairIds, segmentations, settingsList, platformNames, zarrPaths -> + def pairId = pairIds[0].toString() + def segmentation = segmentations[0].toString() + def settings = settingsList[0] + def pathsByPlatform = [:] + platformNames.eachWithIndex { platform, idx -> + pathsByPlatform[platform.toString()] = zarrPaths[idx].toString() + } + tuple( + pairId, + segmentation, + samplesJsonForSegmentation( + pairId, + settings.active_platforms, + pathsByPlatform, + segmentation, + ), + settings, + ) + } + + analysis_samples_ch = + analysis_samples_from_aligned_pairs_ch.mix(analysis_samples_from_dataset_ch) + + visualize_without_compare_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_visualize && !settings.run_compare + } + .map { pairId, segmentation, samplesJson, _settings -> + tuple(pairId, segmentation, samplesJson) + } + + visualize_after_compare_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_visualize && settings.run_compare + } + .map { pairId, segmentation, samplesJson, _settings -> + tuple("${pairId}|${segmentation}", pairId, segmentation, samplesJson) + } + .join(compare_done_ch) + .map { _branchKey, pairId, segmentation, samplesJson, _doneFlag -> + tuple(pairId, segmentation, samplesJson) + } + + visualize_inputs_ch = visualize_without_compare_ch.mix(visualize_after_compare_ch) + + visualize_results_ch = VISUALIZE(visualize_inputs_ch) + + visualize_done_ch = visualize_results_ch.map { pairId, segmentation, _visualizeOutDir -> + tuple("${pairId}|${segmentation}", true) + } + + spatial_gene_analysis_without_visualize_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_spatial_gene_analysis && !settings.run_visualize + } + .map { pairId, segmentation, samplesJson, settings -> + tuple( + pairId, + segmentation, + samplesJson, + settings.spatial_gene_analysis_transcript_analysis_enabled, + ) + } + + spatial_gene_analysis_after_visualize_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_spatial_gene_analysis && settings.run_visualize + } + .map { pairId, segmentation, samplesJson, settings -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + samplesJson, + settings.spatial_gene_analysis_transcript_analysis_enabled, + ) + } + .join(visualize_done_ch) + .map { + _branchKey, + pairId, + segmentation, + samplesJson, + transcriptAnalysisEnabled, + _doneFlag -> + tuple( + pairId, + segmentation, + samplesJson, + transcriptAnalysisEnabled, + ) + } + + spatial_gene_analysis_unannotated_inputs_ch = + spatial_gene_analysis_without_visualize_ch.mix( + spatial_gene_analysis_after_visualize_ch + ) + + spatial_gene_annotation_rows_ch = sample_rows_ch.flatMap { + pairId, row, settings -> + if (!settings.run_spatial_gene_analysis) { + [] + } else { + settings.analysis_segmentations.collect { segmentation -> + tuple("${pairId}|${segmentation}", row) + } + } + } + + spatial_gene_analysis_inputs_ch = spatial_gene_analysis_unannotated_inputs_ch + .map { pairId, segmentation, samplesJson, transcriptAnalysisEnabled -> + tuple( + "${pairId}|${segmentation}", + pairId, + segmentation, + samplesJson, + transcriptAnalysisEnabled, + ) + } + .join(spatial_gene_annotation_rows_ch) + .map { + _key, + pairId, + segmentation, + samplesJson, + transcriptAnalysisEnabled, + row -> + tuple( + pairId, + segmentation, + spatialGeneSamplesJson(samplesJson, row), + transcriptAnalysisEnabled, + ) + } + + spatial_gene_analysis_results_ch = SPATIAL_GENE_ANALYSIS( + spatial_gene_analysis_inputs_ch + ) + + spatial_gene_analysis_done_ch = spatial_gene_analysis_results_ch.map { + pairId, segmentation, _spatialGeneAnalysisOutDir -> + tuple("${pairId}|${segmentation}", true) + } + + clustering_without_dependencies_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_clustering_squidpy && + !settings.run_visualize && + !settings.run_spatial_gene_analysis + } + .map { pairId, segmentation, samplesJson, _settings -> + tuple(pairId, segmentation, samplesJson) + } + + clustering_after_visualize_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_clustering_squidpy && + settings.run_visualize && + !settings.run_spatial_gene_analysis + } + .map { pairId, segmentation, samplesJson, _settings -> + tuple("${pairId}|${segmentation}", pairId, segmentation, samplesJson) + } + .join(visualize_done_ch) + .map { _branchKey, pairId, segmentation, samplesJson, _doneFlag -> + tuple(pairId, segmentation, samplesJson) + } + + clustering_after_spatial_gene_analysis_ch = analysis_samples_ch + .filter { _pairId, _segmentation, _samplesJson, settings -> + settings.run_clustering_squidpy && settings.run_spatial_gene_analysis + } + .map { pairId, segmentation, samplesJson, _settings -> + tuple("${pairId}|${segmentation}", pairId, segmentation, samplesJson) + } + .join(spatial_gene_analysis_done_ch) + .map { _branchKey, pairId, segmentation, samplesJson, _doneFlag -> + tuple(pairId, segmentation, samplesJson) + } + + clustering_inputs_ch = + clustering_without_dependencies_ch + .mix(clustering_after_visualize_ch) + .mix(clustering_after_spatial_gene_analysis_ch) + + clustering_prepared_ch = CLUSTERING_SQUIDPY_PREPARE(clustering_inputs_ch) + clustering_computed_ch = CLUSTERING_SQUIDPY_COMPUTE(clustering_prepared_ch) + clustering_results_ch = CLUSTERING_SQUIDPY_FINALIZE(clustering_computed_ch) + + // Cortical depth runs after clustering so the per-cell broad_class and + // subcluster_label annotations exist and its depth violin plots can be + // produced. The distance-from-object branch consumes its tissue-region + // annotation when both opt-in stages run in the same invocation. + clustering_done_per_pair_ch = clustering_results_ch + .map { pairId, segmentation, _samplesJson, _clusteringOutDir -> + tuple(pairId, segmentation) + } + .groupTuple() + .map { pairId, _segmentations -> tuple(pairId, true) } + + compute_cortical_depth_after_clustering_gate_ch = + sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_compute_cortical_depth && settings.run_clustering_squidpy)) { + [] + } else { + def corticalSegmentations = settings.analysis_segmentations + if (settings.run_distance_from_object) { + corticalSegmentations = ( + corticalSegmentations + settings.distance_from_object_segmentations + ).unique() + } + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", row, corticalSegmentations) + } + } + } + + compute_cortical_depth_after_clustering_ch = downstream_zarrs_ch + .join(compute_cortical_depth_after_clustering_gate_ch) + .map { key, pairId, platform, latestZarr, row, analysisSegmentations -> + tuple(pairId, key, platform, latestZarr, row, analysisSegmentations) + } + .combine(clustering_done_per_pair_ch, by: 0) + .map { pairId, key, platform, latestZarr, row, analysisSegmentations, _done -> + tuple(key, pairId, platform, latestZarr, row, analysisSegmentations) + } + + // Fallback for invocations where cortical depth is enabled but clustering is + // not run in the same pass (e.g. --only_stage compute_cortical_depth after a + // prior full run): read the enriched/quantified zarr directly. Violins are + // produced when a clustering table already exists in that zarr, otherwise + // they are skipped with a logged note. + compute_cortical_depth_without_clustering_gate_ch = + sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_compute_cortical_depth && !settings.run_clustering_squidpy)) { + [] + } else { + def corticalSegmentations = settings.analysis_segmentations + if (settings.run_distance_from_object) { + corticalSegmentations = ( + corticalSegmentations + settings.distance_from_object_segmentations + ).unique() + } + settings.active_platforms.collect { platform -> + tuple("${pairId}|${platform}", row, corticalSegmentations) + } + } + } + + compute_cortical_depth_without_clustering_ch = downstream_zarrs_ch + .join(compute_cortical_depth_without_clustering_gate_ch) + .map { key, pairId, platform, latestZarr, row, analysisSegmentations -> + tuple(key, pairId, platform, latestZarr, row, analysisSegmentations) + } + + compute_cortical_depth_inputs_ch = compute_cortical_depth_after_clustering_ch + .mix(compute_cortical_depth_without_clustering_ch) + .map { key, pairId, platform, latestZarr, row, analysisSegmentations -> + def depthConfig = corticalDepthConfigForPlatform( + row, + pairId, + platform, + analysisSegmentations, + params, + ) + tuple( + key, + pairId, + platform, + groovy.json.JsonOutput.prettyPrint(groovy.json.JsonOutput.toJson(depthConfig)), + latestZarr, + ) + } + + compute_cortical_depth_results_ch = COMPUTE_CORTICAL_DEPTH( + compute_cortical_depth_inputs_ch + ) + + distance_from_object_after_cortical_gate_ch = + sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_distance_from_object && + settings.run_compute_cortical_depth)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple( + "${pairId}|${platform}", + row, + settings.distance_from_object_segmentations, + ) + } + } + } + + distance_from_object_after_cortical_ch = compute_cortical_depth_results_ch + .join(distance_from_object_after_cortical_gate_ch) + .map { + key, + pairId, + platform, + latestZarr, + _corticalDepthOut, + row, + distanceSegmentations -> + tuple( + key, + pairId, + platform, + latestZarr, + row, + distanceSegmentations, + ) + } + + distance_from_object_after_clustering_gate_ch = + sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_distance_from_object && + !settings.run_compute_cortical_depth && + settings.run_clustering_squidpy)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple( + "${pairId}|${platform}", + row, + settings.distance_from_object_segmentations, + ) + } + } + } + + distance_from_object_after_clustering_ch = downstream_zarrs_ch + .join(distance_from_object_after_clustering_gate_ch) + .map { + key, pairId, platform, latestZarr, row, distanceSegmentations -> + tuple( + pairId, + key, + platform, + latestZarr, + row, + distanceSegmentations, + ) + } + .combine(clustering_done_per_pair_ch, by: 0) + .map { + pairId, + key, + platform, + latestZarr, + row, + distanceSegmentations, + _done -> + tuple( + key, + pairId, + platform, + latestZarr, + row, + distanceSegmentations, + ) + } + + distance_from_object_without_terminal_dependency_gate_ch = + sample_rows_ch.flatMap { pairId, row, settings -> + if (!(settings.run_distance_from_object && + !settings.run_compute_cortical_depth && + !settings.run_clustering_squidpy)) { + [] + } else { + settings.active_platforms.collect { platform -> + tuple( + "${pairId}|${platform}", + row, + settings.distance_from_object_segmentations, + ) + } + } + } + + distance_from_object_without_terminal_dependency_ch = downstream_zarrs_ch + .join(distance_from_object_without_terminal_dependency_gate_ch) + .map { + key, pairId, platform, latestZarr, row, distanceSegmentations -> + tuple( + key, + pairId, + platform, + latestZarr, + row, + distanceSegmentations, + ) + } + + distance_from_object_inputs_ch = distance_from_object_after_cortical_ch + .mix(distance_from_object_after_clustering_ch) + .mix(distance_from_object_without_terminal_dependency_ch) + .map { + key, pairId, platform, latestZarr, row, distanceSegmentations -> + def distanceConfig = distanceFromObjectConfigForPlatform( + pairId, + platform, + distanceSegmentations, + params, + ) + def annotationPath = file( + normalizedPath(distanceFromObjectAnnotationPath(row, platform)) + ) + tuple( + key, + pairId, + platform, + distanceSegmentations, + groovy.json.JsonOutput.prettyPrint( + groovy.json.JsonOutput.toJson(distanceConfig) + ), + latestZarr, + annotationPath, + ) + } + + distance_from_object_annotation_results_ch = DISTANCE_FROM_OBJECT_ANNOTATE( + distance_from_object_inputs_ch + ) + + distance_from_object_cohort_inputs_ch = distance_from_object_annotation_results_ch + .map { + _key, _pairId, platform, distanceSegmentations, + _latestZarr, annotationOutputDir -> + tuple(platform, annotationOutputDir, distanceSegmentations) + } + .groupTuple() + .map { platform, annotationOutputDirs, segmentationLists -> + def annotatedSegmentations = segmentationLists + .flatten() + .collect { value -> value.toString() } + .unique() + def cohortConfig = [ + platform: platform, + annotation_output_dirs: ["pair_outputs"], + output_dir: "distance_from_object_cohort_out", + segmentations: annotatedSegmentations, + min_pairs: params.distance_from_object_min_pairs, + n_cpus: params.distance_from_object_n_cpus, + ] + tuple( + platform, + groovy.json.JsonOutput.prettyPrint( + groovy.json.JsonOutput.toJson(cohortConfig) + ), + annotationOutputDirs, + ) + } + + DISTANCE_FROM_OBJECT_COHORT(distance_from_object_cohort_inputs_ch) + + mapmycells_after_clustering_gate_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.run_mapmycells && settings.run_clustering_squidpy)) { + [] + } else { + validateMapMyCellsParams(params) + settings.analysis_segmentations.collect { segmentation -> + tuple("${pairId}|${segmentation}", true) + } + } + } + + mapmycells_from_clustering_ch = clustering_results_ch + .map { pairId, segmentation, samplesJson, clusteringOutDir -> + tuple("${pairId}|${segmentation}", pairId, segmentation, samplesJson, clusteringOutDir) + } + .join(mapmycells_after_clustering_gate_ch) + .map { _branchKey, pairId, segmentation, samplesJson, clusteringOutDir, _runFlag -> + tuple(pairId, segmentation, samplesJson, clusteringOutDir) + } + + mapmycells_published_ch = sample_rows_ch.flatMap { pairId, _row, settings -> + if (!(settings.run_mapmycells && !settings.run_clustering_squidpy)) { + [] + } else { + validateMapMyCellsParams(params) + settings.analysis_segmentations.collect { segmentation -> + def clusteringOut = requireExistingPath( + publishedPairPath( + params.outdir, + pairId, + "${segmentation}/clustering_squidpy/clustering_squidpy_out", + ), + "CLUSTERING_SQUIDPY ${segmentation} output directory for ${pairId}", + ) + tuple( + pairId, + segmentation, + samplesJsonForSegmentation( + pairId, + settings.active_platforms, + [:], + segmentation, + ), + clusteringOut, + ) + } + } + } + + mapmycells_inputs_ch = mapmycells_from_clustering_ch.mix(mapmycells_published_ch) + + MAPMYCELLS(mapmycells_inputs_ch) +} diff --git a/merxen_stub/workflows/modules/alignment.nf b/merxen_stub/workflows/modules/alignment.nf new file mode 100644 index 0000000..c6c7dcd --- /dev/null +++ b/merxen_stub/workflows/modules/alignment.nf @@ -0,0 +1,48 @@ +process ALIGN { + tag "${pair_id}" + + publishDir { "${params.outdir}/${pair_id}/alignment" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(merscope_zarr_path), + val(xenium_zarr_path) + + output: + tuple val(pair_id), + val(merscope_zarr_path), + val(xenium_zarr_path), + path("align_out/alignment_transform.json"), + path("align_out/alignment_coords") + + script: + """ + set -euo pipefail + mkdir -p align_out/alignment_coords + echo '{"stub": "${pair_id} alignment transform"}' > align_out/alignment_transform.json + echo "x,y" > align_out/alignment_coords/${pair_id}_raw_centroids.csv + """ +} + +process ALIGN_QC { + tag "${pair_id}" + + publishDir { "${params.outdir}/${pair_id}/alignment_qc" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(merscope_zarr), + val(xenium_zarr), + path(transform_json), + path(coords_dir) + + output: + tuple val(pair_id), path("alignment_qc_out") + + script: + """ + set -euo pipefail + mkdir -p alignment_qc_out + echo '{"stub": "${pair_id} alignment qc"}' > alignment_qc_out/${pair_id}_alignment_qc.json + """ +} diff --git a/merxen_stub/workflows/modules/clustering_squidpy.nf b/merxen_stub/workflows/modules/clustering_squidpy.nf new file mode 100644 index 0000000..cce3c56 --- /dev/null +++ b/merxen_stub/workflows/modules/clustering_squidpy.nf @@ -0,0 +1,72 @@ +process CLUSTERING_SQUIDPY_PREPARE { + tag "${pair_id}:${segmentation}" + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json) + + output: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path("clustering_squidpy_config.json"), + path("clustering_prepare_out") + + script: + """ + set -euo pipefail + echo '{"stub": "${pair_id} ${segmentation} clustering config"}' > clustering_squidpy_config.json + mkdir -p clustering_prepare_out + echo "stub prepared clustering input" > clustering_prepare_out/prepared.txt + """ +} + +process CLUSTERING_SQUIDPY_COMPUTE { + tag "${pair_id}:${segmentation}" + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path(clustering_config), + path(prepared_dir) + + output: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path("clustering_compute_out") + + script: + """ + set -euo pipefail + mkdir -p clustering_compute_out + echo "stub computed clustering" > clustering_compute_out/computed.txt + """ +} + +process CLUSTERING_SQUIDPY_FINALIZE { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/clustering_squidpy" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path(computed_dir) + + output: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path("clustering_squidpy_out") + + script: + """ + set -euo pipefail + mkdir -p clustering_squidpy_out + echo "stub clustered outputs for ${pair_id} ${segmentation}" > clustering_squidpy_out/clustering_manifest.txt + """ +} diff --git a/merxen_stub/workflows/modules/comparison.nf b/merxen_stub/workflows/modules/comparison.nf new file mode 100644 index 0000000..c14ed9d --- /dev/null +++ b/merxen_stub/workflows/modules/comparison.nf @@ -0,0 +1,23 @@ +process COMPARE { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/comparison" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(merscope_zarr), + val(xenium_zarr), + val(merscope_table_key), + val(xenium_table_key) + + output: + tuple val(pair_id), val(segmentation), path("compare_out") + + script: + """ + set -euo pipefail + mkdir -p compare_out + echo '{"stub": "${pair_id} ${segmentation} comparison metrics"}' > compare_out/${pair_id}_comparison_metrics.json + """ +} diff --git a/merxen_stub/workflows/modules/compute_cortical_depth.nf b/merxen_stub/workflows/modules/compute_cortical_depth.nf new file mode 100644 index 0000000..07d871d --- /dev/null +++ b/merxen_stub/workflows/modules/compute_cortical_depth.nf @@ -0,0 +1,29 @@ +process COMPUTE_CORTICAL_DEPTH { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/compute_cortical_depth" }, mode: "symlink", overwrite: true + + input: + tuple val(key), + val(pair_id), + val(platform), + val(cortical_depth_config_json), + path(latest_zarr) + + output: + tuple val(key), + val(pair_id), + val(platform), + path("latest_input.zarr"), + path("compute_cortical_depth_out") + + script: + """ + set -euo pipefail + if [[ ! -e latest_input.zarr ]]; then + ln -s ${latest_zarr} latest_input.zarr + fi + mkdir -p compute_cortical_depth_out + echo '{"stub": "${pair_id} ${platform} cortical depth qc summary"}' > compute_cortical_depth_out/cortical_depth_qc_summary.json + """ +} diff --git a/merxen_stub/workflows/modules/distance_from_object.nf b/merxen_stub/workflows/modules/distance_from_object.nf new file mode 100644 index 0000000..e8d1094 --- /dev/null +++ b/merxen_stub/workflows/modules/distance_from_object.nf @@ -0,0 +1,55 @@ +process DISTANCE_FROM_OBJECT_ANNOTATE { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/distance_from_object" }, mode: "symlink", overwrite: true + + input: + tuple val(key), + val(pair_id), + val(platform), + val(distance_segmentations), + val(distance_config_json), + path(latest_zarr), + path(object_annotations) + + output: + tuple val(key), + val(pair_id), + val(platform), + val(distance_segmentations), + path("latest_input.zarr"), + path("distance_from_object_out") + + script: + """ + set -euo pipefail + if [[ ! -e latest_input.zarr ]]; then + ln -s ${latest_zarr} latest_input.zarr + fi + mkdir -p distance_from_object_out + echo '{"stub": "${pair_id} ${platform} distance from object summary"}' > distance_from_object_out/distance_from_object_summary.json + echo '{"stub": "registered object annotations"}' > distance_from_object_out/registered_object_annotations.geojson + """ +} + + +process DISTANCE_FROM_OBJECT_COHORT { + tag "${platform}" + + publishDir { "${params.outdir}/distance_from_object/cohort/${platform.toLowerCase()}" }, mode: "symlink", overwrite: true + + input: + tuple val(platform), + val(cohort_config_json), + path(annotation_output_dirs, stageAs: "pair_outputs/dir??/*") + + output: + tuple val(platform), path("distance_from_object_cohort_out") + + script: + """ + set -euo pipefail + mkdir -p distance_from_object_cohort_out + echo '{"stub": "${platform} distance from object cohort summary"}' > distance_from_object_cohort_out/distance_from_object_cohort_summary.json + """ +} diff --git a/merxen_stub/workflows/modules/enrichment.nf b/merxen_stub/workflows/modules/enrichment.nf new file mode 100644 index 0000000..d5f8725 --- /dev/null +++ b/merxen_stub/workflows/modules/enrichment.nf @@ -0,0 +1,19 @@ +process ENRICH { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/enrichment" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(enrich_config_json), path(latest_zarr), path(mask_path), path(nuclei_mask_path) + + output: + tuple val(key), val(pair_id), val(platform), path("latest_input.zarr"), path("enrich_out") + + script: + """ + set -euo pipefail + ln -s ${latest_zarr} latest_input.zarr + mkdir -p enrich_out + echo "shape,transcripts_assigned" > enrich_out/${pair_id}_${platform}_assignment_summary.csv + """ +} diff --git a/merxen_stub/workflows/modules/mapmycells.nf b/merxen_stub/workflows/modules/mapmycells.nf new file mode 100644 index 0000000..93d38ea --- /dev/null +++ b/merxen_stub/workflows/modules/mapmycells.nf @@ -0,0 +1,21 @@ +process MAPMYCELLS { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/mapmycells" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path(clustering_out_dir, stageAs: "clustering_squidpy_input") + + output: + tuple val(pair_id), val(segmentation), path("mapmycells_out") + + script: + """ + set -euo pipefail + mkdir -p mapmycells_out + echo '{"stub": "${pair_id} ${segmentation} mapmycells manifest"}' > mapmycells_out/${pair_id}_mapmycells_manifest.json + """ +} diff --git a/merxen_stub/workflows/modules/mask_image_quantification.nf b/merxen_stub/workflows/modules/mask_image_quantification.nf new file mode 100644 index 0000000..9e12352 --- /dev/null +++ b/merxen_stub/workflows/modules/mask_image_quantification.nf @@ -0,0 +1,30 @@ +process MASK_IMAGE_QUANTIFICATION { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/mask_image_quantification" }, mode: "symlink", overwrite: true + + input: + tuple val(key), + val(pair_id), + val(platform), + val(mask_image_quantification_config_json), + path(latest_zarr), + path(mask_path) + + output: + tuple val(key), + val(pair_id), + val(platform), + path("latest_input.zarr"), + path("mask_image_quantification_out") + + script: + """ + set -euo pipefail + if [[ ! -e latest_input.zarr ]]; then + ln -s ${latest_zarr} latest_input.zarr + fi + mkdir -p mask_image_quantification_out + echo "stub mask image quantification for ${pair_id}_${platform}" > mask_image_quantification_out/${pair_id}_${platform}_mask_image_quantification_summary.json + """ +} diff --git a/merxen_stub/workflows/modules/mecr.nf b/merxen_stub/workflows/modules/mecr.nf new file mode 100644 index 0000000..31828a9 --- /dev/null +++ b/merxen_stub/workflows/modules/mecr.nf @@ -0,0 +1,44 @@ +process MECR_REFERENCE { + tag "WHB-10Xv3" + + publishDir { "${params.outdir}/mecr_reference" }, mode: "copy", overwrite: true + + input: + val(samples_json) + + output: + path("mecr_reference_out") + + script: + """ + set -euo pipefail + mkdir -p mecr_reference_out + echo "gene,marker" > mecr_reference_out/mecr_reference_markers.csv + echo '{"stub": "mecr reference manifest"}' > mecr_reference_out/mecr_reference_manifest.json + """ +} + + +process MECR { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/mecr" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json), + path(reference_out) + + output: + tuple val(pair_id), + val(segmentation), + path("mecr_out") + + script: + """ + set -euo pipefail + mkdir -p mecr_out + echo '{"stub": "${pair_id} ${segmentation} mecr summary"}' > mecr_out/${pair_id}_mecr_summary.csv + """ +} diff --git a/merxen_stub/workflows/modules/proseg_bootstrap.nf b/merxen_stub/workflows/modules/proseg_bootstrap.nf new file mode 100644 index 0000000..681973a --- /dev/null +++ b/merxen_stub/workflows/modules/proseg_bootstrap.nf @@ -0,0 +1,15 @@ +process ENSURE_PROSEG { + tag "proseg" + + input: + val trigger + + output: + path("proseg_path.txt") + + script: + """ + set -euo pipefail + printf '%s\\n' /usr/bin/true > proseg_path.txt + """ +} diff --git a/merxen_stub/workflows/modules/qc.nf b/merxen_stub/workflows/modules/qc.nf new file mode 100644 index 0000000..7bc5581 --- /dev/null +++ b/merxen_stub/workflows/modules/qc.nf @@ -0,0 +1,64 @@ +process VALIDATE_ANALYSIS_LAYER { + tag "${pair_id}:${platform}:${segmentation}" + + input: + tuple val(key), + val(pair_id), + val(platform), + val(segmentation), + path(latest_zarr), + val(table_key), + val(shape_key), + val(settings) + + output: + tuple val(key), + val(pair_id), + val(platform), + val(segmentation), + path(latest_zarr), + val(table_key), + val(shape_key), + val(settings), + path("analysis_layer_validation.json") + + script: + """ + set -euo pipefail + echo '{"stub": "${pair_id}_${platform} ${segmentation} analysis layer valid"}' > analysis_layer_validation.json + """ +} + + +process QC { + tag "${pair_id}:${platform}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/${segmentation}/qc" }, mode: "symlink", overwrite: true + + input: + tuple val(key), + val(pair_id), + val(platform), + val(segmentation), + path(latest_zarr), + val(table_key), + val(shape_key) + + output: + tuple val(key), + val(pair_id), + val(platform), + val(segmentation), + path(latest_zarr), + path("qc_out"), + val(table_key), + val(shape_key) + + script: + """ + set -euo pipefail + mkdir -p qc_out + dataset=\$(echo "${pair_id}_${platform}" | tr '[:upper:]' '[:lower:]') + echo "metric,value" > qc_out/\${dataset}_qc_summary.csv + """ +} diff --git a/merxen_stub/workflows/modules/segmentation.nf b/merxen_stub/workflows/modules/segmentation.nf new file mode 100644 index 0000000..f8e8d25 --- /dev/null +++ b/merxen_stub/workflows/modules/segmentation.nf @@ -0,0 +1,124 @@ +process CELLPOSE_SEGMENT { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/segmentation" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(seg_config_json) + + output: + tuple val(key), val(pair_id), val(platform), val(seg_config_json), path("segment_out/cellpose_masks_tiled.npy"), path("segment_out/cellpose_cellprobs_tiled.npy"), path("segment_out/transcripts_for_proseg.csv"), path("segment_out/cellpose_transforms.json"), path("segment_out/cellpose_stitching_stats.json") + + script: + """ + set -euo pipefail + mkdir -p segment_out + echo "stub cellpose mask for ${pair_id}_${platform}" > segment_out/cellpose_masks_tiled.npy + echo "stub cellpose cellprobs for ${pair_id}_${platform}" > segment_out/cellpose_cellprobs_tiled.npy + echo "cell_id,x,y" > segment_out/transcripts_for_proseg.csv + echo '{"stub": "cellpose transforms"}' > segment_out/cellpose_transforms.json + echo '{"stub": "cellpose stitching stats"}' > segment_out/cellpose_stitching_stats.json + """ +} + + +process CELLPOSE_NUCLEI_SEGMENT { + tag "${pair_id}:${platform}:nuclei" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/segmentation" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(seg_config_json) + + output: + tuple val(key), val(pair_id), val(platform), val(seg_config_json), path("segment_out/cellpose_nuclei_masks_tiled.npy"), path("segment_out/cellpose_nuclei_stitching_stats.json") + + script: + """ + set -euo pipefail + mkdir -p segment_out + echo "stub cellpose nuclei mask for ${pair_id}_${platform}" > segment_out/cellpose_nuclei_masks_tiled.npy + echo '{"stub": "cellpose nuclei stitching stats"}' > segment_out/cellpose_nuclei_stitching_stats.json + """ +} + + +process PROSEG_SEGMENT { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/segmentation" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(seg_config_json), path(cellpose_mask), path(cellpose_cellprob), path(transcripts_csv), path(cellpose_transforms), path(_stitching_stats), path(nuclei_mask), path(_nuclei_stitching_stats), path(proseg_path_file) + + output: + tuple val(key), val(pair_id), val(platform), path("segment_out/proseg_base_latest.zarr"), path(cellpose_mask), path(transcripts_csv), path(nuclei_mask) + + script: + """ + set -euo pipefail + mkdir -p segment_out/proseg_base_latest.zarr + echo '{"stub": "${pair_id}_${platform} proseg base latest spatialdata"}' > segment_out/proseg_base_latest.zarr/.zattrs + """ +} + + +workflow SEGMENT { + take: + segment_inputs + nuclei_results + proseg_path + + main: + cellpose_results = CELLPOSE_SEGMENT(segment_inputs) + combined_cellpose_results = cellpose_results + .join(nuclei_results) + .map { + key, pair_id, platform, seg_config_json, cellpose_mask, + cellpose_cellprob, transcripts_csv, cellpose_transforms, stitching_stats, + nuclei_pair_id, nuclei_platform, _nuclei_seg_config_json, + nuclei_mask, nuclei_stitching_stats -> + if (pair_id != nuclei_pair_id || platform != nuclei_platform) { + error("Cell/nuclei Cellpose channel mismatch for ${key}") + } + tuple( + key, + pair_id, + platform, + seg_config_json, + cellpose_mask, + cellpose_cellprob, + transcripts_csv, + cellpose_transforms, + stitching_stats, + nuclei_mask, + nuclei_stitching_stats, + ) + } + proseg_inputs = combined_cellpose_results + .combine(proseg_path) + .map { + key, pair_id, platform, seg_config_json, cellpose_mask, + cellpose_cellprob, transcripts_csv, cellpose_transforms, stitching_stats, nuclei_mask, + nuclei_stitching_stats, + proseg_path_file -> + tuple( + key, + pair_id, + platform, + seg_config_json, + cellpose_mask, + cellpose_cellprob, + transcripts_csv, + cellpose_transforms, + stitching_stats, + nuclei_mask, + nuclei_stitching_stats, + proseg_path_file, + ) + } + segment_results = PROSEG_SEGMENT(proseg_inputs) + + emit: + segment_results +} diff --git a/merxen_stub/workflows/modules/spatial_gene_analysis.nf b/merxen_stub/workflows/modules/spatial_gene_analysis.nf new file mode 100644 index 0000000..38c0df5 --- /dev/null +++ b/merxen_stub/workflows/modules/spatial_gene_analysis.nf @@ -0,0 +1,23 @@ +process SPATIAL_GENE_ANALYSIS { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/spatial_gene_analysis" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json), + val(transcript_analysis_enabled) + + output: + tuple val(pair_id), + val(segmentation), + path("spatial_gene_analysis_out") + + script: + """ + set -euo pipefail + mkdir -p spatial_gene_analysis_out + echo '{"stub": "${pair_id} ${segmentation} spatial gene analysis", "transcript_analysis_enabled": ${transcript_analysis_enabled}}' > spatial_gene_analysis_out/${pair_id}_spatial_gene_analysis_manifest.json + """ +} diff --git a/merxen_stub/workflows/modules/spatialdata_build.nf b/merxen_stub/workflows/modules/spatialdata_build.nf new file mode 100644 index 0000000..d154098 --- /dev/null +++ b/merxen_stub/workflows/modules/spatialdata_build.nf @@ -0,0 +1,74 @@ +process BUILD_SPATIALDATA { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/spatialdata" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(build_config_json) + + output: + tuple val(key), val(pair_id), val(platform), path("spatialdata_out/source_spatialdata.zarr") + + script: + """ + set -euo pipefail + + cat > build_config.json <<'JSON' +${build_config_json} +JSON + + # Mimic merxen build-spatialdata input resolution (io/builders/pipeline.py): + # input_path must exist and be an accessible directory — either a raw export + # folder or a reusable .zarr cache — otherwise the real stage raises + # FileNotFoundError. On Flow this is what fails if the selected data path is + # not mounted into the process container, so the stub checks it too. + input_path=\$(grep '"input_path"' build_config.json | head -1 | sed -E 's/.*"input_path": *"(.*)".*/\\1/') + + if [ -z "\${input_path}" ]; then + echo "ERROR: no input_path in build config for ${pair_id}:${platform}" >&2 + exit 1 + fi + if [ ! -e "\${input_path}" ]; then + echo "ERROR: ${platform} input_path is not accessible from this process: \${input_path}" >&2 + echo " (on Flow this usually means the data path is not mounted into the container)" >&2 + exit 1 + fi + if [ ! -d "\${input_path}" ]; then + echo "ERROR: ${platform} input_path exists but is not a directory: \${input_path}" >&2 + exit 1 + fi + + echo "[stub] ${platform} input_path=\${input_path}; contents:" >&2 + ls -A "\${input_path}" >&2 + + # A reusable SpatialData cache is just a .zarr directory; a raw export folder + # must contain the platform-specific files the builder reads. + case "\${input_path}" in + *.zarr) + : ;; + *) + missing="" + if [ "${platform}" = "MERSCOPE" ]; then + for req in images micron_to_mosaic_pixel_transform.csv cell_boundaries.parquet; do + [ -e "\${input_path}/\${req}" ] || missing="\${missing} \${req}" + done + if [ ! -e "\${input_path}/detected_transcripts.parquet" ] && [ ! -e "\${input_path}/detected_transcripts.csv" ]; then + missing="\${missing} detected_transcripts.parquet|csv" + fi + elif [ "${platform}" = "XENIUM" ]; then + if [ ! -e "\${input_path}/experiment.xenium" ] && [ ! -e "\${input_path}/specs.json" ] && [ ! -e "\${input_path}/specs/specs.json" ]; then + missing="\${missing} experiment.xenium|specs.json" + fi + [ -e "\${input_path}/transcripts.parquet" ] || missing="\${missing} transcripts.parquet" + fi + if [ -n "\${missing}" ]; then + echo "ERROR: ${platform} raw input \${input_path} is missing expected entries:\${missing}" >&2 + exit 1 + fi + ;; + esac + + mkdir -p spatialdata_out/source_spatialdata.zarr + echo '{"stub": "${pair_id}_${platform} source spatialdata"}' > spatialdata_out/source_spatialdata.zarr/.zattrs + """ +} diff --git a/merxen_stub/workflows/modules/viewer_cache.nf b/merxen_stub/workflows/modules/viewer_cache.nf new file mode 100644 index 0000000..1b72983 --- /dev/null +++ b/merxen_stub/workflows/modules/viewer_cache.nf @@ -0,0 +1,21 @@ +process VIEWER_CACHE { + tag "${pair_id}:${platform}" + + publishDir { "${params.outdir}/${pair_id}/${platform.toLowerCase()}/viewer_cache" }, mode: "symlink", overwrite: true + + input: + tuple val(key), val(pair_id), val(platform), val(viewer_cache_config_json), path(latest_zarr) + + output: + tuple val(key), val(pair_id), val(platform), path("latest_input.zarr"), path("viewer_cache_out") + + script: + """ + set -euo pipefail + if [[ ! -e latest_input.zarr ]]; then + ln -s ${latest_zarr} latest_input.zarr + fi + mkdir -p viewer_cache_out + echo "stub viewer cache for ${pair_id}_${platform}" > viewer_cache_out/viewer_cache_manifest.json + """ +} diff --git a/merxen_stub/workflows/modules/visualization.nf b/merxen_stub/workflows/modules/visualization.nf new file mode 100644 index 0000000..f708af4 --- /dev/null +++ b/merxen_stub/workflows/modules/visualization.nf @@ -0,0 +1,20 @@ +process VISUALIZE { + tag "${pair_id}:${segmentation}" + + publishDir { "${params.outdir}/${pair_id}/${segmentation}/visualization" }, mode: "copy", overwrite: true + + input: + tuple val(pair_id), + val(segmentation), + val(samples_json) + + output: + tuple val(pair_id), val(segmentation), path("visualize_out") + + script: + """ + set -euo pipefail + mkdir -p visualize_out + echo "stub visualization for ${pair_id} ${segmentation}" > visualize_out/${pair_id}_visualization.txt + """ +} diff --git a/merxen_stub/workflows/nextflow.config b/merxen_stub/workflows/nextflow.config new file mode 100644 index 0000000..0331a59 --- /dev/null +++ b/merxen_stub/workflows/nextflow.config @@ -0,0 +1,481 @@ +params { + samplesheet = null + outdir = "./results" + force_spatialdata_build = false + force_proseg_rerun = false + analysis_mode = "paired" + analysis_segmentation = "both" + start_stage = "build_spatialdata" + stop_stage = "clustering_squidpy" + only_stage = null + + // Cellpose defaults + cellpose_model_type = "cyto3" + cellpose_flow_threshold = 0.7 + cellpose_cellprob = -5.0 + // Execution profiles choose whether Cellpose uses a GPU. + cellpose_gpu = false + cellpose_diameter = null + cellpose_tile_overlap = 0.15 + cellpose_bsize = 256 + cellpose_tile_size_candidates = [6144, 4096, 3072, 2048] + cellpose_min_tile_size = 1024 + cellpose_stitch_overlap_px = 256 + cellpose_stitch_status_every_tiles = 10 + cellpose_filter_per_tile = true + cellpose_duplicate_iou_threshold = 0.25 + cellpose_duplicate_overlap_fraction = 0.5 + cellpose_min_remaining_fraction = 0.05 + cellpose_edge_touch_policy = "keep" + cellpose_write_stitching_stats = true + cellpose_final_min_area_um2 = 5.0 + cellpose_final_max_area_um2 = 400.0 + cellpose_final_filter_chunk_mb = 256 + cellpose_use_bfloat16 = false + + // ProSeg defaults + proseg_binary = null + proseg_search_paths = [] + proseg_install_path = null + proseg_auto_install = true + proseg_cargo_package = "proseg" + proseg_version = "3.2.0" + proseg_git_url = "https://github.com/dcjones/proseg.git" + proseg_git_rev = "e7df1eace923ce4c6ec70b2c597c5d126aa3db88" + proseg_samples = 1200 + proseg_voxel_size = 0.5 + proseg_burnin_voxel_size = 1.0 + proseg_nuclear_reassignment_prob = 0.20 + proseg_diffusion_probability = 0.20 + proseg_cell_compactness = 0.04 + proseg_num_threads = null + default_merscope_voxel_layers = 7 + default_xenium_voxel_layers = 2 + + // Transcript-supported Cellpose/ProSeg hybrid branch (enabled by default) + proseg_hybrid_enabled = true + proseg_hybrid_min_transcripts = 10 + proseg_hybrid_outlier_neighbors = 2 + proseg_hybrid_outlier_mad_multiplier = 2.0 + proseg_hybrid_minimum_external_group = 3 + proseg_hybrid_chain_radius_scale = 2.0 + proseg_hybrid_near_surface_radius_fraction = 0.25 + proseg_hybrid_maximum_expansion_radius_fraction = 1.0 + proseg_hybrid_attachment_arc_width_scale = 0.5 + proseg_hybrid_rounding_radius_fraction = 0.15 + proseg_hybrid_smoothing_radius_um = 10.0 + proseg_hybrid_outward_rounding_um = 0.2 + proseg_hybrid_smoothing_quad_segs = 32 + proseg_hybrid_containment_tolerance_um = 1.0e-5 + + // Cellpose mask image quantification defaults + mask_image_quantification_enabled = true + + // Viewer-cache defaults (pre-build the napari viewer's on-the-fly caches). + // These mirror the viewer's own defaults; keep them in lockstep with the + // viewer's runtime settings so a pre-built cache is reused rather than rebuilt. + viewer_cache_enabled = true + viewer_cache_downsample = 4 + viewer_cache_label_chunk_size = 2048 + viewer_cache_contour_width = 1 + viewer_cache_min_size = 4096 + viewer_cache_build_image_pyramid = true + + // Cortical-depth defaults (opt in with --cortical_depth_enabled true) + cortical_depth_enabled = false + cortical_depth_coordinate_unit_um = 1.0 + cortical_depth_raster_resolution_um = 5.0 + cortical_depth_raster_padding_um = null + cortical_depth_boundary_band_um = null + cortical_depth_boundary_smoothing_window = 0 + cortical_depth_streamline_spacing_um = 50.0 + cortical_depth_streamline_step_um = null + cortical_depth_streamline_max_steps = 4000 + cortical_depth_streamline_resample_points = 101 + cortical_depth_side_boundary_distance_um = 25.0 + cortical_depth_contour_levels = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] + cortical_depth_write_spatialdata_table = true + + // Polygon-object distance analysis (opt in after cortical-depth annotation) + distance_from_object_enabled = false + distance_from_object_segmentations = ["proseg", "original", "cellpose"] + distance_from_object_object_types = null + distance_from_object_coordinate_unit_um = 1.0 + distance_from_object_near_distance_um = 50.0 + distance_from_object_far_distance_um = 100.0 + distance_from_object_max_distance_um = 200.0 + distance_from_object_min_cells_per_pseudobulk = 10 + distance_from_object_min_pairs = 2 + distance_from_object_n_cpus = null + distance_from_object_write_spatialdata_table = true + + // Platform defaults + xenium_min_qv = 20.0 + + // Optional cross-section alignment defaults + enable_alignment = false + alignment_spateo_mode = "SN-S" + alignment_device = "cpu" + alignment_dtype = "float32" + alignment_selected_mode = "nonrigid" + alignment_max_iter = 360 + alignment_nonrigid_start_iter = 220 + alignment_beta = 0.005 + alignment_lambda_vf = 3000.0 + alignment_k = 15 + alignment_partial_robust_level = 100 + alignment_allow_flip = true + alignment_svi_mode = false + alignment_n_sampling = 1000 + alignment_sparse_top_k = 512 + alignment_sparse_calculation_mode = true + alignment_use_chunk = true + alignment_chunk_capacity = 1 + alignment_use_hvg = false + alignment_n_top_genes = 100 + alignment_use_pca = true + alignment_n_pcs = 50 + alignment_max_alignment_cells = 35000 + alignment_seed = 21 + alignment_rbf_neighbors = 64 + alignment_rbf_smoothing = 0.0 + alignment_max_nonrigid_anchors = 5000 + alignment_pytorch_cuda_alloc_conf = null + alignment_conda = "${projectDir}/../environment.alignment.yml" + alignment_bootstrap_dependencies = true + alignment_dynamo_requirement = "dynamo-release @ git+https://github.com/aristoteleo/dynamo-release.git@v1.5.3" + alignment_spateo_requirement = "spateo-release @ git+https://github.com/aristoteleo/spateo-release.git@1bd8a35e6e7b0dfc712d76b1c32987dc06652774" + alignment_anndata_requirement = "anndata>=0.12.10" + alignment_qc_grid_rows = 10 + alignment_qc_grid_cols = 10 + + // Mutually exclusive co-expression rate defaults (enabled by default) + mecr_enabled = true + mecr_neurons_h5ad_path = null + mecr_nonneurons_h5ad_path = null + mecr_cell_metadata_path = null + mecr_taxonomy_metadata_path = null + mecr_cluster_membership_path = null + mecr_taxonomy_level = "CCN202210140_SUPC" + mecr_gene_symbol_column = "gene_symbol" + mecr_target_broad_classes = [ + "Neurons", + "Oligodendrocytes", + "Oligodendrocyte precursors", + "Astrocytes", + "Microglia", + "Fibroblasts", + "Vascular cells", + ] + mecr_marker_min_target_fraction = 0.25 + mecr_marker_max_other_fraction = 0.01 + mecr_normalize_target_sum = 10000.0 + mecr_reference_chunk_rows = 5000 + mecr_wilcoxon_tie_correct = true + mecr_figure_dpi = 180 + mecr_barnyard_top_n_pairs = 6 + mecr_barnyard_max_points = 50000 + mecr_barnyard_random_seed = 0 + mecr_barnyard_log1p = false + + // Scanpy/Squidpy clustering defaults + clustering_squidpy_drop_control_features = true + clustering_squidpy_min_counts = 10 + clustering_squidpy_min_cells = 5 + clustering_squidpy_normalize_target_sum = null + clustering_squidpy_normalize_exclude_highly_expressed = false + clustering_squidpy_normalize_max_fraction = 0.05 + clustering_squidpy_n_pcs = 60 + clustering_squidpy_n_neighbors = 30 + clustering_squidpy_leiden_resolution = 0.5 + clustering_squidpy_umap_min_dist = 0.3 + clustering_squidpy_umap_spread = 1.0 + clustering_squidpy_random_seed = 0 + clustering_squidpy_spatial_point_size = 0.5 + clustering_squidpy_spatial_scatter_point_size = 2.0 + clustering_squidpy_figure_dpi = 180 + clustering_squidpy_use_gpu = false + clustering_squidpy_gpu_conda = "${projectDir}/../environment.clustering-gpu.yml" + clustering_squidpy_gpu_container = "file:///nfsdata/apptainer/merxen_cuda12.6_85107691c6b8660f06f4c841a4d0a24b50d766cb.sif" + clustering_squidpy_gpu_vram_monitor = false + clustering_squidpy_gpu_vram_monitor_interval_seconds = 2 + clustering_squidpy_write_spatialdata_table = true + clustering_squidpy_hierarchical_enabled = true + clustering_squidpy_broad_leiden_resolution = 0.2 + clustering_squidpy_subcluster_leiden_resolution = 0.5 + clustering_squidpy_subcluster_resolution_overrides = [:] + clustering_squidpy_neuron_split_leiden_resolution = 0.15 + clustering_squidpy_neuron_subcluster_leiden_resolution = 0.5 + clustering_squidpy_min_branch_cells = 50 + clustering_squidpy_broad_marker_lookup_path = null + clustering_squidpy_broad_taxonomy_metadata_path = null + clustering_squidpy_broad_cluster_membership_path = null + clustering_squidpy_broad_reference_cache_dir = null + clustering_squidpy_broad_marker_level = "CCN202210140_SUPC" + clustering_squidpy_broad_min_marker_overlap = 3 + clustering_squidpy_broad_max_markers_per_label = 80 + clustering_squidpy_broad_score_margin_threshold = 0.0 + clustering_squidpy_broad_unknown_label = "Mixed/Unknown" + + // Cell-level autocorrelation and transcript-coordinate spatial pattern defaults + spatial_gene_analysis_enabled = true + spatial_gene_analysis_drop_control_features = true + spatial_gene_analysis_min_counts = 0 + spatial_gene_analysis_min_cells = 5 + spatial_gene_analysis_normalize_target_sum = null + spatial_gene_analysis_normalize_exclude_highly_expressed = false + spatial_gene_analysis_normalize_max_fraction = 0.05 + spatial_gene_analysis_n_neighbors = 6 + spatial_gene_analysis_top_n = 10 + spatial_gene_analysis_spatial_point_size = 2.0 + spatial_gene_analysis_figure_dpi = 180 + spatial_gene_analysis_transcript_analysis_enabled = true + spatial_gene_analysis_transcript_min_count = 50 + spatial_gene_analysis_paircorr_min_count = 100 + spatial_gene_analysis_paircorr_max_transcripts_per_gene = 5000 + spatial_gene_analysis_paircorr_distance_edges_um = [0.0, 2.0, 5.0, 20.0, 50.0, 200.0] + spatial_gene_analysis_paircorr_permutations = 100 + spatial_gene_analysis_paircorr_seed = 0 + spatial_gene_analysis_transcript_chunk_rows = 500000 + spatial_gene_analysis_pericellular_distance_um = 5.0 + spatial_gene_analysis_membrane_distance_um = 2.0 + spatial_gene_analysis_signed_distance_edges_um = [-200.0, -50.0, -20.0, -5.0, -2.0, 0.0, 2.0, 5.0, 20.0, 50.0, 200.0] + spatial_gene_analysis_transcript_diagnostic_top_n = 3 + spatial_gene_analysis_transcript_diagnostic_max_genes = 30 + spatial_gene_analysis_transcript_diagnostic_window_um = 250.0 + spatial_gene_analysis_transcript_plot_max_points = 20000 + + // Local MapMyCells defaults (opt in with --stop_stage mapmycells) + mapmycells_reference_mode = "both" + mapmycells_marker_lookup_path = null + mapmycells_precomputed_stats_path = null + mapmycells_region_name = "frontal_a44_a45_a46_a32_acc" + mapmycells_region_labels = ["Human A44-A45", "Human A46", "Human A32", "Human ACC"] + mapmycells_region_cache_dir = null + mapmycells_region_min_cells_per_leaf = 10 + mapmycells_region_force_rebuild = false + mapmycells_region_query_markers_n_per_utility = 10 + mapmycells_drop_level = null + mapmycells_normalization = "raw" + mapmycells_bootstrap_factor = 0.9 + mapmycells_bootstrap_iteration = 100 + mapmycells_n_processors = null + mapmycells_chunk_size = null + mapmycells_rng_seed = 0 + mapmycells_max_gb = null + mapmycells_tmp_dir = null + mapmycells_cloud_safe = false + mapmycells_flatten = false + mapmycells_verbose_csv = false + mapmycells_plots_only = false + mapmycells_query_layer = "counts" + mapmycells_gene_id_column = "ensembl_id" + mapmycells_obs_id_column = null + + // Execution profiles provide machine capacity, concurrency, worker counts, + // local reference paths, and GPU policy. + max_ram_gb = null + warn_ram_gb = null + transcript_chunk_rows = null +} + + +workflow { + failOnIgnore = true +} + +// Process-specific Conda environments (currently ALIGN) must be honored even +// when the pipeline is launched without the optional `conda` profile. Processes +// without a `conda` directive continue to use the caller's active environment. +conda.enabled = true + +process { + //errorStrategy = "ignore" + + withName: "BUILD_SPATIALDATA" { + cpus = 8 + memory = "80 GB" + } + withName: "CELLPOSE_SEGMENT" { + cpus = 12 + memory = "212 GB" + queue = { params.cellpose_gpu ? 'gpu' : 'htc' } + clusterOptions = { params.cellpose_gpu ? '--gpus-per-node=1' : '' } + } + withName: "CELLPOSE_NUCLEI_SEGMENT" { + cpus = 12 + memory = "212 GB" + queue = { params.cellpose_gpu ? 'gpu' : 'htc' } + clusterOptions = { params.cellpose_gpu ? '--gpus-per-node=1' : '' } + } + withName: "PROSEG_SEGMENT" { + cpus = 32 + memory = "220 GB" + queue = 'htc' + } + withName: "ENSURE_PROSEG" { + cpus = 2 + memory = "4 GB" + } + withName: "ENRICH" { + cpus = 8 + memory = "300 GB" + queue = 'dynamic' + } + withName: "VIEWER_CACHE" { + cpus = 8 + memory = "60 GB" + queue = 'dynamic' + } + withName: "MASK_IMAGE_QUANTIFICATION" { + cpus = 4 + memory = "160 GB" + } + withName: "COMPUTE_CORTICAL_DEPTH" { + cpus = 24 + memory = "80 GB" + } + withName: "DISTANCE_FROM_OBJECT_ANNOTATE" { + cpus = 8 + memory = "80 GB" + } + withName: "DISTANCE_FROM_OBJECT_COHORT" { + cpus = params.distance_from_object_n_cpus + memory = "48 GB" + } + withName: "QC" { + cpus = 4 + memory = "24 GB" + } + withName: "ALIGN" { + conda = params.alignment_conda + cpus = 12 + memory = "100 GB" + queue = 'gpu' + clusterOptions = '--gpus-per-node=1' + } + withName: "ALIGN_QC" { + cpus = 4 + memory = "32 GB" + } + withName: "COMPARE" { + cpus = 4 + memory = "32 GB" + } + withName: "VISUALIZE" { + cpus = 4 + memory = "32 GB" + } + withName: "CLUSTERING_SQUIDPY_PREPARE" { + cpus = 4 + memory = "32 GB" + } + withName: "CLUSTERING_SQUIDPY_COMPUTE" { + conda = params.clustering_squidpy_gpu_conda + container = params.clustering_squidpy_gpu_container + cpus = 8 + memory = "32 GB" + queue = 'gpu' + clusterOptions = '--gpus-per-node=1' + } + withName: "CLUSTERING_SQUIDPY_FINALIZE" { + cpus = 4 + memory = "32 GB" + } + withName: "SPATIAL_GENE_ANALYSIS" { + cpus = 12 + memory = "120 GB" + } + withName: "MECR_REFERENCE" { + cpus = 16 + memory = "240 GB" + maxForks = 1 + } + withName: "MECR" { + cpus = 4 + memory = "48 GB" + } + withName: "MAPMYCELLS" { + cpus = 8 + memory = "160 GB" + } +} + +profiles { + // Nextflow selects the reserved `standard` profile when no -profile flag is + // supplied. Keep it equivalent to the explicitly named workstation profile. + standard { + includeConfig 'conf/dwight.config' + } + dwight { + includeConfig 'conf/dwight.config' + } + apptainer { + apptainer.enabled = true + apptainer.autoMounts = true + conda.enabled = false + + // mount directories + process.containerOptions = '-B /data,/nfsdata -C --no-home --home $PWD' + + // set global container + process.container = "file:///nfsdata/apptainer/merxen_cuda12.6_85107691c6b8660f06f4c841a4d0a24b50d766cb.sif" + } + conda { + conda.enabled = true + conda.channels = ['conda-forge', 'bioconda'] + apptainer.enabled = false + process.conda = "${projectDir}/../environment.yml" + } + gpu { + // Keep --nv off CPU-only containers so ProSeg can run on non-GPU nodes. + process { + withName: "CELLPOSE_SEGMENT" { + containerOptions = '-B /data,/nfsdata -C --no-home --home $PWD --nv' + } + withName: "ALIGN" { + containerOptions = '-B /data,/nfsdata -C --no-home --home $PWD --nv' + } + // Retain compatibility with the unsplit clustering process. + withName: "CLUSTERING_SQUIDPY" { + containerOptions = '-B /data,/nfsdata -C --no-home --home $PWD --nv' + } + withName: "CLUSTERING_SQUIDPY_COMPUTE" { + containerOptions = '-B /data,/nfsdata -C --no-home --home $PWD --nv' + } + } + } + + // compute profiles + azure_slurm_hpc { + process.executor = 'slurm' + process.queue = 'htc' + process.time = '24h' + } + local { + includeConfig 'conf/dwight.config' + } +} + +report { + enabled = true + file = "${params.outdir}/nextflow/report.html" + overwrite = true +} + +timeline { + enabled = true + file = "${params.outdir}/nextflow/timeline.html" + overwrite = true +} + +trace { + enabled = true + file = "${params.outdir}/nextflow/trace.tsv" + overwrite = true + sep = '\t' + fields = + 'task_id,hash,native_id,process,tag,name,status,exit,attempt,' + + 'cpus,memory,submit,start,complete,duration,realtime,%cpu,' + + 'peak_rss,peak_vmem,rchar,wchar,read_bytes,write_bytes,workdir' +} diff --git a/merxen_stub/workflows/samplesheet.example.csv b/merxen_stub/workflows/samplesheet.example.csv new file mode 100644 index 0000000..d3fdda7 --- /dev/null +++ b/merxen_stub/workflows/samplesheet.example.csv @@ -0,0 +1,2 @@ +pair_id,analysis_mode,enable_alignment,analysis_segmentation,mecr_enabled,start_stage,stop_stage,only_stage,cortical_depth_enabled,distance_from_object_enabled,distance_from_object_segmentations,merscope_distance_object_annotation_geojson,xenium_distance_object_annotation_geojson,merscope_dir,merscope_spatialdata_path,merscope_image_prefix,merscope_z_range,merscope_transform_path,merscope_channels,xenium_dir,xenium_spatialdata_path,xenium_channels,xenium_min_qv,merscope_voxel_layers,xenium_voxel_layers,xenium_spec_path +P7513,paired,false,both,,,,,,,,,,/media/mathieubo/SSD2/MerXen/P7513/MERSCOPE/region_R2,/media/mathieubo/SSD2/MerXen/P7513/cache/p7513_merscope_source.zarr,202509261108_P7513-2_VMSC19502_region_R2,1-6,,"DAPI,PolyT",/media/mathieubo/SSD2/MerXen/Xenium/,/media/mathieubo/SSD2/MerXen/P7513/cache/p7513_xenium_source.zarr,"DAPI,18S",20,7,2,