diff --git a/scripts/create_resources/spatial/process_bruker_cosmx.sh b/scripts/create_resources/spatial/process_bruker_cosmx.sh index 7e4a7ea29..a303b0997 100644 --- a/scripts/create_resources/spatial/process_bruker_cosmx.sh +++ b/scripts/create_resources/spatial/process_bruker_cosmx.sh @@ -14,8 +14,8 @@ cat > /tmp/params.yaml << HERE param_list: - id: "bruker_cosmx/bruker_mouse_brain_cosmx/rep1" - input_raw: "https://smi-public.objects.liquidweb.services/HalfBrain.zip" - input_flat_files: "https://smi-public.objects.liquidweb.services/Half%20%20Brain%20simple%20%20files%20.zip" + input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/HalfBrain.zip" + input_flat_files: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Half Brain simple files.zip" dataset_name: "Bruker CosMx Mouse Brain" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/cosmx-smi-mouse-brain-ffpe-dataset/" dataset_summary: "Bruker CosMx Mouse Brain dataset on FFPE covering a full hemisphere of a mouse brain." @@ -24,7 +24,7 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_liver_cosmx" - input_raw: "https://smi-public.objects.liquidweb.services/NormalLiverFiles.zip" + input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/NormalLiverFiles.zip" dataset_name: "Bruker CosMx Human Liver" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Liver dataset on FFPE." diff --git a/scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh b/scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh index 9c00dfcb7..c12a28b1b 100644 --- a/scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh +++ b/scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh @@ -14,8 +14,8 @@ cat > /tmp/params.yaml << HERE param_list: - id: "bruker_cosmx/bruker_mouse_brain_cosmx/rep1" - input_raw: "https://smi-public.objects.liquidweb.services/HalfBrain.zip" - input_flat_files: "https://smi-public.objects.liquidweb.services/Half%20%20Brain%20simple%20%20files%20.zip" + input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/HalfBrain.zip" + input_flat_files: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Half Brain simple files.zip" dataset_name: "Bruker CosMx Mouse Brain" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/cosmx-smi-mouse-brain-ffpe-dataset/" dataset_summary: "Bruker CosMx Mouse Brain dataset on FFPE covering a full hemisphere of a mouse brain." @@ -24,7 +24,7 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_liver_cosmx" - input_raw: "https://smi-public.objects.liquidweb.services/NormalLiverFiles.zip" + input_raw: "s3://openproblems-data/resources/raw_data/bruker_cosmx/NormalLiverFiles.zip" dataset_name: "Bruker CosMx Human Liver" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Liver dataset on FFPE." diff --git a/src/methods_cell_type_annotation/moscot/script.py b/src/methods_cell_type_annotation/moscot/script.py index ac0b7343b..26be68189 100644 --- a/src/methods_cell_type_annotation/moscot/script.py +++ b/src/methods_cell_type_annotation/moscot/script.py @@ -75,13 +75,13 @@ adata_sp.X = adata_sp.layers["normalized"] adata_sp.obsm["spatial"] = adata_sp.obs[["centroid_x", "centroid_y"]].to_numpy() -sc.pp.pca(adata_sc, n_comps=50) # X is the normalized layer set above +sc.pp.pca(adata_sc, n_comps=30) # X is the normalized layer set above # Define mapping problem mp = MappingProblem(adata_sc=adata_sc, adata_sp=adata_sp) mp = mp.prepare( - sc_attr={"attr": "obsm", "key": "X_pca"}, # <-- 50-dim, not raw genes + sc_attr={"attr": "obsm", "key": "X_pca"}, # <-- 30-dim, not raw genes xy_callback="local-pca", ) diff --git a/src/methods_transcript_assignment/comseg/config.vsh.yaml b/src/methods_transcript_assignment/comseg/config.vsh.yaml index 000e363f8..a92907f43 100644 --- a/src/methods_transcript_assignment/comseg/config.vsh.yaml +++ b/src/methods_transcript_assignment/comseg/config.vsh.yaml @@ -63,6 +63,30 @@ arguments: type: boolean default: true description: "Allow disconnected polygons in segmentation" + - name: --n_workers + type: integer + default: 4 + description: | + Maximum number of dask worker processes used to run ComSeg patches in + parallel. ComSeg builds an in-memory graph per patch, so peak memory + scales with the number of concurrent workers; the value is clamped to + (cpus - 1). Lower this if the job is OOM-killed (exit 137); raise it + (up to the CPU count) to run faster when memory allows. + - name: --worker_memory_limit + type: string + default: "0" + description: | + Per-worker memory limit for the dask LocalCluster. Default "0" (also + "none"/"off"/"disabled") neutralises dask's per-worker memory monitor by + setting each worker's limit to the FULL detected container memory. ComSeg's + per-patch memory is dominated by unmanaged native allocations + (numba/shapely/graphs) that dask cannot spill, so a tighter limit does not + buy spilling — it only makes the nanny kill otherwise-healthy workers + (KilledWorker) once they cross ~95% of the limit. Total memory is instead + bounded by n_workers, and a genuine OOM is caught by the container cgroup + + Nextflow retry. Set an explicit value (e.g. "20GB") only if you want dask + to police each worker; "auto" is discouraged because it divides detected + RAM by CPU count and is easily smaller than a single patch's working set. resources: diff --git a/src/methods_transcript_assignment/comseg/script.py b/src/methods_transcript_assignment/comseg/script.py index 2829b873f..1905fe29b 100644 --- a/src/methods_transcript_assignment/comseg/script.py +++ b/src/methods_transcript_assignment/comseg/script.py @@ -4,6 +4,7 @@ import anndata as ad import pandas as pd import numpy as np +import dask import dask.dataframe as dd ## VIASH START @@ -24,6 +25,8 @@ "gene_column": "feature_name", "norm_vector": False, "allow_disconnected_polygon": True, + "n_workers": 4, + "worker_memory_limit": "0", } meta = { "name": "comseg", @@ -79,12 +82,54 @@ # ComSeg processes each transcript patch independently and is pure-Python, so it # runs single-threaded unless a parallelization backend is enabled. Use the dask -# backend to spread patches across the allocated CPUs (see midcpu/highmem labels). -n_workers = max((meta["cpus"] or os.cpu_count() or 1) - 1, 1) +# backend to spread patches across workers, but CAP concurrency: comseg builds an +# in-memory graph per patch, so peak memory scales with the number of concurrent +# workers. Running one worker per CPU (14 on midcpu) multiplied peak RAM ~14x and +# got the job OOM-killed (exit 137 = cgroup SIGKILL). Bounding the worker count is +# the real memory safeguard (the per-worker memory monitor is neutralised below, +# see the note there). +cpu_cap = max((meta["cpus"] or os.cpu_count() or 1) - 1, 1) +n_workers = max(min(par["n_workers"], cpu_cap), 1) + +# `distributed` defaults worker startup to "spawn", which re-imports this script +# in every worker. As a plain viash script it has no `if __name__ == "__main__"` +# guard, so the re-import re-runs the module-level code (creating yet another +# cluster) and Python's spawn bootstrap check aborts with "An attempt has been +# made to start a new process before ... bootstrapping". Force "fork" so workers +# inherit the already-initialised interpreter instead of re-importing it; this is +# the Linux-native default and the mode comseg ran under before the deps bumped. +dask.config.set({"distributed.worker.multiprocessing-method": "fork"}) + sopa.settings.parallelization_backend = "dask" sopa.settings.dask_client_kwargs["n_workers"] = n_workers sopa.settings.dask_client_kwargs["threads_per_worker"] = 1 # CPU-bound work, avoid GIL contention -print(f"Running ComSeg with dask backend, n_workers={n_workers}", flush=True) + +# ComSeg's per-patch memory is mostly UNMANAGED (native numba/shapely/graph +# allocations, plus the forked interpreter's copy-on-write pages), which dask +# cannot spill. A per-worker memory_limit therefore never spills here; it only +# lets the nanny kill workers that cross ~95% of the limit -> KilledWorker, even +# though nothing is leaking. dask's "auto" is worse: it divides detected RAM by +# the CPU count (not n_workers), so on an 8-core box each worker gets ~1/8 of RAM +# and a normal ~1 GB patch already trips 95%. We therefore neutralise the monitor +# by default, but sopa reads worker.memory_manager.memory_limit and compares it to +# 4 GiB, so the value must be a real number (0/None makes sopa crash). Setting the +# limit to the FULL detected container memory does both: sopa sees a number, and +# the nanny only ever fires near 95% of the whole cgroup — where the OS + Nextflow +# retry already take over — never prematurely for a single patch. Total memory is +# bounded by n_workers instead. An explicit value (e.g. "20GB") is honoured as-is. +_mem_limit_arg = (par["worker_memory_limit"] or "").strip().lower() +if _mem_limit_arg in ("", "0", "none", "off", "disabled"): + from distributed.system import MEMORY_LIMIT as _CONTAINER_MEMORY + worker_memory_limit = _CONTAINER_MEMORY +else: + worker_memory_limit = par["worker_memory_limit"] +sopa.settings.dask_client_kwargs["memory_limit"] = worker_memory_limit + +print( + f"Running ComSeg with dask backend, n_workers={n_workers} " + f"(cap {cpu_cap}), memory_limit={worker_memory_limit!r}", + flush=True, +) sopa.segmentation.comseg(sdata, config) diff --git a/src/methods_transcript_assignment/fastreseg/config.vsh.yaml b/src/methods_transcript_assignment/fastreseg/config.vsh.yaml index 512828392..d5c9985ec 100644 --- a/src/methods_transcript_assignment/fastreseg/config.vsh.yaml +++ b/src/methods_transcript_assignment/fastreseg/config.vsh.yaml @@ -46,24 +46,45 @@ engines: - gdal-bin - libgdal-dev - r-base - - type: r - packages: - - devtools - - terra - - Matrix - - dbscan - - igraph - - matrixStats - - codetools - - data.table - # Add other CRAN packages here - github: - # Add GitHub packages here if needed - - Nanostring-Biostats/FastReseg - #- type: python - # github: - # - openproblems-bio/core#subdirectory=packages/python/openproblems - __merge__: + # FastReseg's R dependencies as prebuilt Debian binaries. Installing + # these via apt avoids compiling terra/igraph/sf/spatstat/stringi/etc. + # from source (extremely slow, especially under amd64 emulation on + # arm64 hosts). Only concaveman + GiottoUtils/GiottoClass + FastReseg + # itself still build from source below; their heavy deps are already + # present as binaries here. + - r-cran-remotes + - r-cran-terra + - r-cran-igraph + - r-cran-sf + - r-cran-matrix + - r-cran-mass + - r-cran-e1071 + - r-cran-dbscan + - r-cran-data.table + - r-cran-dplyr + - r-cran-ggplot2 + - r-cran-reshape2 + - r-cran-stringr + - r-cran-fs + - r-cran-lmtest + - r-cran-matrixstats + - r-cran-geometry + - r-cran-spatstat.geom + - type: docker + # Install FastReseg from GitHub. upgrade = "never" so remotes reuses the + # apt-installed binary deps above instead of re-fetching/recompiling them + # from CRAN; only the unpackaged deps (concaveman, GiottoUtils, + # GiottoClass) and FastReseg itself are compiled from source. + run: + - Rscript -e 'options(warn = 2); remotes::install_github("Nanostring-Biostats/FastReseg", upgrade = "never", repos = "https://cran.rstudio.com")' + - type: python + # `plankton` (pulled in by txsim.run_ssam for ssam cell annotation in + # input.py) still does `from matplotlib.cm import get_cmap`, removed in + # matplotlib 3.9. Pin < 3.9 so run_ssam imports. This is a local setup + # step, so it runs AFTER the __merge__'d spatialdata/txsim partials (which + # otherwise pull matplotlib >= 3.9) and wins. Mirrors methods_cell_type_annotation/ssam. + pypi: [planktonspace, "matplotlib<3.9"] + __merge__: - /src/base/setup_txsim_partial.yaml - /src/base/setup_spatialdata_partial.yaml - type: native diff --git a/src/methods_transcript_assignment/fastreseg/input.py b/src/methods_transcript_assignment/fastreseg/input.py index 1da6f5a48..9eabfe0ab 100644 --- a/src/methods_transcript_assignment/fastreseg/input.py +++ b/src/methods_transcript_assignment/fastreseg/input.py @@ -7,6 +7,61 @@ import txsim as tx import anndata as ad import argparse +from scipy.sparse import csr_matrix + + +def generate_adata(input_spots, cell_id_col="cell", gene_col="Gene"): + """Aggregate per-transcript assignments into a cell x gene AnnData. + + Reimplementation of txsim.preprocessing.generate_adata: txsim's version is + incompatible with anndata>=0.12 (it passes the removed `dtype=` argument to + AnnData() and uses per-row item assignment `adata[cell, :] = ...`, which + anndata no longer supports). This fills the count matrix directly and + produces an identical output structure (X, layers['raw_counts'], + obs/var stats, uns['spots'/'pct_noise']). Kept in sync with the copy in + src/methods_count_aggregation/basic_count_aggregation/script.py. + """ + spots = input_spots.copy() + pct_noise = sum(spots[cell_id_col] <= 0) / len(spots[cell_id_col]) + spots_raw = spots.copy() # kept in uns['spots']; 0 (background) -> None + spots_raw.loc[spots_raw[cell_id_col] == 0, cell_id_col] = None + spots = spots[spots[cell_id_col] > 0] + + cell_ids = pd.unique(spots[cell_id_col]) + genes = pd.unique(spots[gene_col]) + cell_pos = {c: i for i, c in enumerate(cell_ids)} + + # Populate the count matrix + centroids per cell (no AnnData item assignment). + # feature_name may be categorical, so value_counts() can return unobserved + # categories; reindex to `genes` (fill 0) to align to the var order, mirroring + # txsim's `.reindex(var_names, fill_value=0)`. + X = np.zeros((len(cell_ids), len(genes)), dtype=np.float32) + centroid_x = np.zeros(len(cell_ids)) + centroid_y = np.zeros(len(cell_ids)) + for cell_id, grp in spots.groupby(cell_id_col, sort=False): + row = cell_pos[cell_id] + cts = grp[gene_col].value_counts().reindex(genes, fill_value=0) + X[row, :] = cts.values + centroid_x[row] = grp["x"].mean() + centroid_y[row] = grp["y"].mean() + + adata = ad.AnnData(csr_matrix(X)) + adata.obs["cell_id"] = cell_ids + adata.obs_names = [f"{i:d}" for i in cell_ids] + adata.var_names = genes + adata.obs["centroid_x"] = centroid_x + adata.obs["centroid_y"] = centroid_y + + adata.uns["spots"] = spots_raw + adata.uns["pct_noise"] = pct_noise + adata.layers["raw_counts"] = adata.X.copy() + + adata.obs["n_counts"] = np.ravel(adata.layers["raw_counts"].sum(axis=1)) + adata.obs["n_genes"] = adata.layers["raw_counts"].getnnz(axis=1) + adata.var["n_counts"] = np.ravel(adata.layers["raw_counts"].sum(axis=0)) + adata.var["n_cells"] = adata.layers["raw_counts"].getnnz(axis=0) + return adata + def parse_arguments(): parser = argparse.ArgumentParser( @@ -70,8 +125,12 @@ def parse_arguments(): print('Transforming transcripts coordinates', flush=True) transcripts = sd.transform(sdata['transcripts'], to_coordinate_system='global') -transcripts_df = transcripts.compute() -transcripts_df.rename(columns = {'feature_name': 'target', +# .copy() is required: for a single-partition points object, transcripts.compute() +# returns a reference to the dask graph's backing frame, so an in-place rename here +# would corrupt it and the later transcripts.compute() (passed to run_ssam) would +# yield renamed columns that no longer match the dask _meta -> "Metadata mismatch". +transcripts_df = transcripts.compute().copy() +transcripts_df.rename(columns = {'feature_name': 'target', 'transcript_id': 'UMI_transID', 'cell_id': 'UMI_cellID'}, inplace = True) transcripts_df = transcripts_df.loc[:, ['target', 'x', 'y', 'z', 'UMI_transID', 'UMI_cellID']] @@ -82,7 +141,7 @@ def parse_arguments(): df = sdata['transcripts'].compute() df.feature_name = df.feature_name.astype(str) -adata_sp = tx.preprocessing.generate_adata(df, cell_id_col='cell_id', gene_col='feature_name') #TODO: x and y refers to a specific coordinate system. Decide which space we want to use here. (probably should be handled in the previous assignment step) +adata_sp = generate_adata(df, cell_id_col='cell_id', gene_col='feature_name') #TODO: x and y refers to a specific coordinate system. Decide which space we want to use here. (probably should be handled in the previous assignment step) adata_sp.layers['counts'] = adata_sp.layers['raw_counts'] del adata_sp.layers['raw_counts'] adata_sp.var["gene_name"] = adata_sp.var_names diff --git a/src/methods_transcript_assignment/fastreseg/orchestrator.sh b/src/methods_transcript_assignment/fastreseg/orchestrator.sh index 9529d5b51..34d1eff74 100644 --- a/src/methods_transcript_assignment/fastreseg/orchestrator.sh +++ b/src/methods_transcript_assignment/fastreseg/orchestrator.sh @@ -1,5 +1,10 @@ #!/bin/bash +# Fail loudly: without this, a crash in any sub-step (input.py / script.R / +# output.py) is swallowed and Viash only reports the generic "required output +# file is missing", hiding the real Python/R traceback. +set -eo pipefail + ## VIASH START # The following code has been auto-generated by Viash. par_input_ist='resources_test/task_ist_preprocessing/mouse_brain_combined/raw_ist.zarr' diff --git a/src/methods_transcript_assignment/fastreseg/output.py b/src/methods_transcript_assignment/fastreseg/output.py index 7849d351e..dbca407fb 100644 --- a/src/methods_transcript_assignment/fastreseg/output.py +++ b/src/methods_transcript_assignment/fastreseg/output.py @@ -36,6 +36,12 @@ def convert_to_lower_dtype(arr): ##converting to dask transcript output for spatialdata format df = df.loc[:,["x", "y", "z", "target", "updated_cellID", "updated_celltype", "UMI_transID"]] df.rename(columns={"updated_cellID": "cell_id", "target": "feature_name", "UMI_transID": "transcript_id"}, inplace=True) +# Match the raw_ist convention where feature_name is categorical. Rebuilt from +# R's CSV it would otherwise be a pyarrow-backed nullable string that survives +# the zarr round-trip; the downstream count-aggregation then builds var_names +# from it and write_h5ad refuses to serialize the nullable-string index +# ("anndata.settings.allow_write_nullable_strings is None"). +df["feature_name"] = df["feature_name"].astype(str).astype("category") transcripts_dask = dask.dataframe.from_pandas(df, npartitions = 1) diff --git a/src/methods_transcript_assignment/segger/script.py b/src/methods_transcript_assignment/segger/script.py index 4293760d5..610fc203c 100644 --- a/src/methods_transcript_assignment/segger/script.py +++ b/src/methods_transcript_assignment/segger/script.py @@ -34,6 +34,23 @@ } ## VIASH END +# The Nextflow/k8s GPU task hands us an EMPTY CUDA_VISIBLE_DEVICES. An empty (or +# whitespace) value masks every GPU from the CUDA *driver* API, so cudf/cuspatial and +# numba enumerate zero devices and segger dies deep in tiling with +# "cudaErrorNoDevice: no CUDA-capable device is detected" / numba's +# "IndexError: list index out of range" (self.gpus[devnum] on an empty device list). +# torch's is_available() is NVML-based and ignores CVD, so the gate below still passes +# and the failure only surfaces inside the `segger segment` subprocess. Drop an empty +# value so the allocated device (exposed via NVIDIA_VISIBLE_DEVICES) is visible again; +# run_segger's os.environ.copy() then inherits the fix for the subprocess. +if os.environ.get("CUDA_VISIBLE_DEVICES", "x").strip() == "": + print( + "CUDA_VISIBLE_DEVICES was empty; unsetting so the allocated GPU is visible " + "to cudf/numba (torch's NVML check hides this).", + flush=True, + ) + del os.environ["CUDA_VISIBLE_DEVICES"] + # segger runs its tiling / GNN training / prediction end-to-end on the GPU # (cudf, cuspatial, torch kernels). Fail fast with a clear message rather than # deep inside the subprocess when no CUDA device is present.