Fixes

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"

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"

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:

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)
 ]

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]

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"])