diff --git a/scripts/create_resources/combine/process_datasets_bruker_nebius.sh b/scripts/create_resources/combine/process_datasets_bruker_nebius.sh index dbc103819..a55151794 100644 --- a/scripts/create_resources/combine/process_datasets_bruker_nebius.sh +++ b/scripts/create_resources/combine/process_datasets_bruker_nebius.sh @@ -12,7 +12,11 @@ cd "$REPO_ROOT" set -e -input_dir="s3://openproblems-data/resources/datasets" +# Read raw inputs from the scratch PVC (where the bruker spatial scripts publish via +# publish_dir='/scratch/task_ist_preprocessing/raw'), not S3. Both input_sp (bruker_cosmx/…) +# and input_sc (scrnaseq refs) resolve under here; the raw/ layout mirrors the S3 subpaths. +input_dir="/scratch/task_ist_preprocessing/raw" +#input_dir="s3://openproblems-data/resources/datasets" #publish_dir="s3://openproblems-data/resources/task_ist_preprocessing/datasets" publish_dir='/scratch/task_ist_preprocessing/datasets' @@ -41,17 +45,6 @@ param_list: dataset_description: "Bruker CosMx Human Liver + 2022 Andrews scRNAseq" dataset_organism: "homo_sapiens" - - id: "bruker_human_liver_cancer_cosmx_combined" - input_sp: "$input_dir/bruker_cosmx/bruker_human_liver_cancer_cosmx/dataset.zarr" - input_sc: "$input_dir/scrnaseq_for_ist/2022Lu_human_liver_cancer_sc/dataset.h5ad" - dataset_id: "bruker_human_liver_cancer_cosmx_combined" - dataset_name: "Human liver cancer combined Bruker CosMx 2022 Lu scRNAseq" - dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/" - dataset_reference: "https://doi.org/10.1038/s41467-022-32283-3" - dataset_summary: "Bruker CosMx Human Liver Cancer + 2022 Lu scRNAseq" - dataset_description: "Bruker CosMx Human Liver Cancer + 2022 Lu scRNAseq" - dataset_organism: "homo_sapiens" - - id: "bruker_human_lung_cancer_cosmx_combined/lung5_rep1" input_sp: "$input_dir/bruker_cosmx/bruker_human_lung_cancer_cosmx/lung5_rep1/dataset.zarr" input_sc: "$input_dir/zuani_human_nsclc_sc/2024Zuani_human_nsclc_sc/dataset.h5ad" diff --git a/src/base/setup_spatialdata_partial.yaml b/src/base/setup_spatialdata_partial.yaml index 9c34899a3..9155a7ec8 100644 --- a/src/base/setup_spatialdata_partial.yaml +++ b/src/base/setup_spatialdata_partial.yaml @@ -1,6 +1,13 @@ setup: - type: python - # Pinned exactly: an unpinned >=0.7.3 silently pulled 0.8.0 on a rebuild, whose - # bounding_box/filter_table regressions transposed/emptied transcripts in the crop. - # process_dataset's transform-agnostic crop is validated against 0.8.0, so pin here. - pypi: ["spatialdata==0.8.0", "anndata>=0.12.0,<0.13", "zarr>=3.0.0"] + # Pin exactly per python version: spatialdata 0.8.0 requires python>=3.12, but this + # partial is shared by components on older python (e.g. clustermap) that can only take + # 0.7.3. Pinning stops the unpinned >=0.7.3 from silently drifting again (that 0.7.3->0.8.0 + # drift is what transposed/emptied transcripts in process_dataset's crop). py>=3.12 + # components (process_dataset, vizgen loader) get 0.8.0, where the transform-agnostic crop + # is validated; older-python components stay on 0.7.3 (what they already resolved to). + pypi: + - "spatialdata==0.8.0; python_version>='3.12'" + - "spatialdata>=0.7.3; python_version<'3.12'" + - "anndata>=0.12.0,<0.13" + - "zarr>=3.0.0" diff --git a/src/datasets/loaders/allen_brain_cell_atlas/config.vsh.yaml b/src/datasets/loaders/allen_brain_cell_atlas/config.vsh.yaml index e3c57636c..d5d71c2e3 100644 --- a/src/datasets/loaders/allen_brain_cell_atlas/config.vsh.yaml +++ b/src/datasets/loaders/allen_brain_cell_atlas/config.vsh.yaml @@ -100,8 +100,11 @@ engines: image: openproblems/base_python:1 setup: - type: python - packages: - - git+https://github.com/alleninstitute/abc_atlas_access + packages: + # Pin the tag: an unpinned git install drifted to an API where AbcProjectCache + # dropped list_data_files (split into list_expression_matrix_files etc.), breaking + # the loader. v.3.3.0 has list_expression_matrix_files + get_data_path (deprecated). + - git+https://github.com/alleninstitute/abc_atlas_access@v.3.3.0 - type: native runners: diff --git a/src/datasets/loaders/allen_brain_cell_atlas/script.py b/src/datasets/loaders/allen_brain_cell_atlas/script.py index 8a01bcf42..120e2794f 100644 --- a/src/datasets/loaders/allen_brain_cell_atlas/script.py +++ b/src/datasets/loaders/allen_brain_cell_atlas/script.py @@ -73,14 +73,14 @@ obs = obs.sample(n=par["sample_n_obs"]) -# From abc_cache.list_data_files("WMB-10Xv2") +# From abc_cache.list_expression_matrix_files("WMB-10Xv2") # TODO: potentially also load other chemistries (currently only 10Xv2) print("Downloading and reading expression matrices", flush=True) abca_data_subdir="WMB-10Xv2" abca_region_files = [ [file_name, region] - for file_name in abc_cache.list_data_files(abca_data_subdir) + for file_name in abc_cache.list_expression_matrix_files(abca_data_subdir) for region in REGIONS if re.match(f"WMB-10Xv2-{region}[\\-0-9]*/raw", file_name) ] diff --git a/src/datasets/processors/subsample/config.vsh.yaml b/src/datasets/processors/subsample/config.vsh.yaml new file mode 100644 index 000000000..f9e7fac68 --- /dev/null +++ b/src/datasets/processors/subsample/config.vsh.yaml @@ -0,0 +1,90 @@ +name: subsample +namespace: datasets/processors +description: | + Subsample an h5ad file to create a smaller (test) dataset. + + Local copy of openproblems-bio/datasets `processors/subsample` carrying a + single-modality fix: the mod2 output write is guarded on `--input_mod2` + (whether a second modality was actually read) instead of `--output_mod2`. + Under the Nextflow runner the optional `--output_mod2` path is auto-populated, + so the upstream guard tried to write an `adata_output_mod2` that is only + defined for multimodal input, crashing single-modality runs (Wu, Zuani) with + `NameError: name 'adata_output_mod2' is not defined`. + + Retire this component and switch the run scripts back to the `datasets` + dependency once the upstream fix propagates to openproblems-bio/datasets@build/main. + +arguments: + - name: --input + type: file + required: true + direction: input + description: Input h5ad dataset (must contain a `counts` layer). + - name: --input_mod2 + type: file + required: false + direction: input + description: Optional second-modality h5ad input. + - name: --output + type: file + required: true + direction: output + description: Output subsampled h5ad. + example: output.h5ad + - name: --output_mod2 + type: file + required: false + direction: output + description: Optional second-modality h5ad output (only written when --input_mod2 is given). + example: output_mod2.h5ad + - name: --output_compression + type: string + choices: [gzip, lzf] + required: false + default: gzip + - name: --n_obs + type: integer + default: 500 + description: Maximum number of observations to keep. May end up fewer because empty cells/genes are removed. + - name: --n_vars + type: integer + default: 500 + description: Maximum number of variables to keep. May end up fewer because empty cells/genes are removed. + - name: --keep_features + type: string + multiple: true + required: false + description: A list of genes to keep. + - name: --keep_cell_type_categories + type: string + multiple: true + required: false + description: Cell type categories to be selected. + - name: --keep_batch_categories + type: string + multiple: true + required: false + description: Batch categories to be selected. + - name: --even + type: boolean_true + description: Subsample evenly from different batches. + - name: --seed + type: integer + required: false + description: A seed for the subsampling. + example: 123 + +resources: + - type: python_script + path: script.py + +engines: + - type: docker + image: openproblems/base_python:1 + - type: native + +runners: + - type: executable + - type: nextflow + directives: + label: [midtime, highmem, midcpu] diff --git a/src/datasets/processors/subsample/script.py b/src/datasets/processors/subsample/script.py new file mode 100644 index 000000000..d188550f3 --- /dev/null +++ b/src/datasets/processors/subsample/script.py @@ -0,0 +1,146 @@ +import scanpy as sc +import random +import numpy as np + +### VIASH START +par = { + "input": "resources_test/common/scicar_cell_lines/temp_mod1_full.h5ad", + "input_mod2": "resources_test/common/scicar_cell_lines/temp_mod2_full.h5ad", + "n_obs": 600, + "n_vars": 1500, + "keep_cell_type_categories": None, + "keep_batch_categories": None, + "keep_features": None, + "even": False, + "output": "subsample_mod1.h5ad", + "output_mod2": "subsample_mod2.h5ad", + "output_compression": "gzip", + "seed": 123 +} +### VIASH END + +if par["seed"]: + print(f">> Setting seed to {par['seed']}", flush=True) + random.seed(par["seed"]) + +print(">> Load data", flush=True) +adata_input = sc.read_h5ad(par["input"]) + +if par["input_mod2"] is not None: + adata_mod2 = sc.read_h5ad(par["input_mod2"]) + +# copy counts to .X because otherwise filter_genes and filter_cells won't work +adata_input.X = adata_input.layers["counts"] +if par["input_mod2"] is not None: + adata_mod2.X = adata_mod2.layers["counts"] + +print(">> Determining output shape", flush=True) +min_obs_list = [par["n_obs"], adata_input.shape[0]] +if par["input_mod2"] is not None: + min_obs_list.append(adata_mod2.shape[0]) +n_obs = min(min_obs_list) + +min_vars_list = [par["n_vars"], adata_input.shape[1]] +if par["input_mod2"] is not None: + min_vars_list.append(adata_mod2.shape[1]) +n_vars = min(min_vars_list) + +print(">> Subsampling the observations", flush=True) +obs_filt = np.ones(dtype=np.bool_, shape=adata_input.n_obs) + +# subset by cell_type +if par.get("keep_cell_type_categories"): + print(f">> Selecting cell_type_categories {par['keep_cell_type_categories']}") + obs_filt = obs_filt & adata_input.obs["cell_type"].isin(par["keep_cell_type_categories"]) + +# subset by batch +if par.get("keep_batch_categories"): + print(f">> Selecting cell_type_categories {par['keep_batch_categories']}") + obs_filt = obs_filt & adata_input.obs["batch"].isin(par["keep_batch_categories"]) + +# subsample evenly across batches or not +if par.get("even"): + obs_evenly = "batch" + choice_ix = np.where(obs_filt)[0] + choice_batch = adata_input[choice_ix].obs[obs_evenly] + names, counts = np.unique(choice_batch, return_counts=True) + probs = dict(zip(names, 1 / counts / len(names))) + + choice_probs = [ probs[batch] for batch in choice_batch ] + obs_index = np.random.choice(choice_ix, size=n_obs, replace=False, p=choice_probs) +else: + obs_index = np.random.choice(np.where(obs_filt)[0], n_obs, replace=False) + +# subsample obs +adata_output = adata_input[obs_index].copy() +if par["input_mod2"] is not None: + adata_output_mod2 = adata_mod2[obs_index].copy() + +# filter cells and genes +if par["input_mod2"] is not None: + n_cells = adata_output.X.sum(axis=1).A.flatten() + n_cells_mod2 = adata_output_mod2.X.sum(axis=1).A.flatten() + keep_cells = np.minimum(n_cells, n_cells_mod2) > 1 + adata_output = adata_output[keep_cells, :].copy() + adata_output_mod2 = adata_output_mod2[keep_cells, :].copy() + + sc.pp.filter_genes(adata_output, min_cells=1) + sc.pp.filter_genes(adata_output_mod2, min_cells=1) + +else: + # todo: this should not remove features in keep_features! + print(">> Remove empty observations and features", flush=True) + sc.pp.filter_genes(adata_output, min_cells=1) + sc.pp.filter_cells(adata_output, min_counts=2) + +print(">> Subsampling the features", flush=True) +if par.get("keep_features"): + initial_filt = adata_output.var_names.isin(par["keep_features"]) + initial_idx, *_ = initial_filt.nonzero() + remaining_idx, *_ = (~initial_filt).nonzero() + rest_idx = remaining_idx[np.random.choice(len(remaining_idx), n_vars - len(initial_idx), replace=False)] + var_ix = np.concatenate([initial_idx, rest_idx]) +else: + var_ix = np.random.choice(adata_output.shape[1], n_vars, replace=False) + if par["input_mod2"] is not None: + var_ix_mod2 = np.random.choice(adata_output_mod2.shape[1], n_vars, replace=False) + +# subsample vars +adata_output = adata_output[:, var_ix].copy() +if par["input_mod2"] is not None: + adata_output_mod2 = adata_output_mod2[:, var_ix_mod2].copy() + +# filter cells and genes +if par["input_mod2"] is not None: + n_cells = adata_output.X.sum(axis=1).A.flatten() + n_cells_mod2 = adata_output_mod2.X.sum(axis=1).A.flatten() + keep_cells = np.minimum(n_cells, n_cells_mod2) > 1 + adata_output = adata_output[keep_cells, :].copy() + adata_output_mod2 = adata_output_mod2[keep_cells, :].copy() + + sc.pp.filter_genes(adata_output, min_cells=1) + sc.pp.filter_genes(adata_output_mod2, min_cells=1) + + +else: + # todo: this should not remove features in keep_features! + print(">> Remove empty observations and features", flush=True) + sc.pp.filter_genes(adata_output, min_cells=1) + sc.pp.filter_cells(adata_output, min_counts=2) + +print(">> Update dataset_id", flush=True) +adata_output.uns["dataset_id"] = adata_output.uns["dataset_id"] + "_subsample" +if par["input_mod2"] is not None: + adata_output_mod2.uns["dataset_id"] = adata_output_mod2.uns["dataset_id"] + "_subsample" + +# remove previously copied .X +del adata_output.X +if par["input_mod2"] is not None: + del adata_output_mod2.X + +print(">> Writing data", flush=True) +adata_output.write_h5ad(par["output"], compression=par["output_compression"]) +# Guard on input_mod2 (not output_mod2): adata_output_mod2 only exists for multimodal +# input, but the Nextflow runner auto-populates the optional output_mod2 path. +if par["input_mod2"] is not None: + adata_output_mod2.write_h5ad(par["output_mod2"], compression=par["output_compression"])