Fixes

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"

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

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:

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:

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:

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

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 #