From d68d39fac71e756c8fbf4db79014bf6e5472a078 Mon Sep 17 00:00:00 2001 From: NazimMechkouri Date: Wed, 26 Aug 2026 11:54:48 +0200 Subject: [PATCH 01/38] Add geojson processing rule and script --- config/config.yaml | 4 ++ workflow/Snakefile | 14 ++++++ workflow/rules/common.smk | 1 + workflow/rules/process_geojson.smk | 23 +++++++++ workflow/scripts/process_geojson.py | 78 +++++++++++++++++++++++++++++ 5 files changed, 120 insertions(+) create mode 100644 workflow/rules/process_geojson.smk create mode 100644 workflow/scripts/process_geojson.py diff --git a/config/config.yaml b/config/config.yaml index 6bc2735..9ebef73 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -313,6 +313,10 @@ resources: mem_mb: 32768 runtime: 30 threads: 2 + process_geojson: + mem_mb: 4096 + runtime: 30 + threads: 2 qc_sweep: mem_mb: 49152 runtime: 180 diff --git a/workflow/Snakefile b/workflow/Snakefile index 9df8b03..714df4e 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -7,6 +7,7 @@ postprocessing → exploration). Rules are split across rules/*.smk and included directly below; helper functions live in rules/common.smk. """ +import glob import os import sys import pandas as pd @@ -132,6 +133,17 @@ if CORE_SAMPLES is None: ) SAMPLE_IDS = list(SAMPLES.keys()) +# Optional per-sample region GeoJSON (ROI-based analysis). Presence is optional per +# sample; a geojson-like file that matches neither accepted name is a fatal config +# error (see _find_geojson). +SAMPLES_WITH_GEOJSON = [s for s in SAMPLE_IDS if _find_geojson(s) is not None] +for _s in SAMPLE_IDS: + if _s not in SAMPLES_WITH_GEOJSON: + print( + f"[process_geojson] No GeoJSON found for sample '{_s}' — " + "no ROI-based analysis will be performed for this sample." + ) + # Optional per-column colour palettes (sample + design columns) for downstream visualization SAMPLE_COLORS = config.get("sample_colors", {}) or {} @@ -256,6 +268,7 @@ elif MODE == "merscope": include: "rules/validate_input.smk" include: "rules/qc_sweep.smk" include: "rules/preprocess_umap.smk" +include: "rules/process_geojson.smk" include: "rules/leiden_analysis.smk" include: "rules/generate_annotation_template.smk" include: "rules/ingest_ref.smk" @@ -311,6 +324,7 @@ rule run_preprocessing: expand(rules.preprocess_umap.output.adata, sample=SAMPLE_IDS), expand(rules.preprocess_umap.output.metadata, sample=SAMPLE_IDS), expand(rules.preprocess_umap.output.report, sample=SAMPLE_IDS), + expand(rules.process_geojson.output.geojson, sample=SAMPLES_WITH_GEOJSON), *QUPATH_IMAGES, *( expand( diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 2eb1183..593ac34 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -310,6 +310,7 @@ def get_all_targets(wildcards): targets += expand(rules.preprocess_umap.output.adata, sample=SAMPLE_IDS) targets += expand(rules.preprocess_umap.output.metadata, sample=SAMPLE_IDS) targets += expand(rules.preprocess_umap.output.report, sample=SAMPLE_IDS) + targets += expand(rules.process_geojson.output.geojson, sample=SAMPLES_WITH_GEOJSON) targets += list(QUPATH_IMAGES) # head QuPath image(s) for the active mode (or []) if RUN_LEIDEN_ANALYSIS: targets += expand( diff --git a/workflow/rules/process_geojson.smk b/workflow/rules/process_geojson.smk new file mode 100644 index 0000000..33509d5 --- /dev/null +++ b/workflow/rules/process_geojson.smk @@ -0,0 +1,23 @@ +rule process_geojson: + """Load a sample's OPTIONAL region GeoJSON ({sample}_tissue_hires_image.geojson + or {sample}_morphology.geojson) and write it into the sample dir. Only invoked + for samples in SAMPLES_WITH_GEOJSON (see common.smk / _find_geojson).""" + input: + geojson=lambda wc: _find_geojson(wc.sample), + output: + geojson=f"{SAMPLES_DIR}/{{sample}}/{{sample}}_annotations.geojson", + log: + out=f"{LOGDIR}/process_geojson/{{sample}}.out", + err=f"{LOGDIR}/process_geojson/{{sample}}.err", + benchmark: + f"{LOGDIR}/benchmarks/process_geojson/{{sample}}.tsv" + conda: + "../envs/visiumhd.yaml" + threads: get_resource("process_geojson", "threads") + resources: + mem_mb=mem_mb_attempt("process_geojson"), + runtime=get_resource("process_geojson", "runtime"), + params: + sample_id=lambda wc: wc.sample, + script: + "../scripts/process_geojson.py" \ No newline at end of file diff --git a/workflow/scripts/process_geojson.py b/workflow/scripts/process_geojson.py new file mode 100644 index 0000000..a53d8f5 --- /dev/null +++ b/workflow/scripts/process_geojson.py @@ -0,0 +1,78 @@ +""" +process_geojson.py – CORE: load, validate, and stage a sample's region GeoJSON. +================================================================================ +Reads the OPTIONAL region-annotation GeoJSON for a sample (either +{sample}_tissue_hires_image.geojson or {sample}_morphology.geojson, resolved +upstream by _find_geojson in common.smk — this script only ever sees a single, +already-resolved path, since presence/naming is validated before the rule is even +scheduled). Sanity-checks the geometries, tags each feature with the sample ID, +and writes it into the sample's output directory for downstream ROI-based +analyses to consume. +""" +import logging +import sys +import traceback +from pathlib import Path + +import geopandas as gpd + +# ── Logging ────────────────────────────────────────────────────────────────── +log_handlers = [logging.StreamHandler(sys.stderr)] +if hasattr(snakemake, "log"): + if snakemake.log.out: + Path(snakemake.log.out).parent.mkdir(parents=True, exist_ok=True) + log_handlers.append(logging.FileHandler(snakemake.log.out, mode="w")) + if snakemake.log.err: + Path(snakemake.log.err).parent.mkdir(parents=True, exist_ok=True) + sys.stderr = open(snakemake.log.err, "w") +logging.basicConfig(level=logging.INFO, + format="%(asctime)s %(levelname)-8s %(message)s", + datefmt="%Y-%m-%d %H:%M:%S", handlers=log_handlers) +log = logging.getLogger("process_geojson") + +# ── Parameters ─────────────────────────────────────────────────────────────── +in_geojson = str(snakemake.input.geojson) +sample_id = snakemake.params.sample_id +out_geojson = str(snakemake.output.geojson) + +try: + log.info("=" * 70) + log.info("Processing region GeoJSON for sample: %s", sample_id) + log.info(" input: %s", in_geojson) + log.info("=" * 70) + + # ── 1. Load ──────────────────────────────────────────────────────── + gdf = gpd.read_file(in_geojson) + n_total = len(gdf) + log.info("Loaded %d feature(s)", n_total) + + if n_total == 0: + raise ValueError(f"'{in_geojson}' contains no features.") + + # ── 2. Validate / repair geometries ─────────────────────────────────── + invalid_mask = ~gdf.geometry.is_valid + n_invalid = int(invalid_mask.sum()) + if n_invalid: + log.warning("%d/%d invalid geometries found — attempting buffer(0) repair", + n_invalid, n_total) + gdf.loc[invalid_mask, "geometry"] = gdf.loc[invalid_mask, "geometry"].buffer(0) + still_invalid = int((~gdf.geometry.is_valid).sum()) + if still_invalid: + raise ValueError( + f"{still_invalid}/{n_total} geometries remain invalid after repair " + f"— inspect '{in_geojson}'.") + + # ── 3. Tag with sample of origin ───────────────────────────────────── + # Survives later concatenation across samples without losing provenance. + gdf["sample_id"] = sample_id + + # ── 4. Save ──────────────────────────────────────────────────────── + Path(out_geojson).parent.mkdir(parents=True, exist_ok=True) + log.info("Saving %d feature(s) → %s", n_total, out_geojson) + gdf.to_file(out_geojson, driver="GeoJSON") + + log.info("GeoJSON processing complete for %s.", sample_id) + +except Exception: + log.error("FAILED for %s:\n%s", sample_id, traceback.format_exc()) + raise \ No newline at end of file From 5fffd7e97a737c0dedba3966ac44ff1d2f47588c Mon Sep 17 00:00:00 2001 From: NazimMechkouri Date: Fri, 28 Aug 2026 10:14:51 +0200 Subject: [PATCH 02/38] Apply snakefmt formatting and fix trailing newline --- workflow/rules/process_geojson.smk | 2 +- workflow/scripts/process_geojson.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/rules/process_geojson.smk b/workflow/rules/process_geojson.smk index 33509d5..7d1d3cb 100644 --- a/workflow/rules/process_geojson.smk +++ b/workflow/rules/process_geojson.smk @@ -20,4 +20,4 @@ rule process_geojson: params: sample_id=lambda wc: wc.sample, script: - "../scripts/process_geojson.py" \ No newline at end of file + "../scripts/process_geojson.py" diff --git a/workflow/scripts/process_geojson.py b/workflow/scripts/process_geojson.py index a53d8f5..ea5d104 100644 --- a/workflow/scripts/process_geojson.py +++ b/workflow/scripts/process_geojson.py @@ -75,4 +75,4 @@ except Exception: log.error("FAILED for %s:\n%s", sample_id, traceback.format_exc()) - raise \ No newline at end of file + raise From ca894d85d6f26aa7e5332f9d6c3b9b7019022461 Mon Sep 17 00:00:00 2001 From: NazimMechkouri Date: Fri, 28 Aug 2026 10:53:02 +0200 Subject: [PATCH 03/38] Fix: add missing _find_geojson helper in common.smk --- workflow/rules/common.smk | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 593ac34..b85fc40 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -206,6 +206,17 @@ def _ate_prepare_inputs(wildcards): return inputs +def _find_geojson(sample): + """Locate a sample's OPTIONAL region GeoJSON in GEOJ_DIR, trying known filename + patterns in order. Returns the first path that exists, or None when neither is + present (used both to gate SAMPLES_WITH_GEOJSON and as process_geojson's input).""" + for _suffix in ("_tissue_hires_image.geojson", "_morphology.geojson"): + _gj = os.path.join(GEOJ_DIR, f"{sample}{_suffix}") + if os.path.isfile(_gj): + return _gj + return None + + # ── HEAD (MERSCOPE) input helper ───────────────────────────────────────────── def merscope_dir_for(sample): """MERSCOPE region directory for a sample, from the merscope.merscope_dir From 1600a1b448d18c935cf8ff5d98382787b77f1c63 Mon Sep 17 00:00:00 2001 From: NazimMechkouri Date: Tue, 1 Sep 2026 23:05:18 +0200 Subject: [PATCH 04/38] Add pseudobulk aggregate and DE rules --- workflow/rules/pseudobulk_aggregate.smk | 2 +- workflow/rules/pseudobulk_de.smk | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/rules/pseudobulk_aggregate.smk b/workflow/rules/pseudobulk_aggregate.smk index 312f2d0..fe24ebf 100644 --- a/workflow/rules/pseudobulk_aggregate.smk +++ b/workflow/rules/pseudobulk_aggregate.smk @@ -21,7 +21,7 @@ rule pseudobulk_aggregate: benchmark: f"{LOGDIR}/benchmarks/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.tsv" wildcard_constraints: - annot_type="tsv_annotation|ingest_annotation", + annot_type="tsv_annotation|ingest_annotation|external_annotation", #remove pipe analysis_level="by_region|by_celltype_region|by_niche_region", conda: "../envs/pseudobulk_aggregate.yaml" diff --git a/workflow/rules/pseudobulk_de.smk b/workflow/rules/pseudobulk_de.smk index 15c5536..e554e3c 100644 --- a/workflow/rules/pseudobulk_de.smk +++ b/workflow/rules/pseudobulk_de.smk @@ -15,7 +15,7 @@ rule pseudobulk_de: benchmark: f"{LOGDIR}/benchmarks/pseudobulk_de/{{annot_type}}_{{analysis_level}}.tsv" wildcard_constraints: - annot_type="tsv_annotation|refined_annotation|ingest_annotation", + annot_type="tsv_annotation|refined_annotation|ingest_annotation|external_annotation", analysis_level="by_region|by_celltype_region|by_niche_region", conda: "../envs/pseudobulk_de.yaml" From 3d23c5991c0bd0b938710d305848ddc0917b9326 Mon Sep 17 00:00:00 2001 From: Sergio Oterino Sogo <68707367+Sergio-ote@users.noreply.github.com> Date: Fri, 4 Sep 2026 15:28:03 +0200 Subject: [PATCH 05/38] Fix concatenation of AnnData objects with sample_batch --- workflow/scripts/spatial_niches.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/spatial_niches.py b/workflow/scripts/spatial_niches.py index 19b3ae7..2143aa3 100644 --- a/workflow/scripts/spatial_niches.py +++ b/workflow/scripts/spatial_niches.py @@ -268,7 +268,7 @@ def _resolve_precomputed(): # preprocess_umap already sets obs["sample"] = sample_id per cell, so we do # NOT pass label= (which would clash with that column). keys + index_unique # only disambiguate obs_names across samples. - adata = sc.concat(adatas, join="inner", keys=sample_ids, index_unique="-") + adata = sc.concat(adatas, join="inner", label="sample_batch", keys=sample_ids) del adatas log.info("Concatenated: %d cells, %d genes", adata.n_obs, adata.n_vars) From e1ef7301e60c9c597b710cf4cf8060b6bb511a19 Mon Sep 17 00:00:00 2001 From: NazimMechkouri Date: Tue, 8 Sep 2026 21:37:34 +0200 Subject: [PATCH 06/38] scripts: add per-sample external-metadata normaliser (called by prepare_input for x5k and ate; backs up the user's metadata and reformats to the required 2-column shape) --- .gitignore | 1 + .../scripts/normalize_external_metadata.py | 130 ++++++++++++++++++ 2 files changed, 131 insertions(+) create mode 100644 workflow/scripts/normalize_external_metadata.py diff --git a/.gitignore b/.gitignore index 468aab2..f4362d0 100644 --- a/.gitignore +++ b/.gitignore @@ -43,3 +43,4 @@ config/visiumhd_samples_copy.csv # Results (these live on scratch, not in the repo) # Adjust if you prefer to track them +meta.tsv diff --git a/workflow/scripts/normalize_external_metadata.py b/workflow/scripts/normalize_external_metadata.py new file mode 100644 index 0000000..0a7fdda --- /dev/null +++ b/workflow/scripts/normalize_external_metadata.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +""" +normalize_external_metadata.py — normalise SpaceBlocks external-annotation metadata. + +Rewrites metadata_{sample}.tsv into the strict 2-column shape annotate_cells expects: + + barcode +