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'