From f80b26160f3f2f9a914e9659159b1d993f5427bc Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Sun, 19 Jul 2026 23:22:50 +0200
Subject: [PATCH 01/19] Register cellposev4 in benchmark run scripts

The cellposev4 segmentation component exists but was not listed in any
run script, so it never ran. Add it to segmentation_methods alongside
cellpose (active in the seqeracloud scripts, commented in the local/test
scripts, matching the existing convention).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 scripts/run_benchmark/run_full_local.sh       | 1 +
 scripts/run_benchmark/run_full_seqeracloud.sh | 1 +
 scripts/run_benchmark/run_test_local.sh       | 1 +
 scripts/run_benchmark/run_test_seqeracloud.sh | 1 +
 4 files changed, 4 insertions(+)

diff --git a/scripts/run_benchmark/run_full_local.sh b/scripts/run_benchmark/run_full_local.sh
index db9ba0457..b6fe94888 100755
--- a/scripts/run_benchmark/run_full_local.sh
+++ b/scripts/run_benchmark/run_full_local.sh
@@ -33,6 +33,7 @@ default_methods:
 segmentation_methods:
   - custom_segmentation
   # - cellpose
+  # - cellposev4
   - binning
   # - stardist
   # - watershed
diff --git a/scripts/run_benchmark/run_full_seqeracloud.sh b/scripts/run_benchmark/run_full_seqeracloud.sh
index 65d353b09..ae3d5a3f3 100755
--- a/scripts/run_benchmark/run_full_seqeracloud.sh
+++ b/scripts/run_benchmark/run_full_seqeracloud.sh
@@ -25,6 +25,7 @@ default_methods:
 segmentation_methods:
   - custom_segmentation
   - cellpose
+  - cellposev4
   - binning
   - stardist
   - watershed
diff --git a/scripts/run_benchmark/run_test_local.sh b/scripts/run_benchmark/run_test_local.sh
index 75ba1ac4a..b20058497 100755
--- a/scripts/run_benchmark/run_test_local.sh
+++ b/scripts/run_benchmark/run_test_local.sh
@@ -28,6 +28,7 @@ default_methods:
 segmentation_methods:
   - custom_segmentation
   # - cellpose
+  # - cellposev4
   - binning
   # - stardist
   # - watershed
diff --git a/scripts/run_benchmark/run_test_seqeracloud.sh b/scripts/run_benchmark/run_test_seqeracloud.sh
index 86bef47e8..1a3124961 100755
--- a/scripts/run_benchmark/run_test_seqeracloud.sh
+++ b/scripts/run_benchmark/run_test_seqeracloud.sh
@@ -24,6 +24,7 @@ default_methods:
 segmentation_methods:
   - custom_segmentation
   - cellpose
+  - cellposev4
   - binning
   - stardist
   - watershed

From 1a2fa097b73295ad2f88671906ce53605403f1b3 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 00:11:37 +0200
Subject: [PATCH 02/19] fix anndata version mismatch with txsim

---
 .../basic_count_aggregation/script.py         | 60 ++++++++++++++++++-
 1 file changed, 58 insertions(+), 2 deletions(-)

diff --git a/src/methods_count_aggregation/basic_count_aggregation/script.py b/src/methods_count_aggregation/basic_count_aggregation/script.py
index 5be4265e9..a1ac203f8 100644
--- a/src/methods_count_aggregation/basic_count_aggregation/script.py
+++ b/src/methods_count_aggregation/basic_count_aggregation/script.py
@@ -1,5 +1,61 @@
+import numpy as np
+import pandas as pd
+import anndata as ad
 import spatialdata as sd
-import txsim as tx
+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']).
+    """
+    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
+
 
 ## VIASH START
 par = {
@@ -14,7 +70,7 @@
 sdata = sd.read_zarr(par['input'])
 df = sdata['transcripts'].compute() # TODO: Could optimize tx.preprocessing.generate_adata to work on spatialdata
 
-adata = 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 = 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.layers['counts'] = adata.layers['raw_counts']
 del adata.layers['raw_counts']
 adata.var["gene_name"] = adata.var_names

From 82add803b17430c35df08b6b0222b76783bb1bb2 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 00:29:37 +0200
Subject: [PATCH 03/19] add segger to workflow (test)

---
 src/workflows/run_benchmark/config.vsh.yaml | 3 ++-
 src/workflows/run_benchmark/main.nf         | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index dc79d85e2..d0a478f06 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -68,7 +68,7 @@ argument_groups:
           A list of transcript assignment methods to run.
         type: string
         multiple: true
-        default: "basic_transcript_assignment:baysor:clustermap:pciseq:comseg:proseg"
+        default: "basic_transcript_assignment:baysor:clustermap:pciseq:comseg:proseg:segger"
       - name: "--count_aggregation_methods"
         description: |
           A list of count aggregation methods to run.
@@ -160,6 +160,7 @@ dependencies:
   - name: methods_transcript_assignment/pciseq
   - name: methods_transcript_assignment/comseg
   - name: methods_transcript_assignment/proseg
+  - name: methods_transcript_assignment/segger
   - name: methods_count_aggregation/basic_count_aggregation
   - name: methods_qc_filter/basic_qc_filter
   - name: methods_calculate_cell_volume/alpha_shapes
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 64cfc7f82..f5cbd27c2 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -138,7 +138,8 @@ workflow run_wf {
     clustermap,
     pciseq,
     comseg,
-    proseg
+    proseg,
+    segger
   ]
   
   segm_ass_ch = segm_ch

From 53e172839fb060747a63f09d304a9f51ed16df08 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 00:43:15 +0200
Subject: [PATCH 04/19] duplicates when FOV stiching cleaned up

---
 src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py b/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
index f9e4fe979..056bebb2d 100644
--- a/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
+++ b/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
@@ -458,6 +458,10 @@ def parse_polygon_z0(row):
     how="left",
     predicate="within",
 )
+# A spot lying inside overlapping cell polygons yields multiple sjoin rows,
+# making `joined` longer than `spots_df`. Keep one match per spot (the first)
+# and realign to the original spot index so the assignment stays 1:1.
+joined = joined[~joined.index.duplicated(keep="first")].reindex(spots_gdf.index)
 spots_df["cell_id"] = joined["cell_id"].values
 print(
     datetime.now() - t0,

From 1186b7a09dbad16d2d4dc96a40fddb23fde9a978 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 09:48:14 +0200
Subject: [PATCH 05/19] chunks issue atera

---
 src/data_processors/process_dataset/script.py | 35 +++++++++++--------
 1 file changed, 20 insertions(+), 15 deletions(-)

diff --git a/src/data_processors/process_dataset/script.py b/src/data_processors/process_dataset/script.py
index 28f7e57b8..68f32cb8e 100644
--- a/src/data_processors/process_dataset/script.py
+++ b/src/data_processors/process_dataset/script.py
@@ -8,6 +8,11 @@
 import tempfile
 from pathlib import Path
 
+try:
+    from xarray import DataTree
+except ImportError:  # older xarray ships DataTree as a separate package
+    from datatree import DataTree
+
 # The 10x Atera loader (newer spatialdata/zarr) writes image arrays with
 # rectilinear (variable-size) chunk grids, which zarr>=3 gates behind an
 # experimental flag and refuses to read by default. Enable reading them so
@@ -85,21 +90,21 @@ def rechunk_sdata(sdata, CHUNK_SIZE=1024):
         
     """
     
-    for key in list(sdata.images.keys()):
-        image = sdata.images[key]
-        coords = list(image["scale0"].coords.keys())
-        rechunk_strategy = {c: CHUNK_SIZE for c in coords}
-        if "c" in coords:
-            rechunk_strategy["c"] = image["scale0"]["image"].chunks[0][0]
-        image = image.chunk(rechunk_strategy)
-        sdata.images[key] = image
-        
-    for key in list(sdata.labels.keys()):
-        label_image = sdata.labels[key]
-        coords = list(label_image.coords.keys())
-        rechunk_strategy = {c: CHUNK_SIZE for c in coords}
-        label_image = label_image.chunk(rechunk_strategy)
-        sdata.labels[key] = label_image
+    # Chunk only the spatial dims to a uniform size. Reading coords from the
+    # element directly fails for multiscale rasters (DataTree), whose root
+    # exposes no y/x coords -> an empty strategy leaves the coarse pyramid
+    # levels with their irregular (rectilinear) chunks, which zarr>=3 cannot
+    # write. Read coords from scale0 for DataTrees; keep the channel axis whole.
+    spatial = ("z", "y", "x")
+    for group in (sdata.images, sdata.labels):
+        for key in list(group.keys()):
+            elem = group[key]
+            ref = elem["scale0"] if isinstance(elem, DataTree) else elem
+            coords = list(ref.coords.keys())
+            rechunk_strategy = {c: CHUNK_SIZE for c in coords if c in spatial}
+            if "c" in coords:  # keep the channel axis as a single chunk
+                rechunk_strategy["c"] = -1
+            group[key] = elem.chunk(rechunk_strategy)
 
 
 def subsample_adata_group_balanced(adata, group_key, n_samples, seed=0):

From 18644d7a62a0667876aea1c075cb5a25920ab01c Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 09:52:22 +0200
Subject: [PATCH 06/19] segger update image

---
 .../segger/config.vsh.yaml                           | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
index bb8288d8c..a3b72b8f4 100644
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -69,6 +69,18 @@ engines:
       - type: docker
         run: |
           pip install --no-cache-dir --extra-index-url=https://pypi.nvidia.com cuspatial-cu12
+      # segger imports torch_geometric (at import time via _patches), lightning
+      # (segger segment -> Trainer) and torch_scatter (models: scatter_max) but
+      # declares NONE of them in its pyproject, so the github install below does
+      # not pull them. torch_geometric/lightning are pure Python. torch_scatter
+      # has no prebuilt wheel for the NGC image's custom torch build, so it must
+      # compile against the image's torch (--no-build-isolation). The build host
+      # has no GPU, so force a CUDA build for the target GPU archs (A100=8.0,
+      # A10/A40=8.6, L4/L40=8.9, H100=9.0).
+      - type: docker
+        run: |
+          pip install --no-cache-dir torch_geometric lightning
+          FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
       # segger trunk installed last so its (unpinned) deps resolve against the
       # RAPIDS/torch stack already present in the image.
       - type: python

From ecb302d31df9d88d6851d94afbcd64eff4df9282 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 15:01:20 +0200
Subject: [PATCH 07/19] claude fix for segger

---
 .../allen_brain_cell_atlas_merfish/script.py  |  7 +++-
 .../segger/config.vsh.yaml                    | 41 +++++++++++--------
 2 files changed, 31 insertions(+), 17 deletions(-)

diff --git a/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py b/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
index 056bebb2d..48d2984f3 100644
--- a/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
+++ b/src/datasets/loaders/allen_brain_cell_atlas_merfish/script.py
@@ -75,7 +75,12 @@
 def download_if_missing(url, path):
     if not Path(path).exists():
         print(datetime.now() - t0, f"Downloading {url}", flush=True)
-        urllib.request.urlretrieve(url, path)
+        # Route required downloads through the retrying helper so a transient
+        # network hiccup (e.g. "Connection reset by peer" during the TLS
+        # handshake) does not abort the whole run. robust_urlretrieve is
+        # defined below but only referenced at call time, which is fine.
+        if not robust_urlretrieve(url, path):
+            raise RuntimeError(f"Failed to download required file after retries: {url}")
 
 
 def check_url(url):
diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
index a3b72b8f4..3c6a5ecaf 100644
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -45,7 +45,12 @@ engines:
   #       (cudf, cuspatial, torch). Develop/test on a CUDA host; `viash test`
   #       will not run on a machine without an NVIDIA GPU.
   - type: docker
-    image: nvcr.io/nvidia/pytorch:25.06-py3
+    # NGC 26.06 ships numpy 2.1 and a torch (2.13) compiled against numpy 2.x, so
+    # the torch<->numpy bridge (used by the `segger segment` subprocess) stays
+    # alive alongside the numpy-2.x task stack (spatialdata>=0.7.3 hard-requires
+    # numpy>=2 via multiscale-spatial-image -> xarray-dataclass). The older 25.06
+    # image shipped numpy-1.x torch, which cannot coexist with spatialdata>=0.7.3.
+    image: nvcr.io/nvidia/pytorch:26.06-py3
     setup:
       - type: apt
         packages: [procps, git, libxcb1]
@@ -61,30 +66,34 @@ engines:
           - shapely
           - "rasterio>=1.3"
           - scikit-image
-      # cuspatial is required by segger.geometry.conversion at import time and is
-      # NOT shipped in the NGC PyTorch base image (only cudf is). No version pin,
-      # so pip resolves it against whatever cudf the image already has installed
-      # (pinning would risk diverging from the image's cudf). The NVIDIA index is
-      # required: cuspatial's libcudf-cu12/libcuspatial-cu12 deps are not on PyPI.
+      # cuspatial is imported by segger.geometry.conversion and is NOT shipped in
+      # the 26.06 image (unlike 25.06). Install from the NVIDIA index (its
+      # libcudf-cu12/libcuspatial-cu12 deps are not on PyPI); it pulls a matching
+      # cudf-cu12. Unpinned -> resolves to the newest cuda-12 build (25.4, the last
+      # cuda-12 cuspatial-cu12 on the index).
       - type: docker
         run: |
           pip install --no-cache-dir --extra-index-url=https://pypi.nvidia.com cuspatial-cu12
       # segger imports torch_geometric (at import time via _patches), lightning
       # (segger segment -> Trainer) and torch_scatter (models: scatter_max) but
       # declares NONE of them in its pyproject, so the github install below does
-      # not pull them. torch_geometric/lightning are pure Python. torch_scatter
-      # has no prebuilt wheel for the NGC image's custom torch build, so it must
-      # compile against the image's torch (--no-build-isolation). The build host
-      # has no GPU, so force a CUDA build for the target GPU archs (A100=8.0,
-      # A10/A40=8.6, L4/L40=8.9, H100=9.0).
+      # not pull them. torch_geometric/lightning are pure Python. torch_scatter has
+      # no prebuilt wheel for the NGC custom torch build, so compile against it
+      # (--no-build-isolation). The build host has no GPU, so force a CUDA build for
+      # the target archs (A100=8.0, A10/A40=8.6, L4/L40=8.9, H100=9.0).
       - type: docker
-        run: |
-          pip install --no-cache-dir torch_geometric lightning
-          FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
-      # segger trunk installed last so its (unpinned) deps resolve against the
-      # RAPIDS/torch stack already present in the image.
+        run:
+          - pip install --no-cache-dir torch_geometric lightning
+          - FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
       - type: python
         github: dpeerlab/segger
+      # cudf 25.4 needs pandas<2.2.4, but segger (via scanpy 1.12 / anndata 0.13)
+      # pulls pandas 3.0, which removed pandas.api.types.is_interval that cudf
+      # imports -> `import cudf` crashes. Pin pandas back into cudf's supported
+      # range as the LAST step so it wins over the deps above. numpy stays 2.x.
+      - type: docker
+        run:
+          - pip install --no-cache-dir "pandas>=2.0,<2.2.4"
   - type: native
 
 runners:

From d400ebec013786e3b02b274fe1b87008c3f469c7 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 18:10:23 +0200
Subject: [PATCH 08/19] atera version fix

---
 src/datasets/loaders/tenx_atera/config.vsh.yaml | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/datasets/loaders/tenx_atera/config.vsh.yaml b/src/datasets/loaders/tenx_atera/config.vsh.yaml
index ed74ffeac..fd96307de 100644
--- a/src/datasets/loaders/tenx_atera/config.vsh.yaml
+++ b/src/datasets/loaders/tenx_atera/config.vsh.yaml
@@ -60,8 +60,15 @@ engines:
     setup:
       - type: python
         pypi:
-          - spatialdata-io==0.6.0
-          - spatialdata==0.7.2
+          # spatialdata 0.7.2 + zarr>=3.2 writes multiscale rasters (the xenium
+          # cell/nucleus label pyramids) as rectilinear chunk grids, which zarr
+          # gates behind an experimental flag and which spatialdata cannot read
+          # back (rechunking to a uniform grid does NOT help -- dask's to_zarr
+          # passes per-chunk sizes so the grid is always rectilinear). spatialdata
+          # 0.8.0 (with ome-zarr 0.18 + dask>=2026.3) emits a regular grid again
+          # and round-trips cleanly. Verified on the cluster against the real data.
+          - spatialdata-io==0.7.1
+          - spatialdata==0.8.0
   - type: native
 
 runners:

From 64d7b4e534ce4d9100c74af25efe094765de8c05 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 18:11:19 +0200
Subject: [PATCH 09/19] wf for the custom rnaseq scripts

---
 .../process_scrnaseq/config.vsh.yaml          | 59 +++++++++++++++
 .../workflows/process_scrnaseq/main.nf        | 72 +++++++++++++++++++
 2 files changed, 131 insertions(+)
 create mode 100644 src/datasets/workflows/process_scrnaseq/config.vsh.yaml
 create mode 100644 src/datasets/workflows/process_scrnaseq/main.nf

diff --git a/src/datasets/workflows/process_scrnaseq/config.vsh.yaml b/src/datasets/workflows/process_scrnaseq/config.vsh.yaml
new file mode 100644
index 000000000..9841d81ab
--- /dev/null
+++ b/src/datasets/workflows/process_scrnaseq/config.vsh.yaml
@@ -0,0 +1,59 @@
+name: process_scrnaseq
+namespace: datasets/workflows
+
+# Generic single-cell RNA-seq processing workflow.
+#
+# Unlike the per-dataset process_<name>_sc workflows, this one is loader-agnostic:
+# it takes an already-loaded / standardized raw SC dataset as --input and runs the
+# shared processing chain (log_cp -> hvg -> pca -> knn -> extract_uns_metadata) to
+# produce a file_common_scrnaseq dataset (with a 'normalized' layer, HVG, PCA and
+# kNN) usable as the input_sc reference for process_datasets.
+#
+# Use it for custom references that have no dedicated loader (e.g. the MPII human
+# lung annotation standardized by scripts/create_resources/sc/standardize_mpii_human_lung_sc.py).
+
+argument_groups:
+  - name: Inputs
+    arguments:
+      - type: file
+        name: --input
+        required: true
+        description: |
+          A raw single-cell RNA-seq dataset (loader-equivalent h5ad): a 'counts'
+          layer, cell_type annotations in .obs, feature_name (gene symbols) in .var
+          and dataset_* metadata in .uns.
+  - name: Outputs
+    arguments:
+      - name: "--output_dataset"
+        __merge__: /src/api/file_common_scrnaseq.yaml
+        direction: output
+        required: true
+        default: "$id/dataset.h5ad"
+      - name: "--output_meta"
+        direction: output
+        type: file
+        description: "Dataset metadata"
+        default: "$id/dataset_meta.yaml"
+
+resources:
+  - type: nextflow_script
+    path: main.nf
+    entrypoint: run_wf
+  - path: /common/nextflow_helpers/helper.nf
+
+dependencies:
+  - name: normalization/log_cp
+    repository: datasets
+  - name: processors/pca
+    repository: datasets
+  - name: processors/hvg
+    repository: datasets
+  - name: processors/knn
+    repository: datasets
+  - name: utils/extract_uns_metadata
+    repository: openproblems
+
+runners:
+  - type: nextflow
+    directives:
+      label: [midcpu, midmem, hightime]
diff --git a/src/datasets/workflows/process_scrnaseq/main.nf b/src/datasets/workflows/process_scrnaseq/main.nf
new file mode 100644
index 000000000..4db197416
--- /dev/null
+++ b/src/datasets/workflows/process_scrnaseq/main.nf
@@ -0,0 +1,72 @@
+include { findArgumentSchema } from "${meta.resources_dir}/helper.nf"
+
+workflow auto {
+  findStates(params, meta.config)
+    | meta.workflow.run(
+      auto: [publish: "state"]
+    )
+}
+
+workflow run_wf {
+  take:
+  input_ch
+
+  main:
+  output_ch = input_ch
+
+    // Normalize (log CP10k). Reads the raw counts from --input.
+    | log_cp.run(
+      key: "log_cp10k",
+      fromState: [
+        "input": "input"
+      ],
+      args: [
+        "normalization_id": "log_cp10k",
+        "n_cp": 10000
+      ],
+      toState: [
+        "output_normalized": "output"
+      ]
+    )
+    | hvg.run(
+      fromState: ["input": "output_normalized"],
+      toState: ["output_hvg": "output"]
+    )
+
+    | pca.run(
+      fromState: ["input": "output_hvg"],
+      toState: ["output_pca": "output" ]
+    )
+
+    | knn.run(
+      fromState: ["input": "output_pca"],
+      toState: ["output_knn": "output"]
+    )
+    // add synonym
+    | map{ id, state ->
+      [id, state + [output_dataset: state.output_knn]]
+    }
+
+    | extract_uns_metadata.run(
+      fromState: { id, state ->
+        def schema = findArgumentSchema(meta.config, "output_dataset")
+        // workaround: convert GString to String
+        schema = iterateMap(schema, { it instanceof GString ? it.toString() : it })
+        def schemaYaml = tempFile("schema.yaml")
+        writeYaml(schema, schemaYaml)
+        [
+          "input": state.output_dataset,
+          "schema": schemaYaml
+        ]
+      },
+      toState: ["output_meta": "output"]
+    )
+
+    | setState([
+      "output_dataset": "output_dataset",
+      "output_meta": "output_meta"
+    ])
+
+  emit:
+  output_ch
+}

From 3edfbf192629d8b5c63038dc94392d49d1664514 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 23:32:55 +0200
Subject: [PATCH 10/19] adjust the loader image name

---
 src/base/labels_nebius.config              | 22 +++++++++++++---------
 src/datasets/loaders/tenx_atera/script.py  | 12 ++++++++++--
 src/datasets/loaders/tenx_xenium/script.py | 12 ++++++++++--
 3 files changed, 33 insertions(+), 13 deletions(-)

diff --git a/src/base/labels_nebius.config b/src/base/labels_nebius.config
index 741e88a40..f9eb03fc1 100644
--- a/src/base/labels_nebius.config
+++ b/src/base/labels_nebius.config
@@ -33,7 +33,11 @@ process {
   // Default disk space
   disk = 50.GB
 
-  // Always pull the latest image digest so nodes never serve a stale cached image
+  // Always pull the latest image digest so nodes never serve a stale cached image.
+  // NOTE: a `pod` directive in a withLabel block REPLACES this one (Nextflow does not
+  // merge `pod` across selectors) — so every label that sets `pod` (e.g. a nodeSelector)
+  // MUST also include [imagePullPolicy: 'Always'], or its processes silently revert to
+  // IfNotPresent and can run a stale cached image. Keep that in sync below.
   pod = [[imagePullPolicy: 'Always']]
 
   // Retry for exit codes that have something to do with memory issues
@@ -55,13 +59,13 @@ process {
     // Nebius 48vcpu-192gb nodes: 188 GB allocatable
     memory = { get_memory( 25.GB * task.attempt ) }
     disk = { 50.GB * task.attempt }
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00p5wjvb7bc55b2js']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00p5wjvb7bc55b2js'], [imagePullPolicy: 'Always']]
   }
   withLabel: midmem {
     // Nebius 64vcpu-256gb nodes: 251 GB allocatable
     memory = { get_memory( 50.GB * task.attempt ) }
     disk = { 100.GB * task.attempt }
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00wvqfyde1nfwezs0']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00wvqfyde1nfwezs0'], [imagePullPolicy: 'Always']]
   }
   // highmem: 200 GB base. No hard nodeSelector — the 200 GiB request itself
   // restricts placement to the two large cpu-d3 node groups (251 GiB × 20 and
@@ -96,7 +100,7 @@ process {
     accelerator = 1
     memory = 28.GB
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
@@ -107,7 +111,7 @@ process {
     accelerator = 1
     memory = 26.GB
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
@@ -118,7 +122,7 @@ process {
     accelerator = 1
     memory = 26.GB
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
@@ -128,7 +132,7 @@ process {
     accelerator = 1
     memory = 26.GB
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
@@ -139,7 +143,7 @@ process {
     accelerator = 1
     memory = { [ 100.GB * task.attempt, 150.GB ].min() }
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00dhcgx1xqjskycvc']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00dhcgx1xqjskycvc'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
@@ -150,7 +154,7 @@ process {
     accelerator = 1
     memory = { [ 120.GB * task.attempt, 180.GB ].min() }
     disk = 200.GB
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00jp7hyqr094tmy35']]
+    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00jp7hyqr094tmy35'], [imagePullPolicy: 'Always']]
     containerOptions = { workflow.containerEngine == "singularity" ? '--nv':
        ( workflow.containerEngine == "docker" ? '--gpus all': null ) }
   }
diff --git a/src/datasets/loaders/tenx_atera/script.py b/src/datasets/loaders/tenx_atera/script.py
index c594fb26a..5a1a8c89e 100644
--- a/src/datasets/loaders/tenx_atera/script.py
+++ b/src/datasets/loaders/tenx_atera/script.py
@@ -102,8 +102,16 @@ def extract_zip(input_zip: Path, output_dir: Path, strip_root: bool = False):
     cells_as_circles=False,
 )
 
-# remove morphology_focus
-_ = sdata.images.pop("morphology_focus")
+# Keep exactly one morphology raster named "morphology_mip"; process_dataset
+# renames it to the API-required "image" element. Atera output mirrors Xenium
+# Onboard Analysis v4, which ships only a multi-channel "morphology_focus"
+# (channel 0 is DAPI, which the segmentation methods use via image[0]) and no
+# separate "morphology_mip". Prefer the MIP if present, otherwise fall back to
+# morphology_focus so the dataset always has an image.
+if "morphology_mip" not in sdata.images and "morphology_focus" in sdata.images:
+    sdata["morphology_mip"] = sdata["morphology_focus"]
+if "morphology_focus" in sdata.images:
+    del sdata.images["morphology_focus"]
 
 print("Add uns to table", flush=True)
 new_uns = {
diff --git a/src/datasets/loaders/tenx_xenium/script.py b/src/datasets/loaders/tenx_xenium/script.py
index c17e35cf5..8422da04b 100644
--- a/src/datasets/loaders/tenx_xenium/script.py
+++ b/src/datasets/loaders/tenx_xenium/script.py
@@ -59,8 +59,16 @@
         cells_as_circles=False,
     )
 
-    # remove morphology_focus
-    _ = sdata.images.pop("morphology_focus")
+    # Keep exactly one morphology raster named "morphology_mip"; process_dataset
+    # renames it to the API-required "image" element. Older Xenium ship both a
+    # single-channel "morphology_mip" and a multi-channel "morphology_focus";
+    # newer Xenium (Onboard Analysis v2+) ship only "morphology_focus" (channel 0
+    # is DAPI, which the segmentation methods use via image[0]). Prefer the MIP,
+    # otherwise fall back to morphology_focus so the dataset always has an image.
+    if "morphology_mip" not in sdata.images and "morphology_focus" in sdata.images:
+        sdata["morphology_mip"] = sdata["morphology_focus"]
+    if "morphology_focus" in sdata.images:
+        del sdata.images["morphology_focus"]
 
     print("Add uns to table", flush=True)
     new_uns = {

From cbd2f12216921d46a3bc3b41843a5891753af79f Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Mon, 20 Jul 2026 23:55:10 +0200
Subject: [PATCH 11/19] adjust the memory

---
 src/base/labels_nebius.config              | 14 ++++++++++----
 src/datasets/loaders/tenx_atera/script.py  | 13 +++++++++++++
 src/datasets/loaders/tenx_xenium/script.py | 13 +++++++++++++
 3 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/src/base/labels_nebius.config b/src/base/labels_nebius.config
index f9eb03fc1..bfe4bf2a3 100644
--- a/src/base/labels_nebius.config
+++ b/src/base/labels_nebius.config
@@ -58,13 +58,17 @@ process {
   withLabel: lowmem {
     // Nebius 48vcpu-192gb nodes: 188 GB allocatable
     memory = { get_memory( 25.GB * task.attempt ) }
-    disk = { 50.GB * task.attempt }
+    // Scale disk with retries but cap at 200 GB: disk has no maxMemory-style
+    // clamp, so an uncapped scaled retry can exceed node local-disk capacity
+    // and leave the job Pending until it times out.
+    disk = { [ 50.GB * task.attempt, 200.GB ].min() }
     pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00p5wjvb7bc55b2js'], [imagePullPolicy: 'Always']]
   }
   withLabel: midmem {
     // Nebius 64vcpu-256gb nodes: 251 GB allocatable
     memory = { get_memory( 50.GB * task.attempt ) }
-    disk = { 100.GB * task.attempt }
+    // Scale disk with retries but cap at 200 GB: see lowmem note above.
+    disk = { [ 100.GB * task.attempt, 200.GB ].min() }
     pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00wvqfyde1nfwezs0'], [imagePullPolicy: 'Always']]
   }
   // highmem: 200 GB base. No hard nodeSelector — the 200 GiB request itself
@@ -75,7 +79,8 @@ process {
   // units are binary: 200.GB = 200 GiB, which excludes the 195/188 GiB groups.)
   withLabel: highmem {
     memory = { get_memory( 200.GB * task.attempt ) }
-    disk = { 200.GB * task.attempt }
+    // Scale disk with retries but cap at 200 GB: see lowmem note above.
+    disk =  200.GB 
   }
   // veryhighmem: 400 GB base — genuinely larger than highmem (previously both
   // were 200 GB, so veryhighmem bought no extra RAM). A 400 GiB request only
@@ -83,7 +88,8 @@ process {
   // largest jobs without a hard nodeSelector.
   withLabel: veryhighmem {
     memory = { get_memory( 400.GB * task.attempt ) }
-    disk = { 400.GB * task.attempt }
+    // Scale disk with retries but cap at 200 GB: see lowmem note above.
+    disk = 200.GB 
   }
   withLabel: lowsharedmem {
     containerOptions = { workflow.containerEngine != 'singularity' ? "--shm-size ${String.format("%.0f",task.memory.mega * 0.05)}" : ""}
diff --git a/src/datasets/loaders/tenx_atera/script.py b/src/datasets/loaders/tenx_atera/script.py
index 5a1a8c89e..6a2009ddf 100644
--- a/src/datasets/loaders/tenx_atera/script.py
+++ b/src/datasets/loaders/tenx_atera/script.py
@@ -113,6 +113,19 @@ def extract_zip(input_zip: Path, output_dir: Path, strip_root: bool = False):
 if "morphology_focus" in sdata.images:
     del sdata.images["morphology_focus"]
 
+# Log the morphology image channel names so it's visible in the run log which
+# channel is which — segmentation uses channel 0 (expected to be DAPI).
+try:
+    _mip = sdata["morphology_mip"]
+    try:
+        _arr = _mip["scale0"].image  # multiscale raster
+    except (TypeError, KeyError):
+        _arr = _mip  # single-scale DataArray
+    _channels = [str(c) for c in _arr.coords["c"].values] if "c" in _arr.coords else "<no channel axis>"
+    print(f"morphology image channels (channel 0 used for segmentation): {_channels}", flush=True)
+except Exception as e:
+    print(f"Could not read morphology image channel names: {e}", flush=True)
+
 print("Add uns to table", flush=True)
 new_uns = {
     "dataset_id": par["dataset_id"],
diff --git a/src/datasets/loaders/tenx_xenium/script.py b/src/datasets/loaders/tenx_xenium/script.py
index 8422da04b..b1a2d4b7e 100644
--- a/src/datasets/loaders/tenx_xenium/script.py
+++ b/src/datasets/loaders/tenx_xenium/script.py
@@ -70,6 +70,19 @@
     if "morphology_focus" in sdata.images:
         del sdata.images["morphology_focus"]
 
+    # Log the morphology image channel names so it's visible in the run log which
+    # channel is which — segmentation uses channel 0 (expected to be DAPI).
+    try:
+        _mip = sdata["morphology_mip"]
+        try:
+            _arr = _mip["scale0"].image  # multiscale raster
+        except (TypeError, KeyError):
+            _arr = _mip  # single-scale DataArray
+        _channels = [str(c) for c in _arr.coords["c"].values] if "c" in _arr.coords else "<no channel axis>"
+        print(f"morphology image channels (channel 0 used for segmentation): {_channels}", flush=True)
+    except Exception as e:
+        print(f"Could not read morphology image channel names: {e}", flush=True)
+
     print("Add uns to table", flush=True)
     new_uns = {
         "dataset_id": par["dataset_id"],

From 184260e3cbe68629b89b4ab9439a6c44b1cf5eff Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 01:28:53 +0200
Subject: [PATCH 12/19] troubleshootig edges

---
 .../basic_transcript_assignment/script.py        |  5 +++++
 .../baysor/script.py                             |  6 ++++++
 .../proseg/script.py                             |  5 +++++
 .../segger/config.vsh.yaml                       | 16 +++++++++++-----
 4 files changed, 27 insertions(+), 5 deletions(-)

diff --git a/src/methods_transcript_assignment/basic_transcript_assignment/script.py b/src/methods_transcript_assignment/basic_transcript_assignment/script.py
index 20377b8aa..33ff24337 100644
--- a/src/methods_transcript_assignment/basic_transcript_assignment/script.py
+++ b/src/methods_transcript_assignment/basic_transcript_assignment/script.py
@@ -59,6 +59,11 @@
 # Assign cell ids directly to transcripts_reset (clean-index single-partition dask DataFrame).
 # Using sdata[par['transcripts_key']] here would reintroduce the duplicate parquet index,
 # causing the same "cannot reindex on an axis with duplicate labels" error at write time.
+# Clamp coords to the label-image bounds: transcripts at the crop boundary can
+# round a few pixels past the raster edge (see get_crop_coords in
+# process_dataset). Matches segger's handling; edge/background at the border.
+y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
+x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
 cell_ids = label_image[y_coords, x_coords]
 transcripts_reset["cell_id"] = pd.Series(cell_ids)
 
diff --git a/src/methods_transcript_assignment/baysor/script.py b/src/methods_transcript_assignment/baysor/script.py
index 2e6ca0928..748644adb 100644
--- a/src/methods_transcript_assignment/baysor/script.py
+++ b/src/methods_transcript_assignment/baysor/script.py
@@ -80,6 +80,12 @@
 else:
     label_image = sdata_segm["segmentation"].to_numpy()
 
+# Clamp coords to the label-image bounds: transcripts at the crop boundary can
+# round a few pixels past the raster edge (see get_crop_coords in
+# process_dataset). Matches segger's handling; edge/background at the border.
+y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
+x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
+
 cell_id_dask_series = dask.dataframe.from_dask_array(
     dask.array.from_array(
         label_image[y_coords, x_coords], chunks=tuple(sdata[par['transcripts_key']].map_partitions(len).compute())
diff --git a/src/methods_transcript_assignment/proseg/script.py b/src/methods_transcript_assignment/proseg/script.py
index 2c1acd0f2..8de0afe7d 100644
--- a/src/methods_transcript_assignment/proseg/script.py
+++ b/src/methods_transcript_assignment/proseg/script.py
@@ -75,6 +75,11 @@
 # assign transcripts to cells based on x,y coords and segmentation image
 y_coords = transcripts.y.compute().to_numpy(dtype=np.int64)
 x_coords = transcripts.x.compute().to_numpy(dtype=np.int64)
+# Clamp coords to the label-image bounds: transcripts at the crop boundary can
+# round a few pixels past the raster edge (see get_crop_coords in
+# process_dataset). Matches segger's handling; edge/background at the border.
+y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
+x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
 cell_id_dask_series = dask.dataframe.from_dask_array(
     dask.array.from_array(
         label_image[y_coords, x_coords], chunks=tuple(sdata[par['transcripts_key']].map_partitions(len).compute())
diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
index 3c6a5ecaf..532ecb443 100644
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -87,13 +87,19 @@ engines:
           - FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
       - type: python
         github: dpeerlab/segger
-      # cudf 25.4 needs pandas<2.2.4, but segger (via scanpy 1.12 / anndata 0.13)
-      # pulls pandas 3.0, which removed pandas.api.types.is_interval that cudf
-      # imports -> `import cudf` crashes. Pin pandas back into cudf's supported
-      # range as the LAST step so it wins over the deps above. numpy stays 2.x.
+      # cudf 25.4 (cuda-12; there is no cuspatial-cu13, so we are stuck on this
+      # RAPIDS) hard-requires pandas<2.2.4. But segger pulls the modern single-cell
+      # stack — anndata 0.13 and scanpy 1.12 — which BOTH require pandas>=2.3
+      # (anndata 0.13's zarr reader calls pd.StringDtype(na_value=...), an API added
+      # in pandas 2.3), so under the pinned pandas 2.2.3 reading ANY AnnData/zarr
+      # table crashes with "StringDtype.__init__() got an unexpected keyword
+      # argument 'na_value'". Pin pandas back into cudf's range AND hold anndata /
+      # scanpy at their last pandas-2.2-compatible majors (anndata 0.12.x needs
+      # pandas>=2.1; scanpy 1.11.x needs pandas>=1.5) as the LAST step so they win
+      # over the deps above. numpy stays 2.x. See NOTES.md for the full conflict.
       - type: docker
         run:
-          - pip install --no-cache-dir "pandas>=2.0,<2.2.4"
+          - pip install --no-cache-dir "pandas>=2.0,<2.2.4" "anndata>=0.12,<0.13" "scanpy>=1.11,<1.12"
   - type: native
 
 runners:

From 36631c4f61cc71f4da6732be8038197a9d773672 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 10:56:47 +0200
Subject: [PATCH 13/19] segger update

---
 .../segger/config.vsh.yaml                    | 20 +++++++++++++------
 .../segger/script.py                          | 18 ++++++++++-------
 2 files changed, 25 insertions(+), 13 deletions(-)

diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
index 532ecb443..8dd320557 100644
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -24,17 +24,17 @@ arguments:
     required: false
     description: "Number of training epochs for the segger GNN."
     default: 20
-  - name: --prediction_expansion_ratio
+  - name: --prediction_graph_buffer_ratio
     type: double
     required: false
-    description: "Fraction of each polygon's equivalent radius used to expand it during prediction."
-    default: 0.5
+    description: "Fraction of each polygon's equivalent radius used to buffer (expand) it when building the prediction graph. Maps to segger's --prediction-graph-buffer-ratio (renamed from --prediction-expansion-ratio in segger 0.1.0)."
+    default: 0.05
   - name: --prediction_mode
     type: string
     required: false
     choices: [nucleus, cell, uniform]
     description: "Which polygon set drives the prediction graph."
-    default: nucleus
+    default: cell
 
 resources:
   - type: python_script
@@ -52,8 +52,11 @@ engines:
     # image shipped numpy-1.x torch, which cannot coexist with spatialdata>=0.7.3.
     image: nvcr.io/nvidia/pytorch:26.06-py3
     setup:
+      # libgl1 + libglib2.0-0: segger's writer imports segger.io -> cv2 (opencv),
+      # which needs libGL.so.1 / libgthread at import; without them the run crashes
+      # at write time with "libGL.so.1: cannot open shared object file".
       - type: apt
-        packages: [procps, git, libxcb1]
+        packages: [procps, git, libxcb1, libgl1, libglib2.0-0]
       - type: python
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
@@ -85,8 +88,13 @@ engines:
         run:
           - pip install --no-cache-dir torch_geometric lightning
           - FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
+      # Pin to a specific commit: segger's CLI flags AND output schema change on
+      # trunk (0.1.0 renamed --prediction-expansion-ratio -> --prediction-graph-
+      # buffer-ratio and dropped the `keep` column). script.py targets THIS commit's
+      # contract, so bump this deliberately and re-check `segger segment --help`
+      # and data/writer.py before doing so.
       - type: python
-        github: dpeerlab/segger
+        github: dpeerlab/segger@0233cf62eb177b6e44e49318037705550765a010
       # cudf 25.4 (cuda-12; there is no cuspatial-cu13, so we are stuck on this
       # RAPIDS) hard-requires pandas<2.2.4. But segger pulls the modern single-cell
       # stack — anndata 0.13 and scanpy 1.12 — which BOTH require pandas>=2.3
diff --git a/src/methods_transcript_assignment/segger/script.py b/src/methods_transcript_assignment/segger/script.py
index ff43e2662..373aa7deb 100644
--- a/src/methods_transcript_assignment/segger/script.py
+++ b/src/methods_transcript_assignment/segger/script.py
@@ -22,8 +22,8 @@
   'coordinate_system': 'global',
   'output': './temp/methods/segger/segger_assigned_transcripts.zarr',
   'n_epochs': 20,
-  'prediction_expansion_ratio': 0.5,
-  'prediction_mode': 'nucleus',
+  'prediction_graph_buffer_ratio': 0.05,
+  'prediction_mode': 'cell',
 }
 meta = {
   'name': 'segger',
@@ -122,8 +122,10 @@
 else:
     label_image = sdata_segm["segmentation"].to_numpy()
 # Clip to the label image bounds (transcripts can sit just outside after transform).
+n_oob = int(np.count_nonzero((y_coords < 0) | (y_coords >= label_image.shape[0]) | (x_coords < 0) | (x_coords >= label_image.shape[1])))
 y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
 x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
+print(f"Clamped {n_oob}/{len(x_coords)} transcripts outside the {label_image.shape[0]}x{label_image.shape[1]} label image to its edge", flush=True)
 prior_cell_id = label_image[y_coords, x_coords].astype(np.int64)  # 0 == background
 
 # Canonical transcripts frame (native coordinates, clean RangeIndex). Its row
@@ -162,7 +164,7 @@
     "-i", str(XENIUM_DIR),
     "-o", str(SEGGER_OUT_DIR),
     "--n-epochs", str(par["n_epochs"]),
-    "--prediction-expansion-ratio", str(par["prediction_expansion_ratio"]),
+    "--prediction-graph-buffer-ratio", str(par["prediction_graph_buffer_ratio"]),
     "--prediction-mode", par["prediction_mode"],
 ]
 # cudf's `validate_setup` aborts `import cudf` on hosts without a usable NVIDIA
@@ -185,13 +187,15 @@
 print('Reading segger output', flush=True)
 seg = pd.read_parquet(seg_pq)
 
-# Keep only confident assignments. Mirror segger's canonical valid_cell_id
-# predicate (drop null / "-1" / "UNASSIGNED" / "NONE", case-insensitive) so
-# the unassigned sentinels don't leak into the output.
+# Keep only confident assignments. segger 0.1.0 no longer emits a `keep` column;
+# a transcript is confidently assigned when its per-gene similarity clears the
+# learned threshold (segger applies exactly this `segger_similarity >=
+# similarity_threshold` filter when writing its own AnnData). Also drop
+# null / "-1" / "UNASSIGNED" / "NONE" cell ids (case-insensitive).
 _UNASSIGNED = {"-1", "UNASSIGNED", "NONE", "NAN", ""}
 seg_id_str = seg["segger_cell_id"].astype(str).str.strip().str.upper()
 keep_mask = (
-    seg["keep"].astype(bool)
+    (seg["segger_similarity"] >= seg["similarity_threshold"])
     & seg["segger_cell_id"].notna()
     & ~seg_id_str.isin(_UNASSIGNED)
 )

From 06261271a3aa8ed64dd06674497645ef66680858 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 10:57:43 +0200
Subject: [PATCH 14/19] cell type label correction

---
 .../rctd/script.R                             | 21 +++++++++--
 .../split/script.R                            | 35 +++++++++++--------
 2 files changed, 40 insertions(+), 16 deletions(-)

diff --git a/src/methods_cell_type_annotation/rctd/script.R b/src/methods_cell_type_annotation/rctd/script.R
index bcc28cfd2..6e4524e39 100644
--- a/src/methods_cell_type_annotation/rctd/script.R
+++ b/src/methods_cell_type_annotation/rctd/script.R
@@ -44,8 +44,20 @@ filtered_ref <- ref[,colData(ref)$cell_type %in% valid_celltypes]
 ref_counts <- assay(filtered_ref, "counts")
 # factor to drop filtered cell types
 colData(filtered_ref)$cell_type <- factor(colData(filtered_ref)$cell_type) 
-cell_types <- colData(filtered_ref)$cell_type 
+cell_types <- colData(filtered_ref)$cell_type
 names(cell_types) <- colnames(ref_counts)
+
+# spacexr::check_cell_types() rejects cell-type names containing '/'
+# (e.g. "Ciliated/secretory cells", "T/NK lineage"). Sanitize the factor levels
+# before building the Reference, and keep a map (safe name -> original name) so
+# the predicted labels can be restored below to match the scRNAseq reference.
+cell_type_name_map <- levels(cell_types)
+names(cell_type_name_map) <- gsub("/", "_", levels(cell_types))
+if (any(duplicated(names(cell_type_name_map)))) {
+  names(cell_type_name_map) <- make.unique(names(cell_type_name_map))
+}
+levels(cell_types) <- names(cell_type_name_map)
+
 reference <- Reference(ref_counts, cell_types, min_UMI = 0)
 
 # check cores
@@ -75,7 +87,12 @@ names(spatial_cell_types) <- rownames(results$results_df)
 
 #
 colData(sce)$cell_type <- "None_sp"
-colData(sce)[names(spatial_cell_types),"cell_type"] <- as.character(spatial_cell_types)
+# Restore the original cell-type names (undo the '/' sanitization). "None_sp"
+# and anything not in the map are left unchanged.
+predicted <- as.character(spatial_cell_types)
+mapped <- unname(cell_type_name_map[predicted])
+predicted[!is.na(mapped)] <- mapped[!is.na(mapped)]
+colData(sce)[names(spatial_cell_types),"cell_type"] <- predicted
 
 # Write the final object to h5ad format
 # set to 'w', is this ok?
diff --git a/src/methods_expression_correction/split/script.R b/src/methods_expression_correction/split/script.R
index 5210aab16..aa0e06957 100644
--- a/src/methods_expression_correction/split/script.R
+++ b/src/methods_expression_correction/split/script.R
@@ -7,12 +7,12 @@ library(Seurat)
 library(scuttle)
 
 ## VIASH START
-par <- list(
-  "input_spatial_with_cell_types" = "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/spatial_with_celltypes.h5ad",
-  "input_scrnaseq_reference"= "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/scrnaseq_reference.h5ad",
-  "output" = "task_ist_preprocessing/tmp/split_corrected.h5ad",
-  "keep_all_cells" = FALSE,
-)
+par <- list(
+  "input_spatial_with_cell_types" = "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/spatial_with_celltypes.h5ad",
+  "input_scrnaseq_reference"= "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/scrnaseq_reference.h5ad",
+  "output" = "task_ist_preprocessing/tmp/split_corrected.h5ad",
+  "keep_all_cells" = FALSE,
+)
 
 meta <- list(
   'cpus': 4,
@@ -51,6 +51,13 @@ ref_counts <- assay(filtered_ref, "counts")
 colData(filtered_ref)$cell_type <- factor(colData(filtered_ref)$cell_type)
 cell_types <- colData(filtered_ref)$cell_type
 names(cell_types) <- colnames(ref_counts)
+
+# spacexr::check_cell_types() rejects cell-type names containing '/'
+# (e.g. "Ciliated/secretory cells", "T/NK lineage"). Sanitize the factor levels
+# before building the Reference. SPLIT uses RCTD only internally to purify
+# counts and does not emit cell-type labels, so no name restoration is needed.
+levels(cell_types) <- make.unique(gsub("/", "_", levels(cell_types)))
+
 reference <- Reference(ref_counts, cell_types, min_UMI = 0)
 
 # check cores
@@ -83,14 +90,14 @@ res_split <- SPLIT::purify(
 )
 
 
-# create corrected counts layer in original SingleCell object
-cat("Normalizing counts\n")
-
-# Preserve original normalized values before overwriting with corrected normalization
-assay(sce, "normalized_uncorrected") <- assay(sce, "normalized")
-
-# First copy in counts
-assay(sce, "corrected_counts") <- assay(sce, "counts")
+# create corrected counts layer in original SingleCell object
+cat("Normalizing counts\n")
+
+# Preserve original normalized values before overwriting with corrected normalization
+assay(sce, "normalized_uncorrected") <- assay(sce, "normalized")
+
+# First copy in counts
+assay(sce, "corrected_counts") <- assay(sce, "counts")
 
 # Then, replace only the updated cells
 assay(sce, "corrected_counts")[rownames(res_split$purified_counts), colnames(res_split$purified_counts)] <- res_split$purified_counts

From 318643543dd7a0301bf995c1970d83bca5a30340 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 10:58:51 +0200
Subject: [PATCH 15/19] fix boundaries

---
 .../basic_transcript_assignment/script.py     |  2 ++
 .../baysor/script.py                          | 19 +++++++-------
 .../clustermap/script.py                      | 15 +++++------
 .../comseg/script.py                          | 26 ++++++++++++++++++-
 .../proseg/script.py                          | 17 ++++++------
 5 files changed, 53 insertions(+), 26 deletions(-)

diff --git a/src/methods_transcript_assignment/basic_transcript_assignment/script.py b/src/methods_transcript_assignment/basic_transcript_assignment/script.py
index 33ff24337..334e2b529 100644
--- a/src/methods_transcript_assignment/basic_transcript_assignment/script.py
+++ b/src/methods_transcript_assignment/basic_transcript_assignment/script.py
@@ -62,8 +62,10 @@
 # Clamp coords to the label-image bounds: transcripts at the crop boundary can
 # round a few pixels past the raster edge (see get_crop_coords in
 # process_dataset). Matches segger's handling; edge/background at the border.
+n_oob = int(np.count_nonzero((y_coords < 0) | (y_coords >= label_image.shape[0]) | (x_coords < 0) | (x_coords >= label_image.shape[1])))
 y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
 x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
+print(f"Clamped {n_oob}/{len(x_coords)} transcripts outside the {label_image.shape[0]}x{label_image.shape[1]} label image to its edge", flush=True)
 cell_ids = label_image[y_coords, x_coords]
 transcripts_reset["cell_id"] = pd.Series(cell_ids)
 
diff --git a/src/methods_transcript_assignment/baysor/script.py b/src/methods_transcript_assignment/baysor/script.py
index 748644adb..5cca8e4b5 100644
--- a/src/methods_transcript_assignment/baysor/script.py
+++ b/src/methods_transcript_assignment/baysor/script.py
@@ -83,26 +83,27 @@
 # Clamp coords to the label-image bounds: transcripts at the crop boundary can
 # round a few pixels past the raster edge (see get_crop_coords in
 # process_dataset). Matches segger's handling; edge/background at the border.
+n_oob = int(np.count_nonzero((y_coords < 0) | (y_coords >= label_image.shape[0]) | (x_coords < 0) | (x_coords >= label_image.shape[1])))
 y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
 x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
+print(f"Clamped {n_oob}/{len(x_coords)} transcripts outside the {label_image.shape[0]}x{label_image.shape[1]} label image to its edge", flush=True)
 
-cell_id_dask_series = dask.dataframe.from_dask_array(
-    dask.array.from_array(
-        label_image[y_coords, x_coords], chunks=tuple(sdata[par['transcripts_key']].map_partitions(len).compute())
-    ), 
-    index=sdata[par['transcripts_key']].index
-)
-sdata[par['transcripts_key']]["cell_id"] = cell_id_dask_series
+# Assign cell ids to the clean-index single-partition frame (transcripts_reset).
+# Using sdata[par['transcripts_key']] here reintroduces the duplicate parquet index
+# (multi-partition parquet restarts at 0 per partition) and fails downstream with
+# "cannot reindex on an axis with duplicate labels".
+cell_ids = label_image[y_coords, x_coords]
+transcripts_reset["cell_id"] = pd.Series(cell_ids)
 
 
 ########################
 # Run baysor with sopa #
 ########################
 
-# Create reduced sdata
+# Create reduced sdata (use the clean-index frame that carries cell_id)
 sdata_sopa = sd.SpatialData(
     points={
-        "transcripts": sdata[par['transcripts_key']]
+        "transcripts": transcripts_reset
     },
 )
 
diff --git a/src/methods_transcript_assignment/clustermap/script.py b/src/methods_transcript_assignment/clustermap/script.py
index 071b1bf70..87544c1c3 100644
--- a/src/methods_transcript_assignment/clustermap/script.py
+++ b/src/methods_transcript_assignment/clustermap/script.py
@@ -270,13 +270,12 @@ def run_clustermap(
 spots_cmap['clustermap'] += 1 
 
 #assign ClusterMap output (stored in spots_cmap) to the transcripts (i.e. assign transcripts to cells)
-cell_id_dask_series = dask.dataframe.from_dask_array(
-    dask.array.from_array(
-        spots_cmap["clustermap"].values, chunks=tuple(sdata[par['transcripts_key']].map_partitions(len).compute())
-    ), 
-    index=sdata[par['transcripts_key']].index
-)
-sdata[par['transcripts_key']]["cell_id"] = cell_id_dask_series 
+# Assign to the clean-index single-partition frame (transcripts_reset). Using
+# sdata[par['transcripts_key']] here reintroduces the duplicate parquet index
+# (multi-partition parquet restarts at 0 per partition) and fails downstream with
+# "cannot reindex on an axis with duplicate labels".
+cell_ids = spots_cmap["clustermap"].values
+transcripts_reset["cell_id"] = pd.Series(cell_ids)
 
 #create new .obs for cells based on the segmentation output (corresponding with the transcripts 'cell_id')
 unique_cells = np.unique(spots_cmap["clustermap"].values)
@@ -298,7 +297,7 @@ def run_clustermap(
 print('Subsetting to transcripts cell id data', flush=True)
 sdata_transcripts_only = sd.SpatialData(
   points={
-    "transcripts": sdata[par['transcripts_key']]
+    "transcripts": transcripts_reset
   },
   tables={
     "table": ad.AnnData(
diff --git a/src/methods_transcript_assignment/comseg/script.py b/src/methods_transcript_assignment/comseg/script.py
index 00c4330f4..2829b873f 100644
--- a/src/methods_transcript_assignment/comseg/script.py
+++ b/src/methods_transcript_assignment/comseg/script.py
@@ -1,8 +1,10 @@
+import os
 import spatialdata as sd
 import sopa
 import anndata as ad
 import pandas as pd
 import numpy as np
+import dask.dataframe as dd
 
 ## VIASH START
 par = {
@@ -23,6 +25,10 @@
     "norm_vector": False,
     "allow_disconnected_polygon": True,
 }
+meta = {
+    "name": "comseg",
+    "cpus": 15,
+}
 ## VIASH END
 
 
@@ -31,6 +37,16 @@
 sdata = sd.read_zarr(par['input_ist'])
 sdata_segm = sd.read_zarr(par['input_segmentation'])
 
+# Multi-partition parquet restarts its index at 0 per partition, producing duplicate
+# index labels that break sopa's cell_id assignment and the final write with
+# "cannot reindex on an axis with duplicate labels". Rebuild the transcripts as a
+# single-partition frame with a clean RangeIndex before any sopa op touches them.
+_tx = sdata[par['transcripts_key']]
+_tx_reset = dd.from_pandas(_tx.compute().reset_index(drop=True), npartitions=1)
+_tx_reset.attrs.update(_tx.attrs)
+del sdata[par['transcripts_key']]
+sdata[par['transcripts_key']] = _tx_reset
+
 # Convert the prior segmentation to polygons
 sdata["segmentation_boundaries"] = sd.to_polygons(sdata_segm["segmentation"])
 del sdata["segmentation_boundaries"]["label"] # make_transcript_patches will create a new label column and fails if one exists.
@@ -61,7 +77,15 @@
 }
 
 
-# sopa.settings.parallelization_backend = 'dask' #TODO: get parallelization running.
+# 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)
+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)
+
 sopa.segmentation.comseg(sdata, config)
 
 # Assign transcripts to cell ids
diff --git a/src/methods_transcript_assignment/proseg/script.py b/src/methods_transcript_assignment/proseg/script.py
index 8de0afe7d..b11049ce4 100644
--- a/src/methods_transcript_assignment/proseg/script.py
+++ b/src/methods_transcript_assignment/proseg/script.py
@@ -78,15 +78,16 @@
 # Clamp coords to the label-image bounds: transcripts at the crop boundary can
 # round a few pixels past the raster edge (see get_crop_coords in
 # process_dataset). Matches segger's handling; edge/background at the border.
+n_oob = int(np.count_nonzero((y_coords < 0) | (y_coords >= label_image.shape[0]) | (x_coords < 0) | (x_coords >= label_image.shape[1])))
 y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1)
 x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1)
-cell_id_dask_series = dask.dataframe.from_dask_array(
-    dask.array.from_array(
-        label_image[y_coords, x_coords], chunks=tuple(sdata[par['transcripts_key']].map_partitions(len).compute())
-    ), 
-    index=sdata[par['transcripts_key']].index
-)
-sdata[par['transcripts_key']]["cell_id"] = cell_id_dask_series
+print(f"Clamped {n_oob}/{len(x_coords)} transcripts outside the {label_image.shape[0]}x{label_image.shape[1]} label image to its edge", flush=True)
+# Assign cell ids to the clean-index single-partition frame (transcripts_reset).
+# Using sdata[par['transcripts_key']] here reintroduces the duplicate parquet index
+# (multi-partition parquet restarts at 0 per partition) and fails downstream with
+# "cannot reindex on an axis with duplicate labels".
+cell_ids = label_image[y_coords, x_coords]
+transcripts_reset["cell_id"] = pd.Series(cell_ids)
 
 
 ### Run Proseg with sopa
@@ -95,7 +96,7 @@
 print("Creating sopa SpatialData object")
 sdata_sopa = sd.SpatialData(
     points={
-        "transcripts": sdata[par['transcripts_key']]
+        "transcripts": transcripts_reset
     },
 )
 

From 3505718f559c9ad52659460653190d1b47d799c7 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 11:42:02 +0200
Subject: [PATCH 16/19] OOM fixes

---
 .../moscot/config.vsh.yaml                           |  9 +++++++++
 src/methods_cell_type_annotation/moscot/script.py    |  6 ++++++
 .../tangram/config.vsh.yaml                          |  5 ++++-
 src/methods_cell_type_annotation/tangram/script.py   | 12 ++++++++++--
 4 files changed, 29 insertions(+), 3 deletions(-)

diff --git a/src/methods_cell_type_annotation/moscot/config.vsh.yaml b/src/methods_cell_type_annotation/moscot/config.vsh.yaml
index ea79aa1b5..c29da2da5 100644
--- a/src/methods_cell_type_annotation/moscot/config.vsh.yaml
+++ b/src/methods_cell_type_annotation/moscot/config.vsh.yaml
@@ -34,6 +34,15 @@ arguments:
     direction: input
     type: integer
     default: 500 #5000
+  - name: --batch_size
+    required: false
+    direction: input
+    type: integer
+    # Number of cost-matrix rows/cols materialized at once during the solve.
+    # Computes the (n_sc x n_sc) / (n_sp x n_sp) GW costs online in batches instead
+    # of materializing the full matrix, which otherwise exhausts GPU memory on large
+    # references (the 101k-cell MPII lung reference needs ~76 GiB for n_sc x n_sc).
+    default: 1024
   - name: --mapping_mode
     required: false
     direction: input
diff --git a/src/methods_cell_type_annotation/moscot/script.py b/src/methods_cell_type_annotation/moscot/script.py
index 1f5adc847..d69e21e45 100644
--- a/src/methods_cell_type_annotation/moscot/script.py
+++ b/src/methods_cell_type_annotation/moscot/script.py
@@ -41,6 +41,7 @@
    'epsilon': 0.01,
    'tau': 0.3,
    'rank': 500, #5000
+   'batch_size': 1024,
    'mapping_mode': 'max',
 }
 meta = {
@@ -82,12 +83,17 @@
     xy_callback="local-pca",
 )
 
+# batch_size computes the GW cost matrices online in batches instead of
+# materializing the full (n_sc x n_sc) matrix, which otherwise exhausts GPU
+# memory on large references (e.g. the 101k-cell MPII lung reference needs
+# ~76 GiB just for n_sc x n_sc). Trades a little speed for a bounded footprint.
 mp = mp.solve(
     alpha=par['alpha'],
     epsilon=par['epsilon'],
     tau_a=par['tau'],
     tau_b=par['tau'],
     rank=par['rank'],
+    batch_size=par['batch_size'],
 )
 
 # Map annotations
diff --git a/src/methods_cell_type_annotation/tangram/config.vsh.yaml b/src/methods_cell_type_annotation/tangram/config.vsh.yaml
index 903c82c52..f6f9ec48b 100644
--- a/src/methods_cell_type_annotation/tangram/config.vsh.yaml
+++ b/src/methods_cell_type_annotation/tangram/config.vsh.yaml
@@ -15,7 +15,10 @@ arguments:
     required: false
     direction: input
     type: string
-    default: "cells"
+    # 'clusters' maps per-cell-type clusters -> space; 'cells' maps individual
+    # reference cells and exhausts GPU VRAM on large references (e.g. the 101k-cell
+    # MPII lung reference). 'clusters' is the appropriate mode for label projection.
+    default: "clusters"
   - name: --num_epochs
     required: false
     direction: input
diff --git a/src/methods_cell_type_annotation/tangram/script.py b/src/methods_cell_type_annotation/tangram/script.py
index 38f883c55..ffe679ab4 100644
--- a/src/methods_cell_type_annotation/tangram/script.py
+++ b/src/methods_cell_type_annotation/tangram/script.py
@@ -48,14 +48,22 @@
 # 'uniform', when is a ndarray, shape = (number_spots,). 
 # use 'uniform' if the spatial voxels are at single cell resolution (e.g. MERFISH). 'rna_count_based', assumes that 
 # cell density is proportional to the number of RNA molecules.
-adata_map = tg.map_cells_to_space(
+# In 'cells' mode tangram learns a (n_sc_cells x n_spatial_cells) mapping matrix
+# on the GPU, which runs out of VRAM for large references (e.g. the 101k-cell
+# MPII lung reference). 'clusters' mode maps per-cell-type clusters instead,
+# shrinking that matrix by ~n_cells_per_type and fitting comfortably on the GPU;
+# it requires the cluster_label to average cells within.
+map_kwargs = dict(
     adata_sc=adata_sc,
     adata_sp=adata_sp,
     device=device,
     mode=par['mode'],
     num_epochs=par['num_epochs'],
-    density_prior='uniform'
+    density_prior='uniform',
 )
+if par['mode'] == 'clusters':
+    map_kwargs['cluster_label'] = par['celltype_key']
+adata_map = tg.map_cells_to_space(**map_kwargs)
     
 # Spatial prediction dataframe is saved in `obsm` `tangram_ct_pred` of the spatial AnnData
 tg.project_cell_annotations(

From d6e110a71d1430cac56c1167851a5cdb9f12ac6e Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 11:42:15 +0200
Subject: [PATCH 17/19] fix code

---
 .../pciseq/script.py                          | 30 +++++++++++++------
 1 file changed, 21 insertions(+), 9 deletions(-)

diff --git a/src/methods_transcript_assignment/pciseq/script.py b/src/methods_transcript_assignment/pciseq/script.py
index d7aa84d2d..f3001003c 100644
--- a/src/methods_transcript_assignment/pciseq/script.py
+++ b/src/methods_transcript_assignment/pciseq/script.py
@@ -135,15 +135,27 @@ def eta_update_no_assert(self):
 else:
      label_image = sdata_segm["segmentation"].to_numpy()
 
-# There might be spots at the border of the image, pciseq runs into an error if this is the case.
-# Build the mask as a positional numpy array from the transformed coords and apply the SAME mask
-# to both the pciSeq input and the full transcripts, so pciSeq's per-spot assignments line up
-# row-for-row with the transcripts when cell_ids are written back below. (Filtering the original
-# multi-partition dask frame by a mask derived from the reset-index transcripts raises an
-# "Unalignable boolean Series" IndexingError because the two indices no longer match.)
-transcripts_at_border = (x_coords > (label_image.shape[1] - 0.5)) | (y_coords > (label_image.shape[0] - 0.5))
-transcripts_dataframe = transcripts_dataframe[~transcripts_at_border]
-transcripts_full = transcripts_full[~transcripts_at_border]
+# pciSeq internally drops spots that fall OUTSIDE the label image, and returns
+# per-spot assignments only for the kept spots; txsim.run_pciSeq then does
+# `spots['cell'] = assignments`, which raises "Length of values ... does not match
+# length of index" when any spot was dropped. So we must pre-filter to exactly the
+# spots pciSeq keeps. pciSeq rounds coords to pixel indices and keeps round(coord)
+# in [0, size); the equivalent float bounds are (-0.5, size-0.5]. Filter BOTH borders
+# (a crop transform can push transcripts off either the high OR the low/negative side;
+# filtering only the high border left the negative-coord spots in and caused the
+# mismatch on cropped datasets like MPII).
+# Build the mask as a positional numpy array from the transformed coords and apply the
+# SAME mask to both the pciSeq input and the full transcripts so pciSeq's per-spot
+# assignments line up row-for-row when cell_ids are written back below.
+transcripts_out_of_bounds = (
+    (x_coords > (label_image.shape[1] - 0.5)) | (y_coords > (label_image.shape[0] - 0.5)) |
+    (x_coords < -0.5) | (y_coords < -0.5)
+)
+n_oob = int(transcripts_out_of_bounds.sum())
+print(f"Dropping {n_oob}/{len(x_coords)} transcripts outside the "
+      f"{label_image.shape[0]}x{label_image.shape[1]} label image (pciSeq drops these internally)", flush=True)
+transcripts_dataframe = transcripts_dataframe[~transcripts_out_of_bounds]
+transcripts_full = transcripts_full[~transcripts_out_of_bounds]
 
 # Grab all the pciSeq parameters
 opts_keys = [#'exclude_genes',

From 4660f2619a4f8df8dda4127822933bcbd67655d1 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 22:58:52 +0200
Subject: [PATCH 18/19] RCTD

---
 .../rctd/script.R                             | 19 ++++++++++++-
 .../split/config.vsh.yaml                     | 28 +++++++++++++++++++
 .../split/script.R                            | 15 +++++++++-
 3 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/src/methods_cell_type_annotation/rctd/script.R b/src/methods_cell_type_annotation/rctd/script.R
index 6e4524e39..93a9cc345 100644
--- a/src/methods_cell_type_annotation/rctd/script.R
+++ b/src/methods_cell_type_annotation/rctd/script.R
@@ -43,10 +43,27 @@ filtered_ref <- ref[,colData(ref)$cell_type %in% valid_celltypes]
 
 ref_counts <- assay(filtered_ref, "counts")
 # factor to drop filtered cell types
-colData(filtered_ref)$cell_type <- factor(colData(filtered_ref)$cell_type) 
+colData(filtered_ref)$cell_type <- factor(colData(filtered_ref)$cell_type)
 cell_types <- colData(filtered_ref)$cell_type
 names(cell_types) <- colnames(ref_counts)
 
+# --- Diagnostics: confirm what RCTD actually reads. The reference itself has
+#     ~475/477 panel genes with clear cross-cell-type fold-change, so a "fewer
+#     than 10 DE genes" failure means the data reaching RCTD is degenerate:
+#     cell types collapsed to ~1, non-integer/normalized counts, or few genes
+#     shared with the spatial puck. This block surfaces which. ---
+cat("=== RCTD input diagnostics ===\n")
+cat(sprintf("spatial puck: %d genes x %d cells\n", nrow(counts), ncol(counts)))
+cat(sprintf("reference (>=25 cells/type): %d genes x %d cells, %d cell types\n",
+            nrow(ref_counts), ncol(ref_counts), length(unique(as.character(cell_types)))))
+print(utils::head(sort(table(as.character(cell_types)), decreasing = TRUE), 8))
+cat(sprintf("genes shared reference<->spatial: %d\n",
+            length(intersect(rownames(ref_counts), rownames(counts)))))
+.rcv <- if (methods::is(ref_counts, "sparseMatrix")) ref_counts@x else as.numeric(ref_counts)
+if (length(.rcv)) cat(sprintf("reference counts: min=%.4g max=%.4g mean=%.4g all-integer=%s\n",
+            min(.rcv), max(.rcv), mean(.rcv), isTRUE(all.equal(.rcv, round(.rcv)))))
+cat("==============================\n")
+
 # spacexr::check_cell_types() rejects cell-type names containing '/'
 # (e.g. "Ciliated/secretory cells", "T/NK lineage"). Sanitize the factor levels
 # before building the Reference, and keep a map (safe name -> original name) so
diff --git a/src/methods_expression_correction/split/config.vsh.yaml b/src/methods_expression_correction/split/config.vsh.yaml
index d6fe7e14c..9d56d1fba 100644
--- a/src/methods_expression_correction/split/config.vsh.yaml
+++ b/src/methods_expression_correction/split/config.vsh.yaml
@@ -17,6 +17,34 @@ arguments:
     type: boolean
     default: false
     description: Whether to keep cells with 0 counts (may cause errors if set to TRUE)
+  # RCTD's default DE-gene / UMI thresholds are tuned for whole-transcriptome
+  # references; on small iST panels (a few hundred genes, many similar fine-grained
+  # cell types) fold-changes are compressed and RCTD finds "fewer than 10 DE genes".
+  # These relaxed defaults mirror the rctd component so create.RCTD below succeeds.
+  - name: --gene_cutoff
+    type: double
+    default: 0.0
+    description: "Minimum normalized mean expression for a gene to be a platform-effect DE gene (RCTD default 0.000125)."
+  - name: --fc_cutoff
+    type: double
+    default: 0.1
+    description: "Minimum log-fold-change for a gene to be a platform-effect DE gene (RCTD default 0.5)."
+  - name: --gene_cutoff_reg
+    type: double
+    default: 0.0
+    description: "Minimum normalized mean expression for a gene to be a regression DE gene (RCTD default 0.0002)."
+  - name: --fc_cutoff_reg
+    type: double
+    default: 0.1
+    description: "Minimum log-fold-change for a gene to be a regression DE gene (RCTD default 0.75)."
+  - name: --umi_min
+    type: integer
+    default: 20
+    description: "Minimum total UMI per spatial cell to be included (RCTD default 100; iST cells are low-count)."
+  - name: --umi_min_sigma
+    type: integer
+    default: 20
+    description: "Minimum UMI for cells used to fit the platform-effect variance (RCTD default 300)."
 
 resources:
   - type: r_script
diff --git a/src/methods_expression_correction/split/script.R b/src/methods_expression_correction/split/script.R
index aa0e06957..7bcc06f00 100644
--- a/src/methods_expression_correction/split/script.R
+++ b/src/methods_expression_correction/split/script.R
@@ -12,6 +12,12 @@ par <- list(
   "input_scrnaseq_reference"= "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/scrnaseq_reference.h5ad",
   "output" = "task_ist_preprocessing/tmp/split_corrected.h5ad",
   "keep_all_cells" = FALSE,
+  "gene_cutoff" = 0.0,
+  "fc_cutoff" = 0.1,
+  "gene_cutoff_reg" = 0.0,
+  "fc_cutoff_reg" = 0.1,
+  "umi_min" = 20,
+  "umi_min_sigma" = 20
 )
 
 meta <- list(
@@ -67,7 +73,14 @@ cat(sprintf("Number of cores: %s\n", cores))
 
 # Run the algorithm
 cat("Running RCTD\n")
-myRCTD <- create.RCTD(puck, reference, max_cores = cores)
+# Relaxed DE-gene / UMI thresholds (see config): RCTD's defaults find "fewer than
+# 10 DE genes" on small iST panels with many fine-grained cell types.
+myRCTD <- create.RCTD(
+  puck, reference, max_cores = cores,
+  gene_cutoff = par$gene_cutoff, fc_cutoff = par$fc_cutoff,
+  gene_cutoff_reg = par$gene_cutoff_reg, fc_cutoff_reg = par$fc_cutoff_reg,
+  UMI_min = par$umi_min, UMI_min_sigma = par$umi_min_sigma
+)
 myRCTD <- run.RCTD(myRCTD, doublet_mode = "doublet")
 
 # Get the "spot_class" annotation from RCTD

From 5abd651be071216aec6a39669e8c550ef139197a Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Tue, 21 Jul 2026 22:59:19 +0200
Subject: [PATCH 19/19] segger to RAPIDS

---
 .../segger/config.vsh.yaml                    | 97 +++++++++----------
 1 file changed, 45 insertions(+), 52 deletions(-)

diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
index 8dd320557..fdd8872c1 100644
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -45,69 +45,62 @@ engines:
   #       (cudf, cuspatial, torch). Develop/test on a CUDA host; `viash test`
   #       will not run on a machine without an NVIDIA GPU.
   - type: docker
-    # NGC 26.06 ships numpy 2.1 and a torch (2.13) compiled against numpy 2.x, so
-    # the torch<->numpy bridge (used by the `segger segment` subprocess) stays
-    # alive alongside the numpy-2.x task stack (spatialdata>=0.7.3 hard-requires
-    # numpy>=2 via multiscale-spatial-image -> xarray-dataclass). The older 25.06
-    # image shipped numpy-1.x torch, which cannot coexist with spatialdata>=0.7.3.
-    image: nvcr.io/nvidia/pytorch:26.06-py3
+    # segger is fundamentally a RAPIDS application (cudf + cuspatial + cugraph +
+    # cuml). Bolting that full stack onto the NGC PyTorch base made RAPIDS's UCX
+    # (cugraph -> raft_dask) collide with the image's HPC-X UCX and abort with heap
+    # corruption (SIGABRT, exit 134) that no pip pin could fix. So we base on the
+    # official RAPIDS image (ships cudf/cugraph/cuml/cuspatial 25.04 + a consistent,
+    # working UCX/dask stack) and layer torch (GNN) + the spatialdata/anndata I/O
+    # stack + segger on top. numpy 2.0.2 / pandas 2.2.3 come from the base.
+    #
+    # NOTE: this image runs as the NON-root 'rapids' user, so `type: apt` will NOT
+    # work in a viash build (needs root). Everything below installs via conda/pip
+    # only (both write to the rapids-owned conda env) — do NOT add apt steps.
+    image: rapidsai/base:25.04-cuda12.8-py3.12
     setup:
-      # libgl1 + libglib2.0-0: segger's writer imports segger.io -> cv2 (opencv),
-      # which needs libGL.so.1 / libgthread at import; without them the run crashes
-      # at write time with "libGL.so.1: cannot open shared object file".
-      - type: apt
-        packages: [procps, git, libxcb1, libgl1, libglib2.0-0]
+      # git (for the segger + openproblems github installs) via conda: apt is
+      # unavailable as non-root and the RAPIDS base ships no git.
+      - type: docker
+        run: conda install -y -c conda-forge git
       - type: python
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
-      - type: python
-        packages:
-          - "spatialdata>=0.7.3"
-          - "anndata>=0.12.0"
-          - "zarr>=3.0.0"
-          - geopandas
-          - shapely
-          - "rasterio>=1.3"
-          - scikit-image
-      # cuspatial is imported by segger.geometry.conversion and is NOT shipped in
-      # the 26.06 image (unlike 25.06). Install from the NVIDIA index (its
-      # libcudf-cu12/libcuspatial-cu12 deps are not on PyPI); it pulls a matching
-      # cudf-cu12. Unpinned -> resolves to the newest cuda-12 build (25.4, the last
-      # cuda-12 cuspatial-cu12 on the index).
-      - type: docker
-        run: |
-          pip install --no-cache-dir --extra-index-url=https://pypi.nvidia.com cuspatial-cu12
-      # segger imports torch_geometric (at import time via _patches), lightning
-      # (segger segment -> Trainer) and torch_scatter (models: scatter_max) but
-      # declares NONE of them in its pyproject, so the github install below does
-      # not pull them. torch_geometric/lightning are pure Python. torch_scatter has
-      # no prebuilt wheel for the NGC custom torch build, so compile against it
-      # (--no-build-isolation). The build host has no GPU, so force a CUDA build for
-      # the target archs (A100=8.0, A10/A40=8.6, L4/L40=8.9, H100=9.0).
+      # torch for the GNN. Standard PyPI cu124 wheel (numpy-2 compatible) coexists
+      # with the RAPIDS cuda-12 libs. Unlike the NGC custom torch, torch_scatter has
+      # a PREBUILT wheel for this torch (the pyg index) — no CUDA compile needed.
       - type: docker
         run:
+          - pip install --no-cache-dir torch==2.6.0 --index-url https://download.pytorch.org/whl/cu124
           - pip install --no-cache-dir torch_geometric lightning
-          - FORCE_CUDA=1 TORCH_CUDA_ARCH_LIST="8.0;8.6;8.9;9.0" pip install --no-cache-dir --no-build-isolation torch_scatter
-      # Pin to a specific commit: segger's CLI flags AND output schema change on
-      # trunk (0.1.0 renamed --prediction-expansion-ratio -> --prediction-graph-
-      # buffer-ratio and dropped the `keep` column). script.py targets THIS commit's
-      # contract, so bump this deliberately and re-check `segger segment --help`
-      # and data/writer.py before doing so.
+          - pip install --no-cache-dir torch_scatter -f https://data.pyg.org/whl/torch-2.6.0+cu124.html
+      # I/O stack. spatialdata pinned <0.8 (resolves to 0.7.1): 0.7.2+ require
+      # dask>=2026.x, but RAPIDS 25.04 hard-pins dask==2025.2.0. 0.7.1 accepts dask
+      # 2025.2.0; pin dask/distributed so spatialdata can't drag them up (which would
+      # break RAPIDS's dask-cudf/cugraph). We validated sd.transform works on 2025.2.0.
+      - type: docker
+        run:
+          - pip install --no-cache-dir "spatialdata>=0.7,<0.8" "zarr>=3.0.0" geopandas shapely "rasterio>=1.3" scikit-image "dask==2025.2.0" "distributed==2025.2.0"
+      # Pin segger to a commit: its CLI flags AND output schema change on trunk
+      # (0.1.0 renamed --prediction-expansion-ratio -> --prediction-graph-buffer-
+      # ratio and dropped the `keep` column). script.py targets THIS commit; bump
+      # deliberately and re-check `segger segment --help` and data/writer.py first.
       - type: python
         github: dpeerlab/segger@0233cf62eb177b6e44e49318037705550765a010
-      # cudf 25.4 (cuda-12; there is no cuspatial-cu13, so we are stuck on this
-      # RAPIDS) hard-requires pandas<2.2.4. But segger pulls the modern single-cell
-      # stack — anndata 0.13 and scanpy 1.12 — which BOTH require pandas>=2.3
-      # (anndata 0.13's zarr reader calls pd.StringDtype(na_value=...), an API added
-      # in pandas 2.3), so under the pinned pandas 2.2.3 reading ANY AnnData/zarr
-      # table crashes with "StringDtype.__init__() got an unexpected keyword
-      # argument 'na_value'". Pin pandas back into cudf's range AND hold anndata /
-      # scanpy at their last pandas-2.2-compatible majors (anndata 0.12.x needs
-      # pandas>=2.1; scanpy 1.11.x needs pandas>=1.5) as the LAST step so they win
-      # over the deps above. numpy stays 2.x. See NOTES.md for the full conflict.
+      # Final step, must be LAST (after segger's github install pulls its deps):
+      #  - cudf 25.4 needs pandas<2.2.4; but anndata 0.13 / scanpy 1.12 (pulled by
+      #    segger) need pandas>=2.3 (anndata 0.13's zarr reader uses
+      #    pd.StringDtype(na_value=...), a pandas-2.3 API) -> reading any zarr table
+      #    crashes under pandas 2.2.3. Hold anndata/scanpy at their last pandas-2.2
+      #    majors (anndata 0.12.x: pandas>=2.1; scanpy 1.11.x: pandas>=1.5), and
+      #    re-assert spatialdata<0.8 + dask==2025.2.0 in case a dep nudged them.
+      #  - swap opencv-python -> opencv-python-headless: segger pulls opencv-python,
+      #    whose cv2 needs libGL.so.1 (a system lib we cannot apt-install as non-root);
+      #    the headless build provides the same cv2 without it. See NOTES.md.
       - type: docker
         run:
-          - pip install --no-cache-dir "pandas>=2.0,<2.2.4" "anndata>=0.12,<0.13" "scanpy>=1.11,<1.12"
+          - pip install --no-cache-dir "pandas>=2.0,<2.2.4" "anndata>=0.12,<0.13" "scanpy>=1.11,<1.12" "spatialdata>=0.7,<0.8" "dask==2025.2.0" "distributed==2025.2.0"
+          - pip uninstall -y opencv-python || true
+          - pip install --no-cache-dir opencv-python-headless
   - type: native
 
 runners: