Fixes

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

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

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?

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

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(

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

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)
 

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
     },
 )
 

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(

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

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',

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
     },
 )