Skip to content
Merged

Fixes #167

Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 20 additions & 9 deletions src/base/labels_nebius.config
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,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 }
Expand All @@ -58,17 +63,23 @@ process {
disk = { 100.GB * task.attempt }
pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00wvqfyde1nfwezs0']]
}
// 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 {
// Nebius 128vcpu-512gb nodes: 503 GB allocatable — cap well below node max
memory = { task.attempt == 1 ? 200.GB : 450.GB }
disk = 200.GB
pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00hnqhyfdcsy9m09n']]
memory = { get_memory( 200.GB * task.attempt ) }
disk = { 200.GB * task.attempt }
}
// 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 {
// Nebius 128vcpu-512gb nodes: 503 GB allocatable — cap at 450 GB to leave headroom
memory = { [ 200.GB * task.attempt, 450.GB ].min() }
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)}" : ""}
Expand Down Expand Up @@ -170,7 +181,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
Expand Down
31 changes: 31 additions & 0 deletions src/methods_cell_type_annotation/rctd/config.vsh.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
18 changes: 16 additions & 2 deletions src/methods_cell_type_annotation/rctd/script.R
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/methods_segmentation/binning/config.vsh.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -31,4 +31,4 @@ runners:
- type: executable
- type: nextflow
directives:
label: [ midtime, midcpu, highmem ]
label: [ midtime, midcpu, midmem ]
20 changes: 18 additions & 2 deletions src/methods_segmentation/binning/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
7 changes: 7 additions & 0 deletions src/methods_segmentation/cellpose/config.vsh.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 16 additions & 4 deletions src/methods_segmentation/custom_segmentation/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -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[[]]
)
}
Expand Down
10 changes: 9 additions & 1 deletion src/methods_segmentation/stardist/config.vsh.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,16 @@ 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
- 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

runners:
Expand Down
31 changes: 29 additions & 2 deletions src/methods_segmentation/stardist/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,17 @@
import shutil
from pathlib import Path
import numpy as np
# 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
Expand Down Expand Up @@ -77,9 +88,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)
Expand Down
2 changes: 1 addition & 1 deletion src/methods_segmentation/watershed/config.vsh.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -173,5 +173,5 @@ runners:
- type: executable
- type: nextflow
directives:
label: [ hightime, midcpu, highmem ]
label: [ hightime, midcpu, midmem ]

20 changes: 18 additions & 2 deletions src/methods_segmentation/watershed/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
11 changes: 10 additions & 1 deletion src/methods_transcript_assignment/baysor/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
11 changes: 10 additions & 1 deletion src/methods_transcript_assignment/clustermap/script.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
Loading
Loading