From 557379eaf4dcac89b8c973f0aee70bf369a9788a Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Fri, 17 Jul 2026 23:04:23 +0200
Subject: [PATCH 1/5] fix segmentation errors

---
 .../binning/config.vsh.yaml                   |  2 +-
 src/methods_segmentation/binning/script.py    | 20 ++++++++++++--
 .../cellpose/config.vsh.yaml                  |  7 +++++
 .../custom_segmentation/script.py             | 20 +++++++++++---
 .../stardist/config.vsh.yaml                  |  2 +-
 src/methods_segmentation/stardist/script.py   | 26 +++++++++++++++++--
 .../watershed/config.vsh.yaml                 |  2 +-
 src/methods_segmentation/watershed/script.py  | 20 ++++++++++++--
 .../baysor/script.py                          | 11 +++++++-
 .../clustermap/script.py                      | 11 +++++++-
 .../pciseq/script.py                          | 11 +++++++-
 .../proseg/script.py                          | 11 +++++++-
 12 files changed, 126 insertions(+), 17 deletions(-)

diff --git a/src/methods_segmentation/binning/config.vsh.yaml b/src/methods_segmentation/binning/config.vsh.yaml
index e78d96d40..c9474b430 100644
--- a/src/methods_segmentation/binning/config.vsh.yaml
+++ b/src/methods_segmentation/binning/config.vsh.yaml
@@ -31,4 +31,4 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [ midtime, midcpu, highmem ]
+      label: [ midtime, midcpu, midmem ]
diff --git a/src/methods_segmentation/binning/script.py b/src/methods_segmentation/binning/script.py
index 44df8d73a..49ead7e0b 100644
--- a/src/methods_segmentation/binning/script.py
+++ b/src/methods_segmentation/binning/script.py
@@ -48,9 +48,25 @@ def convert_to_lower_dtype(arr):
 data_array = xr.DataArray(image, name=f'segmentation', dims=('y', 'x'))
 parsed_data = Labels2DModel.parse(data_array, transformations=transformation)
 sd_output.labels['segmentation'] = parsed_data
+metadata = sdata.tables["metadata"]
+# cell_id is required downstream. Standard Xenium exports carry an explicit
+# "cell_id" column, but some exports (e.g. the Xenium WTA preview used for the
+# Atera dataset) don't — there the per-cell identifier lives in the table's
+# instance_key column (falling back to the obs index).
+instance_key = metadata.uns.get("spatialdata_attrs", {}).get("instance_key")
+if "cell_id" in metadata.obs.columns:
+    cell_id = metadata.obs["cell_id"].values
+elif instance_key and instance_key in metadata.obs.columns:
+    cell_id = metadata.obs[instance_key].values
+else:
+    cell_id = metadata.obs.index.values
+obs = metadata.obs[[]].copy()
+obs["cell_id"] = cell_id
+if "region" in metadata.obs.columns:
+    obs["region"] = metadata.obs["region"].values
 sd_output.tables['table'] = ad.AnnData(
-      obs=sdata.tables["metadata"].obs[["cell_id", "region"]],
-      var=sdata.tables["metadata"].var[[]]
+      obs=obs,
+      var=metadata.var[[]]
     )
 
 print("Writing output", flush=True)
diff --git a/src/methods_segmentation/cellpose/config.vsh.yaml b/src/methods_segmentation/cellpose/config.vsh.yaml
index 1187c12a9..5d9a69312 100644
--- a/src/methods_segmentation/cellpose/config.vsh.yaml
+++ b/src/methods_segmentation/cellpose/config.vsh.yaml
@@ -84,6 +84,13 @@ engines:
     setup:
       - type: python
         pypi: cellpose<4.0.0
+      # Pre-download the pretrained cyto weights into the image (~/.cellpose)
+      # so runtime tasks never hit the cellpose model server (avoids the
+      # intermittent "HTTP Error 504: Gateway Time-out" seen when several
+      # segmentation tasks fetch the weights concurrently).
+      - type: docker
+        run:
+          - "python -c \"from cellpose import models; models.Cellpose(gpu=False, model_type='cyto')\""
     __merge__: 
       - /src/base/setup_txsim_partial.yaml
       - /src/base/setup_spatialdata_partial.yaml
diff --git a/src/methods_segmentation/custom_segmentation/script.py b/src/methods_segmentation/custom_segmentation/script.py
index d29cb9ab7..b5eb53b92 100644
--- a/src/methods_segmentation/custom_segmentation/script.py
+++ b/src/methods_segmentation/custom_segmentation/script.py
@@ -22,16 +22,28 @@
 
 print(f"Copy segmentation from '{par['labels_key']}'", flush=True)
 metadata = sdata.tables["metadata"]
-# Select only the columns that exist — Xenium provides cell_id and region,
-# Vizgen uses different column names (or an empty obs) so we take what's available.
-obs_cols = [c for c in ["cell_id", "region"] if c in metadata.obs.columns]
+# cell_id is required downstream. Standard Xenium exports carry an explicit
+# "cell_id" column, but some exports (e.g. the Xenium WTA preview used for the
+# Atera dataset) don't — there the per-cell identifier lives in the table's
+# instance_key column (falling back to the obs index).
+instance_key = metadata.uns.get("spatialdata_attrs", {}).get("instance_key")
+if "cell_id" in metadata.obs.columns:
+    cell_id = metadata.obs["cell_id"].values
+elif instance_key and instance_key in metadata.obs.columns:
+    cell_id = metadata.obs[instance_key].values
+else:
+    cell_id = metadata.obs.index.values
+obs = metadata.obs[[]].copy()
+obs["cell_id"] = cell_id
+if "region" in metadata.obs.columns:
+    obs["region"] = metadata.obs["region"].values
 sdata_segmentation_only = sd.SpatialData(
   labels={
     "segmentation": sdata[par["labels_key"]]
   },
   tables={
     "table": ad.AnnData(
-      obs=metadata.obs[obs_cols],
+      obs=obs,
       var=metadata.var[[]]
     )
   }
diff --git a/src/methods_segmentation/stardist/config.vsh.yaml b/src/methods_segmentation/stardist/config.vsh.yaml
index b8a6e7dc6..d361e5d9d 100644
--- a/src/methods_segmentation/stardist/config.vsh.yaml
+++ b/src/methods_segmentation/stardist/config.vsh.yaml
@@ -33,7 +33,7 @@ engines:
         pypi:
           - stardist
           - tensorflow==2.17.0
-          - numpy<2.0.0
+          - "numpy>=1.26.0,<2.0.0"
           - scipy<1.15.0
   - type: native
 
diff --git a/src/methods_segmentation/stardist/script.py b/src/methods_segmentation/stardist/script.py
index 0bb0a1795..1e4c94e41 100644
--- a/src/methods_segmentation/stardist/script.py
+++ b/src/methods_segmentation/stardist/script.py
@@ -2,6 +2,12 @@
 import shutil
 from pathlib import Path
 import numpy as np
+# numpy>=1.24 removed the deprecated `np.long` alias, but stardist/csbdeep
+# still reference it. TensorFlow 2.17 forces numpy>=1.26, so we can't
+# downgrade numpy far enough to get it back — restore the alias instead
+# (`np.long` was always just Python's built-in `int` on py3).
+if not hasattr(np, "long"):
+    np.long = int
 import xarray as xr
 import spatialdata as sd
 import anndata as ad
@@ -77,9 +83,25 @@ def do_after(self):
 parsed_labels = sd.models.Labels2DModel.parse(labels_array, transformations=transformation)
 sd_output.labels['segmentation'] = parsed_labels
 
+metadata = sdata.tables["metadata"]
+# cell_id is required downstream. Standard Xenium exports carry an explicit
+# "cell_id" column, but some exports (e.g. the Xenium WTA preview used for the
+# Atera dataset) don't — there the per-cell identifier lives in the table's
+# instance_key column (falling back to the obs index).
+instance_key = metadata.uns.get("spatialdata_attrs", {}).get("instance_key")
+if "cell_id" in metadata.obs.columns:
+    cell_id = metadata.obs["cell_id"].values
+elif instance_key and instance_key in metadata.obs.columns:
+    cell_id = metadata.obs[instance_key].values
+else:
+    cell_id = metadata.obs.index.values
+obs = metadata.obs[[]].copy()
+obs["cell_id"] = cell_id
+if "region" in metadata.obs.columns:
+    obs["region"] = metadata.obs["region"].values
 sd_output.tables['table'] = ad.AnnData(
-      obs=sdata.tables["metadata"].obs[["cell_id", "region"]],
-      var=sdata.tables["metadata"].var[[]]
+      obs=obs,
+      var=metadata.var[[]]
     )
 
 print("Writing output", flush=True)
diff --git a/src/methods_segmentation/watershed/config.vsh.yaml b/src/methods_segmentation/watershed/config.vsh.yaml
index 07bb35289..cd554a435 100644
--- a/src/methods_segmentation/watershed/config.vsh.yaml
+++ b/src/methods_segmentation/watershed/config.vsh.yaml
@@ -173,5 +173,5 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [ hightime, midcpu, highmem ]
+      label: [ hightime, midcpu, midmem ]
 
diff --git a/src/methods_segmentation/watershed/script.py b/src/methods_segmentation/watershed/script.py
index 9b6aca07a..198e30b3a 100644
--- a/src/methods_segmentation/watershed/script.py
+++ b/src/methods_segmentation/watershed/script.py
@@ -48,9 +48,25 @@ def convert_to_lower_dtype(arr):
 parsed_data = Labels2DModel.parse(data_array, transformations=transformation)
 sd_output.labels['segmentation'] = parsed_data
 
+metadata = sdata.tables["metadata"]
+# cell_id is required downstream. Standard Xenium exports carry an explicit
+# "cell_id" column, but some exports (e.g. the Xenium WTA preview used for the
+# Atera dataset) don't — there the per-cell identifier lives in the table's
+# instance_key column (falling back to the obs index).
+instance_key = metadata.uns.get("spatialdata_attrs", {}).get("instance_key")
+if "cell_id" in metadata.obs.columns:
+    cell_id = metadata.obs["cell_id"].values
+elif instance_key and instance_key in metadata.obs.columns:
+    cell_id = metadata.obs[instance_key].values
+else:
+    cell_id = metadata.obs.index.values
+obs = metadata.obs[[]].copy()
+obs["cell_id"] = cell_id
+if "region" in metadata.obs.columns:
+    obs["region"] = metadata.obs["region"].values
 sd_output.tables['table'] = ad.AnnData(
-      obs=sdata.tables["metadata"].obs[["cell_id", "region"]],
-      var=sdata.tables["metadata"].var[[]]
+      obs=obs,
+      var=metadata.var[[]]
     )
 
 print("Writing output", flush=True)
diff --git a/src/methods_transcript_assignment/baysor/script.py b/src/methods_transcript_assignment/baysor/script.py
index c65215890..2e6ca0928 100644
--- a/src/methods_transcript_assignment/baysor/script.py
+++ b/src/methods_transcript_assignment/baysor/script.py
@@ -57,7 +57,16 @@
 assert par['coordinate_system'] in segmentation_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
 
 print('Transforming transcripts coordinates', flush=True)
-transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
+# Multi-partition parquet files each start with a 0-based index, producing duplicate index
+# values in the combined dask DataFrame. sd.transform() internally creates a pd.Series with
+# index=transformed.index; when that dask index is computed it triggers an assign expression
+# that fails on duplicate/lazy indices. Fix: materialize to pandas and rebuild as a single
+# dask partition with a clean RangeIndex before transforming.
+# The original sdata[transcripts_key] is left unchanged so lines below remain consistent.
+transcripts_input = sdata[par['transcripts_key']]
+transcripts_reset = dask.dataframe.from_pandas(transcripts_input.compute().reset_index(drop=True), npartitions=1)
+transcripts_reset.attrs.update(transcripts_input.attrs)
+transcripts = sd.transform(transcripts_reset, to_coordinate_system=par['coordinate_system'])
 
 # In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
 trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()
diff --git a/src/methods_transcript_assignment/clustermap/script.py b/src/methods_transcript_assignment/clustermap/script.py
index 21afde621..071b1bf70 100644
--- a/src/methods_transcript_assignment/clustermap/script.py
+++ b/src/methods_transcript_assignment/clustermap/script.py
@@ -199,7 +199,16 @@ def run_clustermap(
 assert par['coordinate_system'] in segmentation_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
 
 print('Transforming transcripts coordinates', flush=True)
-transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
+# Multi-partition parquet files each start with a 0-based index, producing duplicate index
+# values in the combined dask DataFrame. sd.transform() internally creates a pd.Series with
+# index=transformed.index; when that dask index is computed it triggers an assign expression
+# that fails on duplicate/lazy indices. Fix: materialize to pandas and rebuild as a single
+# dask partition with a clean RangeIndex before transforming.
+# The original sdata[transcripts_key] is left unchanged so lines below remain consistent.
+transcripts_input = sdata[par['transcripts_key']]
+transcripts_reset = dask.dataframe.from_pandas(transcripts_input.compute().reset_index(drop=True), npartitions=1)
+transcripts_reset.attrs.update(transcripts_input.attrs)
+transcripts = sd.transform(transcripts_reset, to_coordinate_system=par['coordinate_system'])
 
 # In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
 trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()
diff --git a/src/methods_transcript_assignment/pciseq/script.py b/src/methods_transcript_assignment/pciseq/script.py
index f0ece4a66..273810d31 100644
--- a/src/methods_transcript_assignment/pciseq/script.py
+++ b/src/methods_transcript_assignment/pciseq/script.py
@@ -100,7 +100,16 @@ def eta_update_no_assert(self):
 
 # Transform transcript coordinates to the coordinate system
 print('Transforming transcripts coordinates', flush=True)
-transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
+# Multi-partition parquet files each start with a 0-based index, producing duplicate index
+# values in the combined dask DataFrame. sd.transform() internally creates a pd.Series with
+# index=transformed.index; when that dask index is computed it triggers an assign expression
+# that fails on duplicate/lazy indices. Fix: materialize to pandas and rebuild as a single
+# dask partition with a clean RangeIndex before transforming.
+# The original sdata[transcripts_key] is left unchanged so lines below remain consistent.
+transcripts_input = sdata[par['transcripts_key']]
+transcripts_reset = dd.from_pandas(transcripts_input.compute().reset_index(drop=True), npartitions=1)
+transcripts_reset.attrs.update(transcripts_input.attrs)
+transcripts = sd.transform(transcripts_reset, to_coordinate_system=par['coordinate_system'])
 
 # In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
 trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()
diff --git a/src/methods_transcript_assignment/proseg/script.py b/src/methods_transcript_assignment/proseg/script.py
index 419de5135..2c1acd0f2 100644
--- a/src/methods_transcript_assignment/proseg/script.py
+++ b/src/methods_transcript_assignment/proseg/script.py
@@ -51,7 +51,16 @@
 
 # Transform transcript coordinates to the coordinate system
 print('Transforming transcripts coordinates', flush=True)
-transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
+# Multi-partition parquet files each start with a 0-based index, producing duplicate index
+# values in the combined dask DataFrame. sd.transform() internally creates a pd.Series with
+# index=transformed.index; when that dask index is computed it triggers an assign expression
+# that fails on duplicate/lazy indices. Fix: materialize to pandas and rebuild as a single
+# dask partition with a clean RangeIndex before transforming.
+# The original sdata[transcripts_key] is left unchanged so lines below remain consistent.
+transcripts_input = sdata[par['transcripts_key']]
+transcripts_reset = dask.dataframe.from_pandas(transcripts_input.compute().reset_index(drop=True), npartitions=1)
+transcripts_reset.attrs.update(transcripts_input.attrs)
+transcripts = sd.transform(transcripts_reset, to_coordinate_system=par['coordinate_system'])
 
 # In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
 trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()

From 6a361745aa336d9da0425761e61849a8cdf66ab3 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Fri, 17 Jul 2026 23:05:13 +0200
Subject: [PATCH 2/5] adjust memory

---
 src/base/labels_nebius.config | 23 ++++++++++++++++++-----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/src/base/labels_nebius.config b/src/base/labels_nebius.config
index bce93b749..2106809e9 100644
--- a/src/base/labels_nebius.config
+++ b/src/base/labels_nebius.config
@@ -37,7 +37,12 @@ process {
   // always retry once
   errorStrategy = { exitStrat(task) }
   maxRetries = 3
-  maxMemory = null
+  // Cap retry memory so escalated retries stay schedulable. The largest nodes
+  // (highmem node group) have ~503 GiB allocatable, so leave headroom for
+  // daemonsets/system pods. With maxMemory=null the get_memory() clamp is
+  // disabled and highmem's attempt-3 requests 600 GB (200.GB * task.attempt),
+  // which no node can satisfy — the job then sits Pending until it times out.
+  maxMemory = 480.GB
 
   // Resource labels
   withLabel: lowcpu { cpus = 5 }
@@ -51,15 +56,23 @@ process {
     memory = { get_memory( 50.GB * task.attempt ) }
     disk = { 100.GB * task.attempt }
   }
+  // 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
+  // 503 GiB × 10 nodes); no smaller group can hold it. This gives 30 eligible
+  // nodes instead of the 10 a single selector allowed, and lets the 400 GB
+  // retry escalate onto the 503 GiB nodes automatically. (Nextflow memory
+  // 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 }
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00hnqhyfdcsy9m09n']]
   }
+  // veryhighmem: 400 GB base — genuinely larger than highmem (previously both
+  // were 200 GB, so veryhighmem bought no extra RAM). A 400 GiB request only
+  // fits the 503 GiB group (10 nodes), so it self-reserves those nodes for the
+  // largest jobs without a hard nodeSelector.
   withLabel: veryhighmem {
-    memory = { get_memory( 200.GB * task.attempt ) }
+    memory = { get_memory( 400.GB * task.attempt ) }
     disk = { 400.GB * task.attempt }
-    pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00hnqhyfdcsy9m09n']]
   }
   withLabel: lowsharedmem {
     containerOptions = { workflow.containerEngine != 'singularity' ? "--shm-size ${String.format("%.0f",task.memory.mega * 0.05)}" : ""}
@@ -139,7 +152,7 @@ def get_memory(to_compare) {
       return process.maxMemory
     }
     else if (to_compare.compareTo(process.maxMemory as nextflow.util.MemoryUnit) == 1) {
-      return max_memory as nextflow.util.MemoryUnit
+      return process.maxMemory as nextflow.util.MemoryUnit
     }
     else {
       return to_compare

From d6c36adacc75a0bcd7ed9f5e5300fad10a63b8b3 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Sun, 19 Jul 2026 12:15:00 +0200
Subject: [PATCH 3/5] troubleshoot stardist versions, pciseq debug and rctd
 soften thresholds

---
 .../rctd/config.vsh.yaml                      |  31 ++
 .../rctd/script.R                             |  18 +-
 src/methods_segmentation/stardist/script.py   |  17 +-
 .../pciseq/script.py                          |  23 +-
 .../segger/config.vsh.yaml                    | 101 ++++++
 .../segger/script.py                          | 338 ++++++++++++++++++
 6 files changed, 511 insertions(+), 17 deletions(-)
 create mode 100644 src/methods_transcript_assignment/segger/config.vsh.yaml
 create mode 100644 src/methods_transcript_assignment/segger/script.py

diff --git a/src/methods_cell_type_annotation/rctd/config.vsh.yaml b/src/methods_cell_type_annotation/rctd/config.vsh.yaml
index 5c0b94691..bdeccd62b 100644
--- a/src/methods_cell_type_annotation/rctd/config.vsh.yaml
+++ b/src/methods_cell_type_annotation/rctd/config.vsh.yaml
@@ -10,6 +10,37 @@ links:
 references:
   doi: "10.1038/s41587-021-00830-w"
 
+# RCTD's default DE-gene / UMI thresholds are tuned for whole-transcriptome
+# references (~20k genes). Imaging-based ST uses small curated panels (~100-500
+# genes) where fold-changes are compressed and per-cell UMI counts are low, so
+# the defaults yield "fewer than 10 regression differentially expressed genes"
+# and create.RCTD() aborts. These arguments relax the thresholds for iST.
+arguments:
+  - 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
     path: script.R
diff --git a/src/methods_cell_type_annotation/rctd/script.R b/src/methods_cell_type_annotation/rctd/script.R
index ce4fc051f..bcc28cfd2 100644
--- a/src/methods_cell_type_annotation/rctd/script.R
+++ b/src/methods_cell_type_annotation/rctd/script.R
@@ -8,7 +8,13 @@ library(anndataR)
 par <- list(
   "input_spatial_normalized_counts" = "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/spatial_aggregated_counts.h5ad",
   "input_scrnaseq_reference"= "task_ist_preprocessing/resources_test/task_ist_preprocessing/mouse_brain_combined/scrnaseq_reference.h5ad",
-  "output" = "task_ist_preprocessing/tmp/spatial_types.h5ad"
+  "output" = "task_ist_preprocessing/tmp/spatial_types.h5ad",
+  "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(
@@ -48,7 +54,15 @@ if ("cpus" %in% names(meta) && !is.null(meta$cpus)) cores <- meta$cpus
 cat(sprintf("Number of cores: %s\n", cores))
 
 # Run the algorithm
-myRCTD <- create.RCTD(puck, reference, max_cores = cores)
+# NOTE: RCTD's default DE-gene / UMI thresholds are tuned for whole-transcriptome
+# references and produce "fewer than 10 regression differentially expressed genes"
+# on small iST panels. The relaxed thresholds below are passed from the config.
+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")
 
 # Extract results
diff --git a/src/methods_segmentation/stardist/script.py b/src/methods_segmentation/stardist/script.py
index 1e4c94e41..82d72ec4f 100644
--- a/src/methods_segmentation/stardist/script.py
+++ b/src/methods_segmentation/stardist/script.py
@@ -2,12 +2,17 @@
 import shutil
 from pathlib import Path
 import numpy as np
-# numpy>=1.24 removed the deprecated `np.long` alias, but stardist/csbdeep
-# still reference it. TensorFlow 2.17 forces numpy>=1.26, so we can't
-# downgrade numpy far enough to get it back — restore the alias instead
-# (`np.long` was always just Python's built-in `int` on py3).
-if not hasattr(np, "long"):
-    np.long = int
+# numpy>=1.24 removed the deprecated scalar-type aliases (np.bool, np.int,
+# np.float, np.long, ...), but stardist/csbdeep still reference them. TensorFlow
+# 2.17 forces numpy>=1.26, so we can't downgrade numpy far enough to get them
+# back — restore the aliases instead. Each mapped to its Python builtin / numpy
+# type as numpy itself did before removal.
+for _alias, _target in {
+    "bool": bool, "int": int, "float": float, "complex": complex,
+    "object": object, "str": str, "long": int, "unicode": str,
+}.items():
+    if not hasattr(np, _alias):
+        setattr(np, _alias, _target)
 import xarray as xr
 import spatialdata as sd
 import anndata as ad
diff --git a/src/methods_transcript_assignment/pciseq/script.py b/src/methods_transcript_assignment/pciseq/script.py
index 273810d31..d7aa84d2d 100644
--- a/src/methods_transcript_assignment/pciseq/script.py
+++ b/src/methods_transcript_assignment/pciseq/script.py
@@ -122,7 +122,10 @@ def eta_update_no_assert(self):
 
 #Added for pciSeq
 #TODO this will immediately break when the name of the gene isn't feature_name
-transcripts_dataframe = sdata[par['transcripts_key']].compute()[['feature_name']] 
+# Materialize the full transcripts once, in the same row order as the transformed
+# x/y coordinates above, so every downstream filter stays positionally aligned.
+transcripts_full = sdata[par['transcripts_key']].compute()
+transcripts_dataframe = transcripts_full[['feature_name']].copy()
 transcripts_dataframe['x'] = x_coords
 transcripts_dataframe['y'] = y_coords
 
@@ -132,13 +135,15 @@ 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
-transcripts_at_border = transcripts_dataframe['x'] > (label_image.shape[1]-0.5)
-transcripts_at_border = transcripts_at_border | (transcripts_dataframe['y'] > (label_image.shape[0]-0.5))
-transcripts_dataframe = transcripts_dataframe.loc[~transcripts_at_border]
-transcripts_at_border_dask = transcripts.x > (label_image.shape[1]-0.5)
-transcripts_at_border_dask = transcripts_at_border_dask | (transcripts.y > (label_image.shape[0]-0.5))
-sdata[par['transcripts_key']] = sdata[par['transcripts_key']].loc[~transcripts_at_border_dask]
+# 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]
 
 # Grab all the pciSeq parameters
 opts_keys = [#'exclude_genes',
@@ -168,7 +173,7 @@ def eta_update_no_assert(self):
 
 #assign transcript -> cell
 transformations = sd.transformations.get_transformation(sdata[par['transcripts_key']], get_all=True)
-transcripts_pd = sdata[par['transcripts_key']].compute().copy()
+transcripts_pd = transcripts_full.copy()
 transcripts_pd["cell_id"] = assignments['cell'].to_numpy()
 sdata[par['transcripts_key']] = PointsModel.parse(transcripts_pd, transformations=transformations)
 
diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
new file mode 100644
index 000000000..901282c78
--- /dev/null
+++ b/src/methods_transcript_assignment/segger/config.vsh.yaml
@@ -0,0 +1,101 @@
+__merge__: /src/api/comp_method_transcript_assignment.yaml
+
+name: segger
+label: "Segger Transcript Assignment"
+summary: "Assign transcripts to cells using the Segger method"
+description: "Segger is a tool for cell segmentation in single-molecule spatial omics datasets that leverages graph neural networks (GNNs) and heterogeneous graphs."
+links:
+  documentation: "https://elihei2.github.io/segger_dev/"
+  repository: "https://github.com/EliHei2/segger_dev"
+references:
+  doi: "10.1101/2025.03.14.643160"
+
+arguments:
+  - name: --transcripts_key
+    type: string
+    description: The key of the transcripts within the points of the spatial data
+    default: transcripts
+  - name: --coordinate_system
+    type: string
+    description: The key of the pixel space coordinate system within the spatial data
+    default: global
+
+#  - name: --force_2d
+#    type: string
+#    required: false
+#    description: "Ignores z-column in the data if it is provided"
+#    direction: input 
+#    default: "false"
+#
+#  - name: --min_molecules_per_cell
+#    type: integer
+#    required: false
+#    description: "Minimal number of molecules per cell"
+#    direction: input 
+#    default: 50
+#  
+#  - name: --scale
+#    type: double
+#    required: false
+#    description: |
+#      "Scale parameter, which suggest approximate cell radius for the algorithm. Must be in the same units as 
+#      x and y molecule coordinates. Negative values mean it must be estimated from `min_molecules_per_cell`."
+#    direction: input 
+#    default: -1.0
+#
+#  - name: --scale_std
+#    type: string
+#    required: false
+#    description: "Standard deviation of scale across cells relative to `scale`"
+#    direction: input 
+#    default: "25%"
+#
+#  - name: --n_clusters
+#    type: integer
+#    required: false
+#    description: "Number of molecule clusters, i.e. major cell types."
+#    direction: input 
+#    default: 4
+#
+#  - name: --prior_segmentation_confidence
+#    type: double
+#    required: false
+#    description: "Confidence of the prior segmentation"
+#    direction: input 
+#    default: 0.8
+
+resources:
+  - type: python_script
+    path: script.py
+
+engines:
+#  - type: docker
+#    image: openproblems/base_python:1
+#    __merge__: 
+#      - /src/base/setup_spatialdata_partial.yaml
+#    setup:
+#      - type: python
+#        pypi: [torch, pytorch-lightning, lightning]
+#      - type: docker
+#        run:
+#          - git clone https://github.com/EliHei2/segger_dev.git
+#          - cd segger_dev && pip install .[cuda12]
+engines:
+  - type: docker
+    image: danielunyi42/segger_dev:cuda121
+    setup:
+      - type: apt
+        packages: procps
+      - type: python
+        github:
+          - openproblems-bio/core#subdirectory=packages/python/openproblems
+      - type: python
+        packages:
+          - spatialdata
+  - type: native
+
+runners:
+  - type: executable
+  - type: nextflow
+    directives:
+      label: [ midtime, midcpu, midmem, gpu ]
diff --git a/src/methods_transcript_assignment/segger/script.py b/src/methods_transcript_assignment/segger/script.py
new file mode 100644
index 000000000..84cf7e976
--- /dev/null
+++ b/src/methods_transcript_assignment/segger/script.py
@@ -0,0 +1,338 @@
+from pathlib import Path
+import dask
+import torch
+import xarray as xr
+import numpy as np
+import geopandas as gpd
+import spatialdata as sd
+from segger.data.parquet.sample import STSampleParquet
+from segger.training.segger_data_module import SeggerDataModule
+from segger.training.train import LitSegger, Segger
+from torch_geometric.nn import to_hetero
+from lightning.pytorch import Trainer
+#from pytorch_lightning import Trainer
+from lightning.pytorch.loggers import CSVLogger
+from segger.prediction.predict_parquet import segment, load_model
+
+
+## VIASH START
+# Note: this section is auto-generated by viash at runtime. To edit it, make changes
+# in config.vsh.yaml and then run `viash config inject config.vsh.yaml`.
+par = {
+  'input_ist': 'resources_test/task_ist_preprocessing/mouse_brain_combined/raw_ist.zarr',
+  'input_segmentation': 'resources_test/task_ist_preprocessing/mouse_brain_combined/segmentation.zarr',
+  #'input_segmentation': 'temp/methods/cellpose/segmentation.zarr',
+  'transcripts_key': 'transcripts',
+  'coordinate_system': 'global',
+  'output': './temp/methods/segger/segger_assigned_transcripts.zarr',
+
+  #TODO: Add the interesting parameters: in this regard see:
+  # - https://elihei2.github.io/segger_dev/notebooks/segger_tutorial/#4-tune-parameters
+  # - 
+  #'force_2d': 'false',
+  #'min_molecules_per_cell': 50,
+  #'scale': -1.0,
+  #'scale_std': "25%",
+  #'n_clusters': 4,
+  #'prior_segmentation_confidence': 0.8,
+}
+meta = {
+  'name': 'segger_transcript_assignment',
+  'temp_dir': "./temp/methods/segger",
+  'cpus': 4,
+}
+## VIASH END
+
+TMP_DIR = Path(meta["temp_dir"] or "/tmp")
+TMP_DIR.mkdir(parents=True, exist_ok=True)
+
+
+#NOTE: for datafile preparation we follow 
+# https://github.com/EliHei2/segger_dev/blob/generic_config/platform_guides/platform_preparation_guide.ipynb and
+# https://github.com/EliHei2/segger_dev/blob/main/src/segger/data/parquet/_settings/xenium.yaml
+
+POLYGON_PARQUET = TMP_DIR / "nucleus_boundaries.parquet" 
+TRANSCRIPTS_PARQUET = TMP_DIR / "transcripts.parquet"
+SEGGER_DATA_DIR = TMP_DIR / 'data_segger'
+SEGGER_DATA_DIR.mkdir(parents=True, exist_ok=True)
+
+MODELS_DIR = TMP_DIR / 'models'
+MODELS_DIR.mkdir(parents=True, exist_ok=True)
+
+
+print('Checking if CUDA is available:', flush=True)
+print("\t torch.cuda.is_available() = ", torch.cuda.is_available())
+DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
+
+# Load data
+
+sdata = sd.read_zarr(par['input_ist'])
+sdata_segm = sd.read_zarr(par['input_segmentation'])
+
+######################
+# boundaries.parquet #
+#######################
+
+# Compute boundaries from segmentation
+print('Computing boundaries from segmentation', flush=True)
+boundaries = sd.to_polygons(sdata_segm['segmentation'])[['geometry']]
+
+# Bring boundaries into the coodinate system of the transcripts and save as parquet
+print('Transforming boundaries to transcripts coordinate system', flush=True)
+trans_segm_to_global = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']]
+trans_global_to_transcripts = sd.transformations.get_transformation(sdata[par['transcripts_key']], get_all=True)[par['coordinate_system']].inverse()
+trans = sequence = sd.transformations.Sequence([trans_segm_to_global, trans_global_to_transcripts])
+boundaries = sd.transform(boundaries, trans, par['coordinate_system'])
+boundaries.index.name = "cell_id"
+# Transform to dataframe into Xenium parquet format
+boundaries_df = boundaries['geometry'].get_coordinates().rename(columns={'x': 'vertex_x', 'y': 'vertex_y'})
+boundaries_df = boundaries_df.reset_index()
+boundaries_df['cell_id'] = boundaries_df['cell_id'].astype(str) + "_id"
+boundaries_df['cell_id'] = boundaries_df['cell_id'].astype('object')
+boundaries_df['vertex_x'] = boundaries_df['vertex_x'].astype(np.float32)
+boundaries_df['vertex_y'] = boundaries_df['vertex_y'].astype(np.float32)
+boundaries_df.to_parquet(POLYGON_PARQUET)
+del boundaries
+del boundaries_df
+
+#######################
+# transcripts.parquet #
+#######################
+
+# Map segmentation ids to transcripts (basic assignment)
+transcripts_coord_systems = sd.transformations.get_transformation(sdata[par["transcripts_key"]], get_all=True).keys()
+assert par['coordinate_system'] in transcripts_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
+segmentation_coord_systems = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True).keys()
+assert par['coordinate_system'] in segmentation_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
+
+print('Transforming transcripts coordinates', flush=True)
+transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
+
+# In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
+trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()
+transcripts = sd.transform(transcripts, trans, par['coordinate_system'])
+
+print('Assigning transcripts to cell ids', flush=True)
+y_coords = transcripts.y.compute().to_numpy(dtype=np.int64)
+x_coords = transcripts.x.compute().to_numpy(dtype=np.int64)
+if isinstance(sdata_segm["segmentation"], xr.DataTree):
+    label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() 
+else:
+    label_image = sdata_segm["segmentation"].to_numpy()
+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
+)
+transcripts_original = sdata[par['transcripts_key']]
+transcripts_original['cell_id'] = cell_id_dask_series
+transcripts_original["overlaps_nucleus"] = (transcripts_original['cell_id'] != 0)
+#transcripts_original["overlaps_nucleus"] = np.random.randint(0, 2, size=len(transcripts_original)) #TODO: remove this line
+transcripts_original["transcript_id"] = transcripts_original.index
+#TODO: The qv columns of Xenium should not be needed at the end. 
+#      See the issue at https://github.com/EliHei2/segger_dev/issues/125
+#transcripts_original["qv"] = 1.0 
+
+# Generate the same dtypes as in Xenium parquet files
+#transcripts_original['cell_id'] = transcripts_original['cell_id'].astype('object') # NOTE: this could probably work, but just to be sure we convert to strings
+#transcripts_original['cell_id'] = transcripts_original['cell_id'].replace('0', 'UNASSIGNED')
+transcripts_original['cell_id'] = transcripts_original['cell_id'].astype(str) + "_id"
+transcripts_original['cell_id'] = transcripts_original['cell_id'].replace('0_id', 'UNASSIGNED')
+transcripts_original['overlaps_nucleus'] = transcripts_original['overlaps_nucleus'].astype(np.uint8)
+#transcripts_original['qv'] = transcripts_original['qv'].astype(np.float32)
+transcripts_original['transcript_id'] = transcripts_original['transcript_id'].astype(np.uint64)
+
+# Rename columns to match Xenium parquet files
+rename_cols = {'x': 'x_location', 'y': 'y_location', 'z': 'z_location', 'feature_name': 'feature_name'}
+transcripts_original = transcripts_original.rename(columns=rename_cols)
+cols = list(rename_cols.values()) + ['cell_id', 'transcript_id', 'overlaps_nucleus'] #'qv', 
+
+# Convert to pandas dataframe and convert remaining dtypes to match Xenium parquet files
+df_transcripts = transcripts_original[cols].compute()
+df_transcripts['cell_id'] = df_transcripts['cell_id'].astype('object') # Didn't manage to convert to dtype 'object' above
+df_transcripts['feature_name'] = df_transcripts['feature_name'].astype('object')
+
+# Save transcripts to parquet
+df_transcripts.to_parquet(TRANSCRIPTS_PARQUET)
+
+del df_transcripts
+del transcripts_original
+del transcripts
+del y_coords
+del x_coords
+del label_image
+del cell_id_dask_series
+
+
+
+######################
+# Segger Data loader #
+######################
+
+#NOTE: This was just for debugging purposes with Xenium data, can be deleted
+#TMP_DIR2 = Path("temp/datasets/Xenium_V1_hSkin_nondiseased_section_1_FFPE_outs")
+#SEGGER_DATA_DIR2 = TMP_DIR / 'data_segger2'
+#sample = STSampleParquet(base_dir=TMP_DIR2, n_workers=4, sample_type='xenium')
+
+print('Preparing Segger data', flush=True)
+sample = STSampleParquet(base_dir=TMP_DIR, n_workers=4, sample_type='xenium')
+
+sample.save(
+      data_dir=SEGGER_DATA_DIR,
+      k_bd=3,
+      dist_bd=15.0,
+      k_tx=3,
+      dist_tx=5.0,
+      tile_width=120,
+      tile_height=120,
+      neg_sampling_ratio=5.0,
+      frac=1.0,
+      val_prob=0.1,
+      test_prob=0.2,
+)
+
+
+#TODO: optionally add scRNA-seq based embedding https://elihei2.github.io/segger_dev/notebooks/segger_tutorial/#12-using-custom-gene-embeddings
+#      note that the arguments is_token_based and num_tx_tokens need to change in case of such embedding (see below)
+
+
+######################
+# Train Segger model #
+######################
+
+# Base directory to store Pytorch Lightning models
+
+# Initialize the Lightning data module
+print('Initializing Segger data module', flush=True)
+dm = SeggerDataModule(
+    data_dir=SEGGER_DATA_DIR,
+    batch_size=2,
+    num_workers=meta["cpus"] or 1,
+)
+dm.setup()
+
+#is_token_based = True
+num_tx_tokens = 500
+
+# If you use custom (scRNA-seq based) gene embeddings, use the following two lines instead:
+# is_token_based = False
+# num_tx_tokens = dm.train[0].x_dict["tx"].shape[1] # Set the number of tokens to the number of genes
+
+# Initialize the Lightning model
+#NOTE: This part throws an error:
+# deviates from the current documentation, we follow https://github.com/EliHei2/segger_dev/blob/main/docs/notebooks/segger_tutorial.ipynb
+print('Initializing Segger model', flush=True)
+###ls = LitSegger(
+###    #is_token_based = is_token_based,
+###    #num_node_features = {"tx": num_tx_tokens, "bd": num_bd_features},
+###    #init_emb=8,    
+###    #hidden_channels=64,
+###    out_channels=16,
+###    heads=4,
+###    num_mid_layers=1,
+###    aggr='sum',
+###)
+
+model = Segger(
+    # is_token_based=is_token_based,
+    num_tx_tokens=num_tx_tokens,
+    init_emb=8,
+    hidden_channels=64,
+    out_channels=16,
+    heads=4,
+    num_mid_layers=3,
+)
+model = to_hetero(model, (["tx", "bd"], [("tx", "belongs", "bd"), ("tx", "neighbors", "tx")]), aggr="sum")
+
+batch = dm.train[0]
+model.forward(batch.x_dict, batch.edge_index_dict)
+# Wrap the model in LitSegger
+ls = LitSegger(model=model)
+
+
+# Initialize the Lightning trainer
+print('Training Segger model', flush=True)
+trainer = Trainer(
+    accelerator=DEVICE,
+    strategy='auto',
+    precision='16-mixed',
+    devices=1, # set higher number if more gpus are available
+    max_epochs=2,#100,
+    default_root_dir=MODELS_DIR,
+    logger=CSVLogger(MODELS_DIR),
+)
+
+# Fit model
+trainer.fit(
+    model=ls,
+    datamodule=dm
+)
+
+# For debugging:
+"""
+# Evaluate results
+model_version = 0  # 'v_num' from training output above
+model_path = MODELS_DIR / 'lightning_logs' / f'version_{model_version}'
+metrics = pd.read_csv(model_path / 'metrics.csv', index_col=1)
+
+fig, ax = plt.subplots(1,1, figsize=(2,2))
+
+for col in metrics.columns.difference(['epoch']):
+    metric = metrics[col].dropna()
+    ax.plot(metric.index, metric.values, label=col)
+
+ax.legend(loc=(1, 0.33))
+ax.set_ylim(0, 1)
+ax.set_xlabel('Step')
+"""
+
+#############
+# Inference #
+#############
+
+print('Loading Segger model', flush=True)
+model_version = 0
+model_path = MODELS_DIR / "lightning_logs" / f"version_{model_version}"
+model = load_model(model_path / "checkpoints")
+
+#NOTE: Parameters selected according to https://github.com/EliHei2/segger_dev/issues/126#issuecomment-3160389883
+receptive_field = {'k_bd': 4, 'dist_bd': 7.5, 'k_tx': 15, 'dist_tx': 3}
+
+print('Running Segger inference', flush=True)
+#NOTE: Running into this error (also with the latest segger version, 07.08.2025): https://github.com/EliHei2/segger_dev/issues/126
+segment(
+    model,
+    dm,
+    score_cut = .75, 
+    save_dir=TMP_DIR,
+    seg_tag='segger_output',
+    transcript_file=TRANSCRIPTS_PARQUET,
+    receptive_field=receptive_field,
+    min_transcripts=5,
+    cell_id_col='segger_cell_id',
+    use_cc=False,
+    knn_method='kd_tree',
+    verbose=True,
+)
+
+##########################################################
+# Load Segger results and save transcripts with cell ids #
+##########################################################
+
+#TODO: Look at other transcript assignment methods for these last steps
+
+#TODO: Check if this lines are correct
+#
+#sdata = sd.read_zarr(par['input_ist'])
+#segger_df = gpd.read_parquet(TMP_DIR / "segger_output.parquet")
+#
+#transcripts = sdata[par['transcripts_key']]
+#assert transcripts.index == segger_df.index, "Transcripts and Segger results have different indices"
+#transcripts['cell_id'] = segger_df['segger_cell_id']
+
+
+#TODO: Save transcripts with cell ids
+
+
+

From d836587318876f03ef0335db984c2edf816ce543 Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Sun, 19 Jul 2026 12:33:37 +0200
Subject: [PATCH 4/5] Remove segger draft from fixes (moved to segger_newtry
 branch)

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 .../segger/config.vsh.yaml                    | 101 ------
 .../segger/script.py                          | 338 ------------------
 2 files changed, 439 deletions(-)
 delete mode 100644 src/methods_transcript_assignment/segger/config.vsh.yaml
 delete mode 100644 src/methods_transcript_assignment/segger/script.py

diff --git a/src/methods_transcript_assignment/segger/config.vsh.yaml b/src/methods_transcript_assignment/segger/config.vsh.yaml
deleted file mode 100644
index 901282c78..000000000
--- a/src/methods_transcript_assignment/segger/config.vsh.yaml
+++ /dev/null
@@ -1,101 +0,0 @@
-__merge__: /src/api/comp_method_transcript_assignment.yaml
-
-name: segger
-label: "Segger Transcript Assignment"
-summary: "Assign transcripts to cells using the Segger method"
-description: "Segger is a tool for cell segmentation in single-molecule spatial omics datasets that leverages graph neural networks (GNNs) and heterogeneous graphs."
-links:
-  documentation: "https://elihei2.github.io/segger_dev/"
-  repository: "https://github.com/EliHei2/segger_dev"
-references:
-  doi: "10.1101/2025.03.14.643160"
-
-arguments:
-  - name: --transcripts_key
-    type: string
-    description: The key of the transcripts within the points of the spatial data
-    default: transcripts
-  - name: --coordinate_system
-    type: string
-    description: The key of the pixel space coordinate system within the spatial data
-    default: global
-
-#  - name: --force_2d
-#    type: string
-#    required: false
-#    description: "Ignores z-column in the data if it is provided"
-#    direction: input 
-#    default: "false"
-#
-#  - name: --min_molecules_per_cell
-#    type: integer
-#    required: false
-#    description: "Minimal number of molecules per cell"
-#    direction: input 
-#    default: 50
-#  
-#  - name: --scale
-#    type: double
-#    required: false
-#    description: |
-#      "Scale parameter, which suggest approximate cell radius for the algorithm. Must be in the same units as 
-#      x and y molecule coordinates. Negative values mean it must be estimated from `min_molecules_per_cell`."
-#    direction: input 
-#    default: -1.0
-#
-#  - name: --scale_std
-#    type: string
-#    required: false
-#    description: "Standard deviation of scale across cells relative to `scale`"
-#    direction: input 
-#    default: "25%"
-#
-#  - name: --n_clusters
-#    type: integer
-#    required: false
-#    description: "Number of molecule clusters, i.e. major cell types."
-#    direction: input 
-#    default: 4
-#
-#  - name: --prior_segmentation_confidence
-#    type: double
-#    required: false
-#    description: "Confidence of the prior segmentation"
-#    direction: input 
-#    default: 0.8
-
-resources:
-  - type: python_script
-    path: script.py
-
-engines:
-#  - type: docker
-#    image: openproblems/base_python:1
-#    __merge__: 
-#      - /src/base/setup_spatialdata_partial.yaml
-#    setup:
-#      - type: python
-#        pypi: [torch, pytorch-lightning, lightning]
-#      - type: docker
-#        run:
-#          - git clone https://github.com/EliHei2/segger_dev.git
-#          - cd segger_dev && pip install .[cuda12]
-engines:
-  - type: docker
-    image: danielunyi42/segger_dev:cuda121
-    setup:
-      - type: apt
-        packages: procps
-      - type: python
-        github:
-          - openproblems-bio/core#subdirectory=packages/python/openproblems
-      - type: python
-        packages:
-          - spatialdata
-  - type: native
-
-runners:
-  - type: executable
-  - type: nextflow
-    directives:
-      label: [ midtime, midcpu, midmem, gpu ]
diff --git a/src/methods_transcript_assignment/segger/script.py b/src/methods_transcript_assignment/segger/script.py
deleted file mode 100644
index 84cf7e976..000000000
--- a/src/methods_transcript_assignment/segger/script.py
+++ /dev/null
@@ -1,338 +0,0 @@
-from pathlib import Path
-import dask
-import torch
-import xarray as xr
-import numpy as np
-import geopandas as gpd
-import spatialdata as sd
-from segger.data.parquet.sample import STSampleParquet
-from segger.training.segger_data_module import SeggerDataModule
-from segger.training.train import LitSegger, Segger
-from torch_geometric.nn import to_hetero
-from lightning.pytorch import Trainer
-#from pytorch_lightning import Trainer
-from lightning.pytorch.loggers import CSVLogger
-from segger.prediction.predict_parquet import segment, load_model
-
-
-## VIASH START
-# Note: this section is auto-generated by viash at runtime. To edit it, make changes
-# in config.vsh.yaml and then run `viash config inject config.vsh.yaml`.
-par = {
-  'input_ist': 'resources_test/task_ist_preprocessing/mouse_brain_combined/raw_ist.zarr',
-  'input_segmentation': 'resources_test/task_ist_preprocessing/mouse_brain_combined/segmentation.zarr',
-  #'input_segmentation': 'temp/methods/cellpose/segmentation.zarr',
-  'transcripts_key': 'transcripts',
-  'coordinate_system': 'global',
-  'output': './temp/methods/segger/segger_assigned_transcripts.zarr',
-
-  #TODO: Add the interesting parameters: in this regard see:
-  # - https://elihei2.github.io/segger_dev/notebooks/segger_tutorial/#4-tune-parameters
-  # - 
-  #'force_2d': 'false',
-  #'min_molecules_per_cell': 50,
-  #'scale': -1.0,
-  #'scale_std': "25%",
-  #'n_clusters': 4,
-  #'prior_segmentation_confidence': 0.8,
-}
-meta = {
-  'name': 'segger_transcript_assignment',
-  'temp_dir': "./temp/methods/segger",
-  'cpus': 4,
-}
-## VIASH END
-
-TMP_DIR = Path(meta["temp_dir"] or "/tmp")
-TMP_DIR.mkdir(parents=True, exist_ok=True)
-
-
-#NOTE: for datafile preparation we follow 
-# https://github.com/EliHei2/segger_dev/blob/generic_config/platform_guides/platform_preparation_guide.ipynb and
-# https://github.com/EliHei2/segger_dev/blob/main/src/segger/data/parquet/_settings/xenium.yaml
-
-POLYGON_PARQUET = TMP_DIR / "nucleus_boundaries.parquet" 
-TRANSCRIPTS_PARQUET = TMP_DIR / "transcripts.parquet"
-SEGGER_DATA_DIR = TMP_DIR / 'data_segger'
-SEGGER_DATA_DIR.mkdir(parents=True, exist_ok=True)
-
-MODELS_DIR = TMP_DIR / 'models'
-MODELS_DIR.mkdir(parents=True, exist_ok=True)
-
-
-print('Checking if CUDA is available:', flush=True)
-print("\t torch.cuda.is_available() = ", torch.cuda.is_available())
-DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
-
-# Load data
-
-sdata = sd.read_zarr(par['input_ist'])
-sdata_segm = sd.read_zarr(par['input_segmentation'])
-
-######################
-# boundaries.parquet #
-#######################
-
-# Compute boundaries from segmentation
-print('Computing boundaries from segmentation', flush=True)
-boundaries = sd.to_polygons(sdata_segm['segmentation'])[['geometry']]
-
-# Bring boundaries into the coodinate system of the transcripts and save as parquet
-print('Transforming boundaries to transcripts coordinate system', flush=True)
-trans_segm_to_global = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']]
-trans_global_to_transcripts = sd.transformations.get_transformation(sdata[par['transcripts_key']], get_all=True)[par['coordinate_system']].inverse()
-trans = sequence = sd.transformations.Sequence([trans_segm_to_global, trans_global_to_transcripts])
-boundaries = sd.transform(boundaries, trans, par['coordinate_system'])
-boundaries.index.name = "cell_id"
-# Transform to dataframe into Xenium parquet format
-boundaries_df = boundaries['geometry'].get_coordinates().rename(columns={'x': 'vertex_x', 'y': 'vertex_y'})
-boundaries_df = boundaries_df.reset_index()
-boundaries_df['cell_id'] = boundaries_df['cell_id'].astype(str) + "_id"
-boundaries_df['cell_id'] = boundaries_df['cell_id'].astype('object')
-boundaries_df['vertex_x'] = boundaries_df['vertex_x'].astype(np.float32)
-boundaries_df['vertex_y'] = boundaries_df['vertex_y'].astype(np.float32)
-boundaries_df.to_parquet(POLYGON_PARQUET)
-del boundaries
-del boundaries_df
-
-#######################
-# transcripts.parquet #
-#######################
-
-# Map segmentation ids to transcripts (basic assignment)
-transcripts_coord_systems = sd.transformations.get_transformation(sdata[par["transcripts_key"]], get_all=True).keys()
-assert par['coordinate_system'] in transcripts_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
-segmentation_coord_systems = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True).keys()
-assert par['coordinate_system'] in segmentation_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data."
-
-print('Transforming transcripts coordinates', flush=True)
-transcripts = sd.transform(sdata[par['transcripts_key']], to_coordinate_system=par['coordinate_system'])
-
-# In case of a translation transformation of the segmentation (e.g. crop of the data), we need to adjust the transcript coordinates
-trans = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True)[par['coordinate_system']].inverse()
-transcripts = sd.transform(transcripts, trans, par['coordinate_system'])
-
-print('Assigning transcripts to cell ids', flush=True)
-y_coords = transcripts.y.compute().to_numpy(dtype=np.int64)
-x_coords = transcripts.x.compute().to_numpy(dtype=np.int64)
-if isinstance(sdata_segm["segmentation"], xr.DataTree):
-    label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() 
-else:
-    label_image = sdata_segm["segmentation"].to_numpy()
-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
-)
-transcripts_original = sdata[par['transcripts_key']]
-transcripts_original['cell_id'] = cell_id_dask_series
-transcripts_original["overlaps_nucleus"] = (transcripts_original['cell_id'] != 0)
-#transcripts_original["overlaps_nucleus"] = np.random.randint(0, 2, size=len(transcripts_original)) #TODO: remove this line
-transcripts_original["transcript_id"] = transcripts_original.index
-#TODO: The qv columns of Xenium should not be needed at the end. 
-#      See the issue at https://github.com/EliHei2/segger_dev/issues/125
-#transcripts_original["qv"] = 1.0 
-
-# Generate the same dtypes as in Xenium parquet files
-#transcripts_original['cell_id'] = transcripts_original['cell_id'].astype('object') # NOTE: this could probably work, but just to be sure we convert to strings
-#transcripts_original['cell_id'] = transcripts_original['cell_id'].replace('0', 'UNASSIGNED')
-transcripts_original['cell_id'] = transcripts_original['cell_id'].astype(str) + "_id"
-transcripts_original['cell_id'] = transcripts_original['cell_id'].replace('0_id', 'UNASSIGNED')
-transcripts_original['overlaps_nucleus'] = transcripts_original['overlaps_nucleus'].astype(np.uint8)
-#transcripts_original['qv'] = transcripts_original['qv'].astype(np.float32)
-transcripts_original['transcript_id'] = transcripts_original['transcript_id'].astype(np.uint64)
-
-# Rename columns to match Xenium parquet files
-rename_cols = {'x': 'x_location', 'y': 'y_location', 'z': 'z_location', 'feature_name': 'feature_name'}
-transcripts_original = transcripts_original.rename(columns=rename_cols)
-cols = list(rename_cols.values()) + ['cell_id', 'transcript_id', 'overlaps_nucleus'] #'qv', 
-
-# Convert to pandas dataframe and convert remaining dtypes to match Xenium parquet files
-df_transcripts = transcripts_original[cols].compute()
-df_transcripts['cell_id'] = df_transcripts['cell_id'].astype('object') # Didn't manage to convert to dtype 'object' above
-df_transcripts['feature_name'] = df_transcripts['feature_name'].astype('object')
-
-# Save transcripts to parquet
-df_transcripts.to_parquet(TRANSCRIPTS_PARQUET)
-
-del df_transcripts
-del transcripts_original
-del transcripts
-del y_coords
-del x_coords
-del label_image
-del cell_id_dask_series
-
-
-
-######################
-# Segger Data loader #
-######################
-
-#NOTE: This was just for debugging purposes with Xenium data, can be deleted
-#TMP_DIR2 = Path("temp/datasets/Xenium_V1_hSkin_nondiseased_section_1_FFPE_outs")
-#SEGGER_DATA_DIR2 = TMP_DIR / 'data_segger2'
-#sample = STSampleParquet(base_dir=TMP_DIR2, n_workers=4, sample_type='xenium')
-
-print('Preparing Segger data', flush=True)
-sample = STSampleParquet(base_dir=TMP_DIR, n_workers=4, sample_type='xenium')
-
-sample.save(
-      data_dir=SEGGER_DATA_DIR,
-      k_bd=3,
-      dist_bd=15.0,
-      k_tx=3,
-      dist_tx=5.0,
-      tile_width=120,
-      tile_height=120,
-      neg_sampling_ratio=5.0,
-      frac=1.0,
-      val_prob=0.1,
-      test_prob=0.2,
-)
-
-
-#TODO: optionally add scRNA-seq based embedding https://elihei2.github.io/segger_dev/notebooks/segger_tutorial/#12-using-custom-gene-embeddings
-#      note that the arguments is_token_based and num_tx_tokens need to change in case of such embedding (see below)
-
-
-######################
-# Train Segger model #
-######################
-
-# Base directory to store Pytorch Lightning models
-
-# Initialize the Lightning data module
-print('Initializing Segger data module', flush=True)
-dm = SeggerDataModule(
-    data_dir=SEGGER_DATA_DIR,
-    batch_size=2,
-    num_workers=meta["cpus"] or 1,
-)
-dm.setup()
-
-#is_token_based = True
-num_tx_tokens = 500
-
-# If you use custom (scRNA-seq based) gene embeddings, use the following two lines instead:
-# is_token_based = False
-# num_tx_tokens = dm.train[0].x_dict["tx"].shape[1] # Set the number of tokens to the number of genes
-
-# Initialize the Lightning model
-#NOTE: This part throws an error:
-# deviates from the current documentation, we follow https://github.com/EliHei2/segger_dev/blob/main/docs/notebooks/segger_tutorial.ipynb
-print('Initializing Segger model', flush=True)
-###ls = LitSegger(
-###    #is_token_based = is_token_based,
-###    #num_node_features = {"tx": num_tx_tokens, "bd": num_bd_features},
-###    #init_emb=8,    
-###    #hidden_channels=64,
-###    out_channels=16,
-###    heads=4,
-###    num_mid_layers=1,
-###    aggr='sum',
-###)
-
-model = Segger(
-    # is_token_based=is_token_based,
-    num_tx_tokens=num_tx_tokens,
-    init_emb=8,
-    hidden_channels=64,
-    out_channels=16,
-    heads=4,
-    num_mid_layers=3,
-)
-model = to_hetero(model, (["tx", "bd"], [("tx", "belongs", "bd"), ("tx", "neighbors", "tx")]), aggr="sum")
-
-batch = dm.train[0]
-model.forward(batch.x_dict, batch.edge_index_dict)
-# Wrap the model in LitSegger
-ls = LitSegger(model=model)
-
-
-# Initialize the Lightning trainer
-print('Training Segger model', flush=True)
-trainer = Trainer(
-    accelerator=DEVICE,
-    strategy='auto',
-    precision='16-mixed',
-    devices=1, # set higher number if more gpus are available
-    max_epochs=2,#100,
-    default_root_dir=MODELS_DIR,
-    logger=CSVLogger(MODELS_DIR),
-)
-
-# Fit model
-trainer.fit(
-    model=ls,
-    datamodule=dm
-)
-
-# For debugging:
-"""
-# Evaluate results
-model_version = 0  # 'v_num' from training output above
-model_path = MODELS_DIR / 'lightning_logs' / f'version_{model_version}'
-metrics = pd.read_csv(model_path / 'metrics.csv', index_col=1)
-
-fig, ax = plt.subplots(1,1, figsize=(2,2))
-
-for col in metrics.columns.difference(['epoch']):
-    metric = metrics[col].dropna()
-    ax.plot(metric.index, metric.values, label=col)
-
-ax.legend(loc=(1, 0.33))
-ax.set_ylim(0, 1)
-ax.set_xlabel('Step')
-"""
-
-#############
-# Inference #
-#############
-
-print('Loading Segger model', flush=True)
-model_version = 0
-model_path = MODELS_DIR / "lightning_logs" / f"version_{model_version}"
-model = load_model(model_path / "checkpoints")
-
-#NOTE: Parameters selected according to https://github.com/EliHei2/segger_dev/issues/126#issuecomment-3160389883
-receptive_field = {'k_bd': 4, 'dist_bd': 7.5, 'k_tx': 15, 'dist_tx': 3}
-
-print('Running Segger inference', flush=True)
-#NOTE: Running into this error (also with the latest segger version, 07.08.2025): https://github.com/EliHei2/segger_dev/issues/126
-segment(
-    model,
-    dm,
-    score_cut = .75, 
-    save_dir=TMP_DIR,
-    seg_tag='segger_output',
-    transcript_file=TRANSCRIPTS_PARQUET,
-    receptive_field=receptive_field,
-    min_transcripts=5,
-    cell_id_col='segger_cell_id',
-    use_cc=False,
-    knn_method='kd_tree',
-    verbose=True,
-)
-
-##########################################################
-# Load Segger results and save transcripts with cell ids #
-##########################################################
-
-#TODO: Look at other transcript assignment methods for these last steps
-
-#TODO: Check if this lines are correct
-#
-#sdata = sd.read_zarr(par['input_ist'])
-#segger_df = gpd.read_parquet(TMP_DIR / "segger_output.parquet")
-#
-#transcripts = sdata[par['transcripts_key']]
-#assert transcripts.index == segger_df.index, "Transcripts and Segger results have different indices"
-#transcripts['cell_id'] = segger_df['segger_cell_id']
-
-
-#TODO: Save transcripts with cell ids
-
-
-

From 3b9a1fb8c28201c7be15f6d413c01a78ab21d50f Mon Sep 17 00:00:00 2001
From: dariarom94 <romanovskaiadaria@gmail.com>
Date: Sun, 19 Jul 2026 12:51:19 +0200
Subject: [PATCH 5/5] adjusting stardist versions

---
 src/methods_segmentation/stardist/config.vsh.yaml | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/methods_segmentation/stardist/config.vsh.yaml b/src/methods_segmentation/stardist/config.vsh.yaml
index d361e5d9d..0330b270a 100644
--- a/src/methods_segmentation/stardist/config.vsh.yaml
+++ b/src/methods_segmentation/stardist/config.vsh.yaml
@@ -31,9 +31,15 @@ engines:
     setup:
       - type: python
         pypi:
+          # TensorFlow must be >=2.18 so it allows numpy>=2.0. anndata/spatialdata
+          # reference np.bool (readded in numpy 2.0), so a numpy<2.0 pin makes
+          # `import anndata` crash at runtime; TF 2.17 caps numpy<2.0, so we bump TF.
           - stardist
-          - tensorflow==2.17.0
-          - "numpy>=1.26.0,<2.0.0"
+          - tensorflow==2.18.0
+          # csbdeep/stardist use the Keras 2 API via tf_keras, which must match
+          # the tensorflow version (the base image ships tf_keras 2.17).
+          - tf_keras==2.18.0
+          - "numpy>=2.0.0,<2.2.0"
           - scipy<1.15.0
   - type: native