diff --git a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh index b695b0f08..4d4a3aebe 100644 --- a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh +++ b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh @@ -19,7 +19,7 @@ # SCRATCH_DIR=/path/to/big/scratch ./mirror_bruker_to_s3.sh set -euo pipefail - +SCRATCH_DIR="/Volumes/SeagateHHD" # --- Config ----------------------------------------------------------------- SRC_BASE="https://smi-public.objects.liquidweb.services" @@ -33,11 +33,11 @@ SCRATCH_DIR="${SCRATCH_DIR:-$PWD/bruker_mirror_scratch}" # - SOURCE is either a full https URL, or a name relative to SRC_BASE (the mouse/liver host). # The URL-encoded names are what the liquidweb server serves. # - LOCAL_NAME is the (decoded) name to store under on S3. -FILES=( - "HalfBrain.zip|HalfBrain.zip" - "Half%20%20Brain%20simple%20%20files%20.zip|Half Brain simple files.zip" - "NormalLiverFiles.zip|NormalLiverFiles.zip" -) +#FILES=( +# "HalfBrain.zip|HalfBrain.zip" +# "Half%20%20Brain%20simple%20%20files%20.zip|Half Brain simple files.zip" +# "NormalLiverFiles.zip|NormalLiverFiles.zip" +#) # NSCLC lung-cancer samples: each ships a flat-files+cell-labels archive and a # raw-morphology-images archive. The bruker_cosmx_nsclc loader streams both from S3. @@ -132,7 +132,7 @@ for entry in "${FILES[@]}"; do # Upload to S3 echo " Uploading to s3://$bucket/$key ..." - aws s3 cp "$local_path" "s3://$bucket/$key" + aws s3 cp "$local_path" "s3://$bucket/$key" --profile op # Free scratch space before the next (much larger) file echo " Removing local copy to free space" diff --git a/src/base/setup_spatialdata_partial.yaml b/src/base/setup_spatialdata_partial.yaml index 7e3f4be48..c887a67f0 100644 --- a/src/base/setup_spatialdata_partial.yaml +++ b/src/base/setup_spatialdata_partial.yaml @@ -1,3 +1,3 @@ setup: - type: python - pypi: ["spatialdata>=0.7.3", "anndata>=0.12.0", "zarr>=3.0.0"] + pypi: ["spatialdata>=0.7.3", "anndata>=0.12.0,<0.13", "zarr>=3.0.0"] diff --git a/src/methods_expression_correction/resolvi_correction/config.vsh.yaml b/src/methods_expression_correction/resolvi_correction/config.vsh.yaml index d494ebcc0..543c379aa 100644 --- a/src/methods_expression_correction/resolvi_correction/config.vsh.yaml +++ b/src/methods_expression_correction/resolvi_correction/config.vsh.yaml @@ -48,7 +48,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0", scvi-tools] + pypi: ["anndata>=0.12.0,<0.13", scvi-tools] - type: native runners: diff --git a/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml b/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml index 6f9ea9013..5be364d9c 100644 --- a/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml +++ b/src/methods_gene_efficiency_correction/gene_efficiency_correction/config.vsh.yaml @@ -31,7 +31,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0"] + pypi: ["anndata>=0.12.0,<0.13"] - type: native runners: diff --git a/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml b/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml index fca3c2bd6..624672965 100644 --- a/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml +++ b/src/methods_qc_filter/basic_qc_filter/config.vsh.yaml @@ -33,7 +33,7 @@ engines: - /src/base/setup_txsim_partial.yaml setup: - type: python - pypi: ["anndata>=0.12.0"] + pypi: ["anndata>=0.12.0,<0.13"] - type: native runners: diff --git a/src/methods_transcript_assignment/fastreseg/input.py b/src/methods_transcript_assignment/fastreseg/input.py index 8279a4d65..967d9ed41 100644 --- a/src/methods_transcript_assignment/fastreseg/input.py +++ b/src/methods_transcript_assignment/fastreseg/input.py @@ -132,6 +132,20 @@ def parse_arguments(): label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() else: label_image = sdata_segm["segmentation"].to_numpy() +# Guard against an empty segmentation. If the label image has no cells (all zeros) +# every transcript is assigned to background (cell_id 0), generate_adata drops them +# all, and the downstream tacco annotation then divides by a zero count norm and +# fails with a cryptic "divide by zero" instead of a useful message. Observed for the +# vizgen lung-cancer merscope dataset, whose polygon->label rasterization produced an +# all-zero image (for custom_segmentation this means the dataset's `cell_labels` is +# empty; re-process the dataset). Mirrors the guard in methods_transcript_assignment/pciseq. +if label_image.max() == 0: + raise ValueError( + f"Segmentation '{input_segmentation_path}' contains no cells: the label image " + f"(shape={label_image.shape}, dtype={label_image.dtype}) is entirely zero. FastReseg " + f"cannot assign transcripts to cells " + f"(for custom_segmentation this means the dataset's `cell_labels` is empty)." + ) # Clamp coords to the label-image bounds: transcripts at the crop boundary can # round a few pixels past the raster edge, or the transcript field of view can # extend beyond the segmented raster (real data), so label_image[y, x] would diff --git a/src/methods_transcript_assignment/segger/script.py b/src/methods_transcript_assignment/segger/script.py index 040417644..6bbb28ad6 100644 --- a/src/methods_transcript_assignment/segger/script.py +++ b/src/methods_transcript_assignment/segger/script.py @@ -76,6 +76,25 @@ segmentation_coord_systems = sd.transformations.get_transformation(sdata_segm["segmentation"], get_all=True).keys() assert par['coordinate_system'] in segmentation_coord_systems, f"Coordinate system '{par['coordinate_system']}' not found in input data." +# Guard against an empty (all-background) segmentation. sd.to_polygons() on an all-zero +# label image returns an EMPTY GeoDataFrame, which then crashes deep inside spatialdata's +# ShapesModel.parse (`data["geometry"].iloc[0]` -> IndexError: "single positional indexer +# is out-of-bounds"). This happens for datasets whose cell_labels are all-zero from the +# stale pre-0.8.0 spatialdata rasterize bug -- e.g. 2022_vizgen_human_lung_cancer_merfish_ +# combined via custom_segmentation (which just copies cell_labels). The real fix is to +# re-process that dataset under spatialdata>=0.8.0; fail fast here with an actionable +# message rather than the cryptic pandas IndexError. Mirrors the guard in pciseq/script.py. +_seg_elem = sdata_segm["segmentation"] +_seg_arr = _seg_elem["scale0"].image if isinstance(_seg_elem, xr.DataTree) else _seg_elem +if int(_seg_arr.data.max()) == 0: + raise ValueError( + "Segmentation label image is empty (all background, no cells) -- segger cannot " + "assign transcripts to zero cells. This typically means the dataset's cell_labels " + "are all-zero from the stale pre-0.8.0 spatialdata rasterize bug (e.g. " + "2022_vizgen_human_lung_cancer_merfish_combined via custom_segmentation). " + "Re-process that dataset under spatialdata>=0.8.0 to fix it at the source." + ) + ########################################## # boundaries.parquet (from segmentation) # ########################################## @@ -141,15 +160,35 @@ label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() else: label_image = sdata_segm["segmentation"].to_numpy() -# Clip to the label image bounds (transcripts can sit just outside after transform). -n_oob = int(np.count_nonzero((y_coords < 0) | (y_coords >= label_image.shape[0]) | (x_coords < 0) | (x_coords >= label_image.shape[1]))) -y_coords = np.clip(y_coords, 0, label_image.shape[0] - 1) -x_coords = np.clip(x_coords, 0, label_image.shape[1] - 1) -print(f"Clamped {n_oob}/{len(x_coords)} transcripts outside the {label_image.shape[0]}x{label_image.shape[1]} label image to its edge", flush=True) -prior_cell_id = label_image[y_coords, x_coords].astype(np.int64) # 0 == background - -# Canonical transcripts frame (native coordinates, clean RangeIndex). Its row -# order matches prior_cell_id and the row_index segger reports back below. +# Transcripts can land outside the label image after the transform (e.g. the combine step +# crops the image/labels to ~20000x20000 but transcripts span a larger field -> here ~38% +# of them fall outside). Those OOB transcripts must NOT be fed to segger: segger tiles the +# transcript field into a heterogeneous graph with 'tx' (transcript) and 'bd' (boundary/ +# cell) node types, and a tile covering a transcripts-only region (no boundaries) produces a +# training mini-batch with ZERO 'bd' nodes. segger's positional encoder then runs +# `torch.zeros((batch.max()+1, 2))` on an empty batch -> `RuntimeError: max(): Expected +# reduction dim to be specified for input.numel() == 0`, which kills training MID-RUN (only +# when the DataLoader's shuffle happens to draw an all-OOB mini-batch, so it can survive +# several epochs first). A transcript outside the segmentation can't belong to any segmented +# cell anyway, so we EXCLUDE it from the segger input and leave it unassigned (cell_id 0) in +# the final output. Datasets whose segmentation covers the whole transcript field are +# unaffected (in_bounds all-True, n_oob == 0). +in_bounds = ( + (y_coords >= 0) & (y_coords < label_image.shape[0]) + & (x_coords >= 0) & (x_coords < label_image.shape[1]) +) +n_oob = int((~in_bounds).sum()) +y_look = np.clip(y_coords, 0, label_image.shape[0] - 1) +x_look = np.clip(x_coords, 0, label_image.shape[1] - 1) +prior_cell_id = label_image[y_look, x_look].astype(np.int64) # 0 == background +prior_cell_id[~in_bounds] = 0 # OOB -> background (excluded from segger regardless) +print(f"{n_oob}/{len(in_bounds)} transcripts fall outside the " + f"{label_image.shape[0]}x{label_image.shape[1]} label image; excluding them from the " + f"segger input (they stay unassigned in the output).", flush=True) + +# Canonical transcripts frame (native coordinates, clean RangeIndex). Its row order matches +# prior_cell_id / in_bounds; the IN-BOUNDS subset's order matches the parquet we write and +# hence segger's reported row_index (mapped back through seg_orig_idx below). tx_pd = transcripts_reset.compute() n_tx = len(tx_pd) @@ -158,19 +197,30 @@ else: transcript_id = np.arange(n_tx, dtype=np.uint64) +# Full-frame positions of the transcripts actually handed to segger. segger's output +# row_index indexes into the (in-bounds) parquet row order; seg_orig_idx maps it back. +seg_orig_idx = np.nonzero(in_bounds)[0] +if seg_orig_idx.size == 0: + raise ValueError( + "No transcripts fall within the segmentation label image, so segger has nothing to " + "assign. This usually means the transcripts and segmentation are in mismatched " + "coordinate frames (check the dataset's crop / transforms)." + ) + +n_in = int(in_bounds.sum()) tx_out = pd.DataFrame({ - "transcript_id": transcript_id, - "x_location": tx_pd["x"].to_numpy().astype(np.float32), - "y_location": tx_pd["y"].to_numpy().astype(np.float32), - "feature_name": tx_pd["feature_name"].astype(str).to_numpy(), + "transcript_id": transcript_id[in_bounds], + "x_location": tx_pd["x"].to_numpy().astype(np.float32)[in_bounds], + "y_location": tx_pd["y"].to_numpy().astype(np.float32)[in_bounds], + "feature_name": tx_pd["feature_name"].astype(str).to_numpy()[in_bounds], # segger's Xenium loader keys transcripts to boundaries by cell_id string; # background (0) becomes the UNASSIGNED sentinel. - "cell_id": np.where(prior_cell_id > 0, prior_cell_id.astype(str), "UNASSIGNED"), - "qv": np.full(n_tx, 40.0, dtype=np.float32), - "overlaps_nucleus": (prior_cell_id > 0).astype(np.int8), + "cell_id": np.where(prior_cell_id[in_bounds] > 0, prior_cell_id[in_bounds].astype(str), "UNASSIGNED"), + "qv": np.full(n_in, 40.0, dtype=np.float32), + "overlaps_nucleus": (prior_cell_id[in_bounds] > 0).astype(np.int8), }) if "z" in tx_pd.columns: - tx_out.insert(3, "z_location", tx_pd["z"].to_numpy().astype(np.float32)) + tx_out.insert(3, "z_location", tx_pd["z"].to_numpy().astype(np.float32)[in_bounds]) tx_out.to_parquet(XENIUM_DIR / "transcripts.parquet", index=False) del transcripts, y_coords, x_coords, label_image @@ -291,14 +341,15 @@ def run_segger(node_dim): seg = seg[keep_mask] print(f"segger kept {len(seg)} assignments of {n_tx} transcripts", flush=True) -# segger reports row_index into the transcripts.parquet we wrote (i.e. tx_pd -# row order). Factorize the (string) segger_cell_id into contiguous positive -# integers; unassigned transcripts stay 0. +# segger reports row_index into the transcripts.parquet we wrote, which is the IN-BOUNDS +# subset (OOB transcripts were excluded above). Map that back to the full-frame position via +# seg_orig_idx, then factorize the (string) segger_cell_id into contiguous positive integers; +# unassigned / out-of-bounds transcripts stay 0. cell_id_per_tx = np.zeros(n_tx, dtype=np.int64) if len(seg): row_idx = seg["row_index"].to_numpy().astype(np.int64) codes, _ = pd.factorize(seg["segger_cell_id"].astype(str), sort=True) - cell_id_per_tx[row_idx] = codes.astype(np.int64) + 1 + cell_id_per_tx[seg_orig_idx[row_idx]] = codes.astype(np.int64) + 1 ############################################# # Build transcript-assignment output object #