new rule : Process geojson

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

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(

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
@@ -310,6 +321,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(

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"

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