Fixes

scripts/create_resources/spatial/process_bruker_cosmx.sh

@@ -14,8 +14,8 @@ cat > /tmp/params.yaml << HERE
 param_list:
 
   - id: "bruker_cosmx/bruker_mouse_brain_cosmx/rep1"
-    input_raw: "https://smi-public.objects.liquidweb.services/HalfBrain.zip"
-    input_flat_files: "https://smi-public.objects.liquidweb.services/Half%20%20Brain%20simple%20%20files%20.zip"
+    input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/HalfBrain.zip"
+    input_flat_files: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Half Brain simple files.zip"
     dataset_name: "Bruker CosMx Mouse Brain"
     dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/cosmx-smi-mouse-brain-ffpe-dataset/"
     dataset_summary: "Bruker CosMx Mouse Brain dataset on FFPE covering a full hemisphere of a mouse brain."
@@ -24,7 +24,7 @@ param_list:
     segmentation_id: ["cell"]
 
   - id: "bruker_cosmx/bruker_human_liver_cosmx"
-    input_raw: "https://smi-public.objects.liquidweb.services/NormalLiverFiles.zip"
+    input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/NormalLiverFiles.zip"
     dataset_name: "Bruker CosMx Human Liver"
     dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/"
     dataset_summary: "Bruker CosMx Human Liver dataset on FFPE."

scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh

@@ -14,8 +14,8 @@ cat > /tmp/params.yaml << HERE
 param_list:
 
   - id: "bruker_cosmx/bruker_mouse_brain_cosmx/rep1"
-    input_raw: "https://smi-public.objects.liquidweb.services/HalfBrain.zip"
-    input_flat_files: "https://smi-public.objects.liquidweb.services/Half%20%20Brain%20simple%20%20files%20.zip"
+    input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/HalfBrain.zip"
+    input_flat_files: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Half Brain simple files.zip"
     dataset_name: "Bruker CosMx Mouse Brain"
     dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/cosmx-smi-mouse-brain-ffpe-dataset/"
     dataset_summary: "Bruker CosMx Mouse Brain dataset on FFPE covering a full hemisphere of a mouse brain."
@@ -24,7 +24,7 @@ param_list:
     segmentation_id: ["cell"]
 
   - id: "bruker_cosmx/bruker_human_liver_cosmx"
-    input_raw: "https://smi-public.objects.liquidweb.services/NormalLiverFiles.zip"
+    input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/NormalLiverFiles.zip"
     dataset_name: "Bruker CosMx Human Liver"
     dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/"
     dataset_summary: "Bruker CosMx Human Liver dataset on FFPE."

src/methods_cell_type_annotation/moscot/script.py

@@ -75,13 +75,13 @@
 adata_sp.X = adata_sp.layers["normalized"]
 adata_sp.obsm["spatial"] = adata_sp.obs[["centroid_x", "centroid_y"]].to_numpy()
 
-sc.pp.pca(adata_sc, n_comps=50)   # X is the normalized layer set above
+sc.pp.pca(adata_sc, n_comps=30)   # X is the normalized layer set above
 
 # Define mapping problem
 mp = MappingProblem(adata_sc=adata_sc, adata_sp=adata_sp)
 
 mp = mp.prepare(
-    sc_attr={"attr": "obsm", "key": "X_pca"},   # <-- 50-dim, not raw genes
+    sc_attr={"attr": "obsm", "key": "X_pca"},   # <-- 30-dim, not raw genes
     xy_callback="local-pca",
 )
 

src/methods_transcript_assignment/comseg/config.vsh.yaml

@@ -63,6 +63,30 @@ arguments:
     type: boolean
     default: true
     description: "Allow disconnected polygons in segmentation"
+  - name: --n_workers
+    type: integer
+    default: 4
+    description: |
+      Maximum number of dask worker processes used to run ComSeg patches in
+      parallel. ComSeg builds an in-memory graph per patch, so peak memory
+      scales with the number of concurrent workers; the value is clamped to
+      (cpus - 1). Lower this if the job is OOM-killed (exit 137); raise it
+      (up to the CPU count) to run faster when memory allows.
+  - name: --worker_memory_limit
+    type: string
+    default: "0"
+    description: |
+      Per-worker memory limit for the dask LocalCluster. Default "0" (also
+      "none"/"off"/"disabled") neutralises dask's per-worker memory monitor by
+      setting each worker's limit to the FULL detected container memory. ComSeg's
+      per-patch memory is dominated by unmanaged native allocations
+      (numba/shapely/graphs) that dask cannot spill, so a tighter limit does not
+      buy spilling — it only makes the nanny kill otherwise-healthy workers
+      (KilledWorker) once they cross ~95% of the limit. Total memory is instead
+      bounded by n_workers, and a genuine OOM is caught by the container cgroup +
+      Nextflow retry. Set an explicit value (e.g. "20GB") only if you want dask
+      to police each worker; "auto" is discouraged because it divides detected
+      RAM by CPU count and is easily smaller than a single patch's working set.
 
 
 resources:

src/methods_transcript_assignment/comseg/script.py

@@ -4,6 +4,7 @@
 import anndata as ad
 import pandas as pd
 import numpy as np
+import dask
 import dask.dataframe as dd
 
 ## VIASH START
@@ -24,6 +25,8 @@
     "gene_column": "feature_name",
     "norm_vector": False,
     "allow_disconnected_polygon": True,
+    "n_workers": 4,
+    "worker_memory_limit": "0",
 }
 meta = {
     "name": "comseg",
@@ -79,12 +82,54 @@
 
 # ComSeg processes each transcript patch independently and is pure-Python, so it
 # runs single-threaded unless a parallelization backend is enabled. Use the dask
-# backend to spread patches across the allocated CPUs (see midcpu/highmem labels).
-n_workers = max((meta["cpus"] or os.cpu_count() or 1) - 1, 1)
+# backend to spread patches across workers, but CAP concurrency: comseg builds an
+# in-memory graph per patch, so peak memory scales with the number of concurrent
+# workers. Running one worker per CPU (14 on midcpu) multiplied peak RAM ~14x and
+# got the job OOM-killed (exit 137 = cgroup SIGKILL). Bounding the worker count is
+# the real memory safeguard (the per-worker memory monitor is neutralised below,
+# see the note there).
+cpu_cap = max((meta["cpus"] or os.cpu_count() or 1) - 1, 1)
+n_workers = max(min(par["n_workers"], cpu_cap), 1)
+
+# `distributed` defaults worker startup to "spawn", which re-imports this script
+# in every worker. As a plain viash script it has no `if __name__ == "__main__"`
+# guard, so the re-import re-runs the module-level code (creating yet another
+# cluster) and Python's spawn bootstrap check aborts with "An attempt has been
+# made to start a new process before ... bootstrapping". Force "fork" so workers
+# inherit the already-initialised interpreter instead of re-importing it; this is
+# the Linux-native default and the mode comseg ran under before the deps bumped.
+dask.config.set({"distributed.worker.multiprocessing-method": "fork"})
+
 sopa.settings.parallelization_backend = "dask"
 sopa.settings.dask_client_kwargs["n_workers"] = n_workers
 sopa.settings.dask_client_kwargs["threads_per_worker"] = 1  # CPU-bound work, avoid GIL contention
-print(f"Running ComSeg with dask backend, n_workers={n_workers}", flush=True)
+
+# ComSeg's per-patch memory is mostly UNMANAGED (native numba/shapely/graph
+# allocations, plus the forked interpreter's copy-on-write pages), which dask
+# cannot spill. A per-worker memory_limit therefore never spills here; it only
+# lets the nanny kill workers that cross ~95% of the limit -> KilledWorker, even
+# though nothing is leaking. dask's "auto" is worse: it divides detected RAM by
+# the CPU count (not n_workers), so on an 8-core box each worker gets ~1/8 of RAM
+# and a normal ~1 GB patch already trips 95%. We therefore neutralise the monitor
+# by default, but sopa reads worker.memory_manager.memory_limit and compares it to
+# 4 GiB, so the value must be a real number (0/None makes sopa crash). Setting the
+# limit to the FULL detected container memory does both: sopa sees a number, and
+# the nanny only ever fires near 95% of the whole cgroup — where the OS + Nextflow
+# retry already take over — never prematurely for a single patch. Total memory is
+# bounded by n_workers instead. An explicit value (e.g. "20GB") is honoured as-is.
+_mem_limit_arg = (par["worker_memory_limit"] or "").strip().lower()
+if _mem_limit_arg in ("", "0", "none", "off", "disabled"):
+    from distributed.system import MEMORY_LIMIT as _CONTAINER_MEMORY
+    worker_memory_limit = _CONTAINER_MEMORY
+else:
+    worker_memory_limit = par["worker_memory_limit"]
+sopa.settings.dask_client_kwargs["memory_limit"] = worker_memory_limit
+
+print(
+    f"Running ComSeg with dask backend, n_workers={n_workers} "
+    f"(cap {cpu_cap}), memory_limit={worker_memory_limit!r}",
+    flush=True,
+)
 
 sopa.segmentation.comseg(sdata, config)
 

src/methods_transcript_assignment/fastreseg/config.vsh.yaml

@@ -46,24 +46,45 @@ engines:
           - gdal-bin
           - libgdal-dev
           - r-base
-      - type: r
-        packages:
-          - devtools
-          - terra
-          - Matrix
-          - dbscan
-          - igraph
-          - matrixStats
-          - codetools
-          - data.table
-          # Add other CRAN packages here
-        github:
-          # Add GitHub packages here if needed
-          - Nanostring-Biostats/FastReseg
-      #- type: python
-      #  github:
-      #    - openproblems-bio/core#subdirectory=packages/python/openproblems
-    __merge__: 
+          # FastReseg's R dependencies as prebuilt Debian binaries. Installing
+          # these via apt avoids compiling terra/igraph/sf/spatstat/stringi/etc.
+          # from source (extremely slow, especially under amd64 emulation on
+          # arm64 hosts). Only concaveman + GiottoUtils/GiottoClass + FastReseg
+          # itself still build from source below; their heavy deps are already
+          # present as binaries here.
+          - r-cran-remotes
+          - r-cran-terra
+          - r-cran-igraph
+          - r-cran-sf
+          - r-cran-matrix
+          - r-cran-mass
+          - r-cran-e1071
+          - r-cran-dbscan
+          - r-cran-data.table
+          - r-cran-dplyr
+          - r-cran-ggplot2
+          - r-cran-reshape2
+          - r-cran-stringr
+          - r-cran-fs
+          - r-cran-lmtest
+          - r-cran-matrixstats
+          - r-cran-geometry
+          - r-cran-spatstat.geom
+      - type: docker
+        # Install FastReseg from GitHub. upgrade = "never" so remotes reuses the
+        # apt-installed binary deps above instead of re-fetching/recompiling them
+        # from CRAN; only the unpackaged deps (concaveman, GiottoUtils,
+        # GiottoClass) and FastReseg itself are compiled from source.
+        run:
+          - Rscript -e 'options(warn = 2); remotes::install_github("Nanostring-Biostats/FastReseg", upgrade = "never", repos = "https://cran.rstudio.com")'
+      - type: python
+        # `plankton` (pulled in by txsim.run_ssam for ssam cell annotation in
+        # input.py) still does `from matplotlib.cm import get_cmap`, removed in
+        # matplotlib 3.9. Pin < 3.9 so run_ssam imports. This is a local setup
+        # step, so it runs AFTER the __merge__'d spatialdata/txsim partials (which
+        # otherwise pull matplotlib >= 3.9) and wins. Mirrors methods_cell_type_annotation/ssam.
+        pypi: [planktonspace, "matplotlib<3.9"]
+    __merge__:
       - /src/base/setup_txsim_partial.yaml
       - /src/base/setup_spatialdata_partial.yaml
   - type: native

src/methods_transcript_assignment/fastreseg/input.py

@@ -7,6 +7,61 @@
 import txsim as tx
 import anndata as ad
 import argparse
+from scipy.sparse import csr_matrix
+
+
+def generate_adata(input_spots, cell_id_col="cell", gene_col="Gene"):
+    """Aggregate per-transcript assignments into a cell x gene AnnData.
+
+    Reimplementation of txsim.preprocessing.generate_adata: txsim's version is
+    incompatible with anndata>=0.12 (it passes the removed `dtype=` argument to
+    AnnData() and uses per-row item assignment `adata[cell, :] = ...`, which
+    anndata no longer supports). This fills the count matrix directly and
+    produces an identical output structure (X, layers['raw_counts'],
+    obs/var stats, uns['spots'/'pct_noise']). Kept in sync with the copy in
+    src/methods_count_aggregation/basic_count_aggregation/script.py.
+    """
+    spots = input_spots.copy()
+    pct_noise = sum(spots[cell_id_col] <= 0) / len(spots[cell_id_col])
+    spots_raw = spots.copy()  # kept in uns['spots']; 0 (background) -> None
+    spots_raw.loc[spots_raw[cell_id_col] == 0, cell_id_col] = None
+    spots = spots[spots[cell_id_col] > 0]
+
+    cell_ids = pd.unique(spots[cell_id_col])
+    genes = pd.unique(spots[gene_col])
+    cell_pos = {c: i for i, c in enumerate(cell_ids)}
+
+    # Populate the count matrix + centroids per cell (no AnnData item assignment).
+    # feature_name may be categorical, so value_counts() can return unobserved
+    # categories; reindex to `genes` (fill 0) to align to the var order, mirroring
+    # txsim's `.reindex(var_names, fill_value=0)`.
+    X = np.zeros((len(cell_ids), len(genes)), dtype=np.float32)
+    centroid_x = np.zeros(len(cell_ids))
+    centroid_y = np.zeros(len(cell_ids))
+    for cell_id, grp in spots.groupby(cell_id_col, sort=False):
+        row = cell_pos[cell_id]
+        cts = grp[gene_col].value_counts().reindex(genes, fill_value=0)
+        X[row, :] = cts.values
+        centroid_x[row] = grp["x"].mean()
+        centroid_y[row] = grp["y"].mean()
+
+    adata = ad.AnnData(csr_matrix(X))
+    adata.obs["cell_id"] = cell_ids
+    adata.obs_names = [f"{i:d}" for i in cell_ids]
+    adata.var_names = genes
+    adata.obs["centroid_x"] = centroid_x
+    adata.obs["centroid_y"] = centroid_y
+
+    adata.uns["spots"] = spots_raw
+    adata.uns["pct_noise"] = pct_noise
+    adata.layers["raw_counts"] = adata.X.copy()
+
+    adata.obs["n_counts"] = np.ravel(adata.layers["raw_counts"].sum(axis=1))
+    adata.obs["n_genes"] = adata.layers["raw_counts"].getnnz(axis=1)
+    adata.var["n_counts"] = np.ravel(adata.layers["raw_counts"].sum(axis=0))
+    adata.var["n_cells"] = adata.layers["raw_counts"].getnnz(axis=0)
+    return adata
+
 
 def parse_arguments():
     parser = argparse.ArgumentParser(
@@ -70,8 +125,12 @@ def parse_arguments():
 print('Transforming transcripts coordinates', flush=True)
 transcripts = sd.transform(sdata['transcripts'], to_coordinate_system='global')
 
-transcripts_df = transcripts.compute()
-transcripts_df.rename(columns = {'feature_name': 'target', 
+# .copy() is required: for a single-partition points object, transcripts.compute()
+# returns a reference to the dask graph's backing frame, so an in-place rename here
+# would corrupt it and the later transcripts.compute() (passed to run_ssam) would
+# yield renamed columns that no longer match the dask _meta -> "Metadata mismatch".
+transcripts_df = transcripts.compute().copy()
+transcripts_df.rename(columns = {'feature_name': 'target',
 'transcript_id': 'UMI_transID', 'cell_id': 'UMI_cellID'}, inplace = True)
 
 transcripts_df = transcripts_df.loc[:, ['target', 'x', 'y', 'z', 'UMI_transID', 'UMI_cellID']]
@@ -82,7 +141,7 @@ def parse_arguments():
 df = sdata['transcripts'].compute()
 df.feature_name = df.feature_name.astype(str)
 
-adata_sp = tx.preprocessing.generate_adata(df, cell_id_col='cell_id', gene_col='feature_name') #TODO: x and y refers to a specific coordinate system. Decide which space we want to use here. (probably should be handled in the previous assignment step)
+adata_sp = generate_adata(df, cell_id_col='cell_id', gene_col='feature_name') #TODO: x and y refers to a specific coordinate system. Decide which space we want to use here. (probably should be handled in the previous assignment step)
 adata_sp.layers['counts'] = adata_sp.layers['raw_counts']
 del adata_sp.layers['raw_counts']
 adata_sp.var["gene_name"] = adata_sp.var_names

src/methods_transcript_assignment/fastreseg/orchestrator.sh

@@ -1,5 +1,10 @@
 #!/bin/bash
 
+# Fail loudly: without this, a crash in any sub-step (input.py / script.R /
+# output.py) is swallowed and Viash only reports the generic "required output
+# file is missing", hiding the real Python/R traceback.
+set -eo pipefail
+
 ## VIASH START
 # The following code has been auto-generated by Viash.
 par_input_ist='resources_test/task_ist_preprocessing/mouse_brain_combined/raw_ist.zarr'

src/methods_transcript_assignment/fastreseg/output.py

@@ -36,6 +36,12 @@ def convert_to_lower_dtype(arr):
 ##converting to dask transcript output for spatialdata format
 df = df.loc[:,["x", "y", "z", "target", "updated_cellID", "updated_celltype", "UMI_transID"]]
 df.rename(columns={"updated_cellID": "cell_id", "target": "feature_name", "UMI_transID": "transcript_id"}, inplace=True)
+# Match the raw_ist convention where feature_name is categorical. Rebuilt from
+# R's CSV it would otherwise be a pyarrow-backed nullable string that survives
+# the zarr round-trip; the downstream count-aggregation then builds var_names
+# from it and write_h5ad refuses to serialize the nullable-string index
+# ("anndata.settings.allow_write_nullable_strings is None").
+df["feature_name"] = df["feature_name"].astype(str).astype("category")
 transcripts_dask = dask.dataframe.from_pandas(df, npartitions = 1)
 
 

src/methods_transcript_assignment/segger/script.py

@@ -34,6 +34,23 @@
 }
 ## VIASH END
 
+# The Nextflow/k8s GPU task hands us an EMPTY CUDA_VISIBLE_DEVICES. An empty (or
+# whitespace) value masks every GPU from the CUDA *driver* API, so cudf/cuspatial and
+# numba enumerate zero devices and segger dies deep in tiling with
+# "cudaErrorNoDevice: no CUDA-capable device is detected" / numba's
+# "IndexError: list index out of range" (self.gpus[devnum] on an empty device list).
+# torch's is_available() is NVML-based and ignores CVD, so the gate below still passes
+# and the failure only surfaces inside the `segger segment` subprocess. Drop an empty
+# value so the allocated device (exposed via NVIDIA_VISIBLE_DEVICES) is visible again;
+# run_segger's os.environ.copy() then inherits the fix for the subprocess.
+if os.environ.get("CUDA_VISIBLE_DEVICES", "x").strip() == "":
+    print(
+        "CUDA_VISIBLE_DEVICES was empty; unsetting so the allocated GPU is visible "
+        "to cudf/numba (torch's NVML check hides this).",
+        flush=True,
+    )
+    del os.environ["CUDA_VISIBLE_DEVICES"]
+
 # segger runs its tiling / GNN training / prediction end-to-end on the GPU
 # (cudf, cuspatial, torch kernels). Fail fast with a clear message rather than
 # deep inside the subprocess when no CUDA device is present.