diff --git a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh index 4a02bf7db..b695b0f08 100644 --- a/scripts/create_resources/spatial/mirror_bruker_to_s3.sh +++ b/scripts/create_resources/spatial/mirror_bruker_to_s3.sh @@ -23,17 +23,30 @@ set -euo pipefail # --- Config ----------------------------------------------------------------- SRC_BASE="https://smi-public.objects.liquidweb.services" -S3_DEST="${S3_DEST:-s3://openproblems-data/resources_raw/bruker_cosmx}" +# NSCLC lung-cancer samples live on a different NanoString host, one folder per sample. +NSCLC_BASE="https://nanostring-public-share.s3.us-west-2.amazonaws.com/SMI-Compressed" +# Where the loaders read from (see the process_bruker_cosmx*.sh run scripts). +S3_DEST="${S3_DEST:-s3://openproblems-data/resources/raw_data/bruker_cosmx}" SCRATCH_DIR="${SCRATCH_DIR:-$PWD/bruker_mirror_scratch}" -# Files to mirror. The URL-encoded names are what the server serves; the second -# column is the (decoded) name to store under on S3. +# Files to mirror. Each entry is "SOURCE|LOCAL_NAME": +# - 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" ) +# 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. +NSCLC_SAMPLES=(Lung5_Rep1 Lung5_Rep2 Lung5_Rep3 Lung6 Lung9_Rep1 Lung9_Rep2 Lung12 Lung13) +for s in "${NSCLC_SAMPLES[@]}"; do + FILES+=("$NSCLC_BASE/$s/${s}+SMI+Flat+data.tar.gz|${s}+SMI+Flat+data.tar.gz") + FILES+=("$NSCLC_BASE/$s/${s}+RawMorphologyImages.tar.gz|${s}+RawMorphologyImages.tar.gz") +done + # --- Preflight -------------------------------------------------------------- command -v curl >/dev/null || { echo "ERROR: curl not found" >&2; exit 1; } @@ -64,7 +77,11 @@ s3_size() { for entry in "${FILES[@]}"; do url_name="${entry%%|*}" local_name="${entry##*|}" - url="$SRC_BASE/$url_name" + if [[ "$url_name" == http*://* ]]; then + url="$url_name" + else + url="$SRC_BASE/$url_name" + fi local_path="$SCRATCH_DIR/$local_name" # S3 key = everything after s3://bucket/ @@ -126,9 +143,11 @@ done echo "================================================================" echo "All files mirrored to $S3_DEST" echo -echo "Next: update the input_raw / input_flat_files URLs in" -echo " scripts/create_resources/spatial/process_bruker_cosmx_nebius.sh" -echo "to point at the S3 paths, e.g.:" -echo " input_raw: $S3_DEST/HalfBrain.zip" -echo " input_flat_files: $S3_DEST/Half Brain simple files.zip" -echo " input_raw: $S3_DEST/NormalLiverFiles.zip" +echo "Next: the process_bruker_cosmx*.sh run scripts already point at these S3 paths, e.g.:" +echo " mouse/liver (process_bruker_cosmx_nebius.sh):" +echo " input_raw: $S3_DEST/HalfBrain.zip" +echo " input_flat_files: $S3_DEST/Half Brain simple files.zip" +echo " input_raw: $S3_DEST/NormalLiverFiles.zip" +echo " nsclc (process_bruker_cosmx_nsclc_nebius.sh), per sample:" +echo " input_raw: $S3_DEST/+SMI+Flat+data.tar.gz" +echo " input_morphology: $S3_DEST/+RawMorphologyImages.tar.gz" diff --git a/scripts/create_resources/spatial/process_bruker_cosmx_nsclc.sh b/scripts/create_resources/spatial/process_bruker_cosmx_nsclc.sh index 210540b88..bce4bc243 100644 --- a/scripts/create_resources/spatial/process_bruker_cosmx_nsclc.sh +++ b/scripts/create_resources/spatial/process_bruker_cosmx_nsclc.sh @@ -10,11 +10,16 @@ set -e publish_dir="s3://openproblems-data/resources/datasets" +# Raw NSCLC archives mirrored to S3 (see scripts/create_resources/spatial/mirror_bruker_to_s3.sh). +# Each sample ships two tar.gz archives: the flat files + cell labels, and the raw morphology images. +raw_dir="s3://openproblems-data/resources/raw_data/bruker_cosmx" + cat > /tmp/params.yaml << HERE param_list: - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung5_rep1" - sample: "Lung5_Rep1" + input_raw: "$raw_dir/Lung5_Rep1+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung5_Rep1+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung5 Rep1" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung5 Rep1 dataset on FFPE." @@ -23,16 +28,18 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung5_rep2" - sample: "Lung5_Rep2" + input_raw: "$raw_dir/Lung5_Rep2+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung5_Rep2+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung5 Rep2" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung5 Rep2 dataset on FFPE." dataset_description: "Bruker CosMx Human Lung Cancer Lung5 Rep2 dataset on FFPE. Adenocarcinoma, G1, T2aN2M0, IIIA, 75% tumour content." dataset_organism: "homo_sapiens" segmentation_id: ["cell"] - + - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung5_rep3" - sample: "Lung5_Rep3" + input_raw: "$raw_dir/Lung5_Rep3+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung5_Rep3+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung5 Rep3" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung5 Rep3 dataset on FFPE." @@ -41,7 +48,8 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung6" - sample: "Lung6" + input_raw: "$raw_dir/Lung6+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung6+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung6" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung6 dataset on FFPE." @@ -50,7 +58,8 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung9_rep1" - sample: "Lung9_Rep1" + input_raw: "$raw_dir/Lung9_Rep1+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung9_Rep1+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung9 Rep1" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung9 Rep1 dataset on FFPE." @@ -59,7 +68,8 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung9_rep2" - sample: "Lung9_Rep2" + input_raw: "$raw_dir/Lung9_Rep2+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung9_Rep2+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung9 Rep2" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung9 Rep2 dataset on FFPE." @@ -68,7 +78,8 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung12" - sample: "Lung12" + input_raw: "$raw_dir/Lung12+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung12+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung12" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung12 dataset on FFPE." @@ -77,7 +88,8 @@ param_list: segmentation_id: ["cell"] - id: "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung13" - sample: "Lung13" + input_raw: "$raw_dir/Lung13+SMI+Flat+data.tar.gz" + input_morphology: "$raw_dir/Lung13+RawMorphologyImages.tar.gz" dataset_name: "Bruker CosMx Human Lung Cancer Lung13" dataset_url: "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/" dataset_summary: "Bruker CosMx Human Lung Cancer Lung13 dataset on FFPE." @@ -98,4 +110,4 @@ tw launch https://github.com/openproblems-bio/task_ist_preprocessing.git \ --workspace 53907369739130 \ --params-file /tmp/params.yaml \ --config common/nextflow_helpers/labels_tw.config \ - --labels datasets,bruker_cosmx_nsclc \ No newline at end of file + --labels datasets,bruker_cosmx_nsclc diff --git a/src/base/labels_nebius.config b/src/base/labels_nebius.config index 67ab80153..93ab8d178 100644 --- a/src/base/labels_nebius.config +++ b/src/base/labels_nebius.config @@ -173,7 +173,20 @@ process { // cannot write Nextflow's root-owned task scratch dir (.command.log / // .command.begin / .exitcode -> "Permission denied"). Harmless for already-root // GPU images. If segger moves to another gpu label, carry this with it. - pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00jp7hyqr094tmy35'], [imagePullPolicy: 'Always'], [runAsUser: 0]] + // + // emptyDir Memory -> /dev/shm — Kubernetes defaults a container's /dev/shm to 64 MB, + // but segger's Lightning training uses a torch_geometric multi-worker DataLoader that + // moves each collated graph Batch to the main process through /dev/shm (torch's + // file-descriptor sharing strategy, _new_using_fd_cpu). 64 MB overflows on the very + // first batch -> "No space left on device (28)" writing /torch_* / "Unexpected bus + // error ... insufficient shared memory (shm)", killing the DataLoader workers in the + // pre-train sanity-check. NB: the `*sharedmem` labels' `--shm-size` containerOptions do + // NOT help here — containerOptions is a docker/singularity concept the k8s executor + // ignores; on k8s /dev/shm must be enlarged via a pod volume. The tmpfs only consumes + // RAM for pages actually written (idle co-tenants on gpuh100 pay ~0) and counts against + // the pod memory limit, so sizeLimit=16Gi sits comfortably under the 120-180 GB request. + // Carry this with runAsUser:0 if segger moves labels. + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00jp7hyqr094tmy35'], [imagePullPolicy: 'Always'], [runAsUser: 0], [emptyDir: [medium: 'Memory', sizeLimit: '16Gi'], mountPath: '/dev/shm']] containerOptions = { workflow.containerEngine == "singularity" ? '--nv': ( workflow.containerEngine == "docker" ? '--gpus all': null ) } } diff --git a/src/datasets/loaders/bruker_cosmx_nsclc/config.vsh.yaml b/src/datasets/loaders/bruker_cosmx_nsclc/config.vsh.yaml index 71c52fced..9b822cbd2 100644 --- a/src/datasets/loaders/bruker_cosmx_nsclc/config.vsh.yaml +++ b/src/datasets/loaders/bruker_cosmx_nsclc/config.vsh.yaml @@ -5,9 +5,19 @@ argument_groups: - name: Inputs arguments: - type: string - name: --sample - example: "Lung9_Rep1" - description: "Sample id" + name: --input_raw + example: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+SMI+Flat+data.tar.gz" + description: | + URL/path to the flat-files + cell-labels archive (`+SMI+Flat+data.tar.gz`). + Passed as a string (not a staged file) on purpose: the script streams it and extracts + only the members sopa needs. Supports s3:// (via s3fs) or a local path. + - type: string + name: --input_morphology + example: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+RawMorphologyImages.tar.gz" + description: | + URL/path to the raw morphology-images archive (`+RawMorphologyImages.tar.gz`). + Streamed in place; only the 4th z-plane (Z004) of each FOV is extracted. Supports + s3:// (via s3fs) or a local path. - type: string name: --segmentation_id default: ["cell"] @@ -60,9 +70,10 @@ engines: __merge__: - /src/base/setup_spatialdata_partial.yaml setup: - - type: python + - type: python pypi: - sopa + - s3fs # stream the raw tar.gz archives from s3:// and extract only what sopa needs (see --input_raw) - type: native runners: diff --git a/src/datasets/loaders/bruker_cosmx_nsclc/script.py b/src/datasets/loaders/bruker_cosmx_nsclc/script.py index 5085a5b2b..124203b12 100644 --- a/src/datasets/loaders/bruker_cosmx_nsclc/script.py +++ b/src/datasets/loaders/bruker_cosmx_nsclc/script.py @@ -1,57 +1,51 @@ """ -Notes about data structure: - -The directory structure that we need, looks as follows (expected by the sopa.io.cosmx function): -├── CellStatsDir ( = `DATA_DIR`) -│ ├── CellLabels -│ │ ├── CellLabels_F001.tif -│ │ ├── ... -│ │ └── CellLabels_F.tif -│ ├── Morphology2D -│ │ ├── _F001.tif -│ │ ├── ... -│ │ └── _F.tif -│ ├── _exprMat_file.csv -│ ├── _fov_positions_file.csv -│ ├── _metadata_file.csv -│ ├── _tx_file.csv -│ └── -polygons.csv -├── (AnalysisResults) -└── (RunSummary) - -The NSCLC samples download two files: -1. +SMI+Flat+data.tar.gz --> decompressed: - -└── /-Flat_files_and_images ( = `DATA_DIR`) - ├── CellLabels - │ ├── CellLabels_F001.tif - │ ├── ... - │ └── CellLabels_F.tif - ├── _exprMat_file.csv - ├── _fov_positions_file.csv - ├── _metadata_file.csv - ├── _tx_file.csv - └── NOTE: missing: -polygons.csv - -2. +RawMorphologyImages.tar.gz --> decompressed: - -└── 2/-RawMorphologyImages ( = `MORPHOLOGY_DIR`) - ├── _F001_Z001.TIF - ├── _F001_Z002.TIF - ├── ... - ├── _F001_Z008.TIF - ├── ... - └── _F_Z008.TIF - -The morphology images need to be moved to DATA_DIR / "Morphology2D". Also, we only take the 4th plane (Z004) of each image -and rename e.g. _F001_Z004.TIF to _F001.tif. - - +Notes about supported data structure (CosMx NSCLC / lung cancer, "Version 3"): + +The general `bruker_cosmx` loader handles the mouse-brain / liver CosMx releases (one big zip, +optionally a separate flat-files zip). The NSCLC lung-cancer release is shaped differently and +is handled here, but this loader now follows the SAME pattern as the general one: the raw +archives are mirrored to S3 once (see scripts/create_resources/spatial/mirror_bruker_to_s3.sh) +and streamed in place at load time — only the members sopa actually needs are extracted, so the +enormous per-z-plane morphology stacks never touch scratch. The archives are passed as strings +(not staged Viash files) so Nextflow does not copy them first. + +Each NSCLC sample ships TWO gzipped tarballs on the NanoString share: + +1. +SMI+Flat+data.tar.gz (par["input_raw"]) --> decompressed: + + └── /-Flat_files_and_images ( = `DATA_DIR`) + ├── CellLabels + │ ├── CellLabels_F001.tif + │ ├── ... + │ └── CellLabels_F.tif + ├── _exprMat_file.csv + ├── _fov_positions_file.csv + ├── _metadata_file.csv + ├── _tx_file.csv + └── NOTE: missing: -polygons.csv + +2. +RawMorphologyImages.tar.gz (par["input_morphology"]) --> decompressed: + + └── 2/-RawMorphologyImages ( = morphology dir) + ├── _F001_Z001.TIF + ├── ... + ├── _F001_Z008.TIF + ├── ... + └── _F_Z008.TIF + +The morphology images have one TIF per z-plane; sopa expects a single 2D image per FOV under +`DATA_DIR / "Morphology2D"`. We keep only the 4th plane (Z004) — this is the important filter, +it drops ~7/8 of the (very large) morphology archive — and rename e.g. +`_F001_Z004.TIF` to `_F001.tif` under `Morphology2D`. + +The NSCLC flat files carry no `-polygons.csv`, so `cell_boundaries` is derived from the stitched +label image (`sd.to_polygons`), matching the original NSCLC loader behaviour. """ import os +import re import shutil import tarfile from pathlib import Path @@ -61,10 +55,14 @@ ## VIASH START par = { - "sample": "Lung9_Rep1", + # Mirrored from + # https://nanostring-public-share.s3.us-west-2.amazonaws.com/SMI-Compressed/Lung9_Rep1/Lung9_Rep1+SMI+Flat+data.tar.gz + "input_raw": "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+SMI+Flat+data.tar.gz", + # https://nanostring-public-share.s3.us-west-2.amazonaws.com/SMI-Compressed/Lung9_Rep1/Lung9_Rep1+RawMorphologyImages.tar.gz + "input_morphology": "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+RawMorphologyImages.tar.gz", "segmentation_id": ["cell"], "output": "output.zarr", - "dataset_id": "bruker_cosmx/bruker_human_nsclc_cosmx/lung9_rep1", + "dataset_id": "bruker_cosmx/bruker_human_lung_cancer_cosmx/lung9_rep1", "dataset_name": "value", "dataset_url": "https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/nsclc-ffpe-dataset/", "dataset_reference": "value", @@ -73,71 +71,185 @@ "dataset_organism": "human", } meta = { - #"temp_dir": "./temp/datasets/bruker_cosmx", - "temp_dir": "/Volumes/T7/G3_temp/datasets/cosmx/test_folder" + "temp_dir": "./temp/datasets/bruker_cosmx_nsclc/temp", } - ## VIASH END assert ("cell" in par["segmentation_id"]) and (len(par["segmentation_id"]) == 1), "Currently cell labels are definitely assigned in this script. And cosmx does not provide other segmentations." t0 = datetime.now() -# Download urls for sample -URL = f"https://nanostring-public-share.s3.us-west-2.amazonaws.com/SMI-Compressed/{par['sample']}/{par['sample']}+SMI+Flat+data.tar.gz" -URL_MORPHOLOGY = f"https://nanostring-public-share.s3.us-west-2.amazonaws.com/SMI-Compressed/{par['sample']}/{par['sample']}+RawMorphologyImages.tar.gz" +def log(msg): + print(datetime.now() - t0, msg, flush=True) -# Define temp dir and file names +# Define temp dir TMP_DIR = Path(meta["temp_dir"] or "/tmp") TMP_DIR.mkdir(parents=True, exist_ok=True) -FILE_NAME = TMP_DIR / (par["sample"] + "+SMI+Flat+data.tar.gz") -FILE_NAME_MORPHOLOGY = TMP_DIR / (par["sample"] + "+RawMorphologyImages.tar.gz") -DATA_DIR = TMP_DIR / par["sample"] / (par["sample"] + "-Flat_files_and_images") -MORPHOLOGY_DIR = TMP_DIR / par["sample"] / (par["sample"] + "-RawMorphologyImages") - -# Download flat files and cell labels -print(datetime.now() - t0, "Download flat files and cell labels", flush=True) -os.system(f"wget {URL} -O '{FILE_NAME}'") - -# Extract tar.gz files -print(datetime.now() - t0, "Extract tar.gz of flat files and cell labels", flush=True) -with tarfile.open(FILE_NAME, "r:gz") as tar: - tar.extractall(TMP_DIR) - -# Check that flat files are present -FLAT_FILES_ENDINGS = ["_exprMat_file.csv", "_fov_positions_file.csv", "_metadata_file.csv", "_tx_file.csv"] #, "polygons.csv"] + +# Only the 4th morphology z-plane is kept (see module docstring). The regex is +# case-insensitive because the raw files use `.TIF`. +_Z_PLANE_RE = re.compile(r"_Z004\.tif$", re.I) + + +def _open_source(src): + """Seekable binary handle for a local path or an s3:// URL (streamed, never fully downloaded).""" + s = str(src) + if s.startswith("s3://"): + import s3fs + # openproblems-data is public; anonymous read is deterministic and needs no creds. + fs = s3fs.S3FileSystem(anon=True, default_block_size=32 * 2**20) + return fs.open(s, "rb") + return open(s, "rb") + + +def _member_wanted_default(name: str) -> bool: + parts = name.split("/") + return not (name.startswith("__MACOSX/") or parts[-1] == ".DS_Store") + + +def _member_wanted_morphology(name: str) -> bool: + return _member_wanted_default(name) and bool(_Z_PLANE_RE.search(name)) + + +def extract_targz(input_targz, output_dir: Path, wanted=_member_wanted_default): + """ + Stream a .tar.gz (local path or s3:// URL) into output_dir, writing only the members for + which `wanted(name)` is True. The archive is read once, sequentially (gzip is not seekable), + but only wanted members are decoded to disk — so for the morphology tarball only the ~1/8 of + files that are Z004 planes ever hit scratch. The archive itself is never fully downloaded to + a local file first. + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + kept = skipped = 0 + kept_bytes = 0 + with _open_source(input_targz) as fh, tarfile.open(fileobj=fh, mode="r|gz") as tar: + for member in tar: + if not member.isfile() or not wanted(member.name): + skipped += 1 + continue + # Flatten: strip any leading directory components, keep a stable relative path so + # globbing below is layout-independent. + target = output_dir.joinpath(*Path(member.name).parts) + target.parent.mkdir(parents=True, exist_ok=True) + extracted = tar.extractfile(member) + if extracted is None: + skipped += 1 + continue + with extracted as src, open(target, "wb") as dst: + shutil.copyfileobj(src, dst) + kept += 1 + kept_bytes += member.size + log(f"Extracted {kept} files ({kept_bytes / 1e9:.2f} GB); skipped {skipped} members") + + +######################################### +# Extract raw + morphology archives # +######################################### + +log("Extract tar.gz of flat files and cell labels") +FLAT_EXTRACTED = TMP_DIR / "input_raw" +extract_targz(par["input_raw"], FLAT_EXTRACTED) + +# DATA_DIR is the directory that holds the flat CSVs + CellLabels (the "*-Flat_files_and_images" +# folder). Locate it via the tx flat file rather than hard-coding the sample-name nesting. +tx_files = sorted(FLAT_EXTRACTED.rglob("*_tx_file.csv")) +assert tx_files, ( + f"No '*_tx_file.csv' found under {FLAT_EXTRACTED}. " + f"Contents: {[p.name for p in FLAT_EXTRACTED.rglob('*')][:50]}" +) +DATA_DIR = tx_files[0].parent +log(f"Detected flat-files directory: {DATA_DIR}") + +# Check that all flat files are present +FLAT_FILES_ENDINGS = ["_exprMat_file.csv", "_fov_positions_file.csv", "_metadata_file.csv", "_tx_file.csv"] for ending in FLAT_FILES_ENDINGS: - if any(f.endswith(ending) for f in os.listdir(DATA_DIR)): - print(f"Flat file with ending {ending} is present", flush=True) - else: - print(f"Flat file with ending {ending} is missing", flush=True) - -# Download image files -print(datetime.now() - t0, "Download image files", flush=True) -os.system(f"wget {URL_MORPHOLOGY} -O '{FILE_NAME_MORPHOLOGY}'") - -# Extract tar.gz files -print(datetime.now() - t0, "Extract tar.gz of image files", flush=True) -with tarfile.open(FILE_NAME_MORPHOLOGY, "r:gz") as tar: - tar.extractall(TMP_DIR) - -# Move image files to DATA_DIR / "Morphology2D" -print(datetime.now() - t0, f"Move image files of plane Z004 to {DATA_DIR / 'Morphology2D'} and rename.", flush=True) -(DATA_DIR / "Morphology2D").mkdir(parents=True, exist_ok=True) -for image_file in MORPHOLOGY_DIR.glob("*_Z004.TIF"): - shutil.move(image_file, DATA_DIR / "Morphology2D" / image_file.name.replace("_Z004.TIF", ".tif")) + present = any(f.name.endswith(ending) for f in DATA_DIR.iterdir()) + log(f"Flat file with ending {ending} is {'present' if present else 'MISSING'}") +log("Extract tar.gz of morphology images (Z004 plane only)") +MORPH_EXTRACTED = TMP_DIR / "input_morphology" +extract_targz(par["input_morphology"], MORPH_EXTRACTED, wanted=_member_wanted_morphology) +# Move the Z004 morphology tifs into DATA_DIR / "Morphology2D", renamed so sopa reads one 2D +# image per FOV (e.g. _F001_Z004.TIF -> _F001.tif). +log(f"Move Z004 morphology images to {DATA_DIR / 'Morphology2D'} and rename") +(DATA_DIR / "Morphology2D").mkdir(parents=True, exist_ok=True) +z004_tifs = sorted(MORPH_EXTRACTED.rglob("*_Z004.[Tt][Ii][Ff]")) +assert z004_tifs, f"No '*_Z004.TIF' morphology images found under {MORPH_EXTRACTED}" +for image_file in z004_tifs: + dest_name = _Z_PLANE_RE.sub(".tif", image_file.name) + shutil.move(str(image_file), DATA_DIR / "Morphology2D" / dest_name) +log(f"Moved {len(z004_tifs)} morphology images") + +# sopa expects a CellLabels/ folder. The NSCLC flat archive ships one directly, but some CosMx +# exports only ship per-FOV label tifs inside FOV* folders — handle both, matching the general +# bruker_cosmx loader. See https://github.com/gustaveroussy/sopa/issues/285 +labels_dir = DATA_DIR / "CellLabels" +fov_label_tifs = sorted(DATA_DIR.glob("FOV*/CellLabels_F*.tif")) +if labels_dir.is_dir(): + log("CellLabels folder already present in raw data") +elif fov_label_tifs: + log(f"Symlink {len(fov_label_tifs)} per-FOV CellLabels tifs into {labels_dir}") + labels_dir.mkdir(parents=True, exist_ok=True) + for tif in fov_label_tifs: + os.symlink(tif.resolve(), labels_dir / tif.name) +else: + raise AssertionError(f"No CellLabels folder and no per-FOV CellLabels tifs found in {DATA_DIR}") + +############################################# +# Memory optimization: lean transcript read # +############################################# +# sopa's CosMx reader loads the ENTIRE `*_tx_file.csv` into a pandas DataFrame in one +# `pd.read_csv` (sopa/io/reader/cosmx.py::_CosMXReader.read_transcripts). The string columns +# ('target', 'CellComp', ...) come in as object dtype (~60 B/row) — the peak of the loader. +# We wrap the `read_csv` used inside sopa's cosmx module so that, *only for the transcripts +# file*, we (a) read the 'target' feature column as categorical and (b) keep just the columns +# sopa's stitched-FOV path (read_transcripts + _get_global_cell_id) and the downstream schema +# actually use. Every other sopa read is left untouched. (Same shim as the general loader.) +from sopa.io.reader import cosmx as _cosmx_mod + +_real_pd = _cosmx_mod.pd +_real_read_csv = _real_pd.read_csv + +_TX_KEEP = ["fov", "cell_ID", "target", "x_global_px", "y_global_px", "z"] +_TX_REQUIRED = {"target", "x_global_px", "y_global_px", "fov", "cell_ID"} + + +def _lean_read_csv(filepath_or_buffer, *args, **kwargs): + name = str(filepath_or_buffer) + if "_tx_file.csv" in name and "usecols" not in kwargs: + header = _real_read_csv( + filepath_or_buffer, nrows=0, + compression=kwargs.get("compression", "infer"), + ) + cols = set(header.columns) + if _TX_REQUIRED.issubset(cols): + kwargs["usecols"] = [c for c in _TX_KEEP if c in cols] + dtype = dict(kwargs.get("dtype") or {}) + dtype.setdefault("target", "category") + kwargs["dtype"] = dtype + log(f"Lean transcript read of {Path(name).name}: usecols={kwargs['usecols']}, target=category") + return _real_read_csv(filepath_or_buffer, *args, **kwargs) + + +class _PandasReadCsvProxy: + """Forwards every attribute to pandas, but leans the transcripts read_csv.""" + def __init__(self, real): + self.__dict__["_real"] = real + def read_csv(self, *args, **kwargs): + return _lean_read_csv(*args, **kwargs) + def __getattr__(self, item): + return getattr(self._real, item) + + +_cosmx_mod.pd = _PandasReadCsvProxy(_real_pd) ######################################### # Convert raw files to spatialdata zarr # ######################################### -# from pathlib import Path -# import sopa -# DATA_DIR = Path("/Volumes/T7/G3_temp/datasets/cosmx/test_folder/Lung9_Rep1/Lung9_Rep1-Flat_files_and_images") - -print(datetime.now() - t0, "Convert raw files to spatialdata zarr", flush=True) +log("Convert raw files to spatialdata zarr") # The tif images of the NSCLC samples miss metadata information, so we hardcode the channels here. def fixed_get_morphology_coords(images_dir: Path) -> list[str]: @@ -146,27 +258,26 @@ def fixed_get_morphology_coords(images_dir: Path) -> list[str]: sopa.io.reader.cosmx._get_morphology_coords = fixed_get_morphology_coords sdata = sopa.io.cosmx( - DATA_DIR, - dataset_id=None, - fov=None, - read_proteins=False, - cells_labels=True, - cells_table=True, - cells_polygons=False, - flip_image=False + DATA_DIR, + dataset_id=None, + fov=None, + read_proteins=False, + cells_labels=True, + cells_table=True, + cells_polygons=False, # NSCLC flat files ship no -polygons.csv; derived from labels below + flip_image=False, ) # Retrieve polygons from segmentation (no polygons file in flat files present for NSCLC samples) +log("Derive cell polygons from stitched labels") sdata["cells_polygons"] = sd.to_polygons(sdata["stitched_labels"]) - ############### # Rename keys # ############### -print(datetime.now() - t0, "Rename keys", flush=True) - +log("Rename keys") elements_renaming_map = { - "stitched_image" : "morphology_mip", + "stitched_image" : "morphology_mip", "stitched_labels" : "cell_labels", "points" : "transcripts", "cells_polygons" : "cell_boundaries", @@ -174,38 +285,45 @@ def fixed_get_morphology_coords(images_dir: Path) -> list[str]: } for old_key, new_key in elements_renaming_map.items(): + if old_key not in sdata: + log(f"Element '{old_key}' not present, skipping rename to '{new_key}'") + continue sdata[new_key] = sdata[old_key] del sdata[old_key] -# Rename transcript column (somehow overwriting the 'target' column leads to an error, so instead we add a duplicate with the right name) -#sdata['transcripts'] = sdata['transcripts'].rename(columns={"global_cell_id":"cell_id", "target":"feature_name"}) -sdata['transcripts'] = sdata['transcripts'].rename(columns={"global_cell_id":"cell_id"}) +# Add transcript columns with the names expected downstream. +sdata['transcripts']["cell_id"] = sdata['transcripts']["global_cell_id"] sdata['transcripts']["feature_name"] = sdata['transcripts']["target"] ######################################### # Throw out all channels except of DAPI # ######################################### -print(datetime.now() - t0, "Throw out all channels except of 'DAPI'", flush=True) +log("Throw out all channels except of 'DAPI'") # TODO: in the future we want to keep PolyT and Cellbound1/2/3 stains. Note however, that somehow saving or plotting the sdata fails when # these channels aren't excluded, not sure why... sdata["morphology_mip"] = sdata["morphology_mip"].sel(c=["DAPI"]) - ############################## # Add info to metadata table # ############################## -print(datetime.now() - t0, "Add metadata to table", flush=True) - -#TODO: values as input variables +log("Add metadata to table") +# The common iST API requires a per-cell unique `cell_id` in obs. sopa keys the table by the +# global cell id (the index), so expose that as `cell_id`. But the CosMx metadata already +# carries a per-FOV-*local* `cell_ID` column, and obs cannot hold two keys that differ only in +# case ('cell_ID' vs 'cell_id') — spatialdata's zarr writer rejects that as an invalid name. +# Rename the vendor local id out of the way first. +if "cell_ID" in sdata["table"].obs.columns: + sdata["table"].obs.rename(columns={"cell_ID": "cell_ID_local"}, inplace=True) +sdata["table"].obs["cell_id"] = sdata["table"].obs.index.astype(str) for key in ["dataset_id", "dataset_name", "dataset_url", "dataset_reference", "dataset_summary", "dataset_description", "dataset_organism", "segmentation_id"]: sdata["table"].uns[key] = par[key] ######### # Write # ######### -print(datetime.now() - t0, f"Writing to {par['output']}", flush=True) +log(f"Writing to {par['output']}") sdata.write(par["output"]) -print(datetime.now() - t0, "Done", flush=True) \ No newline at end of file +log(f"Done in {datetime.now() - t0}") diff --git a/src/datasets/workflows/process_bruker_cosmx_nsclc/config.vsh.yaml b/src/datasets/workflows/process_bruker_cosmx_nsclc/config.vsh.yaml index 572ebd6eb..9d864cb66 100644 --- a/src/datasets/workflows/process_bruker_cosmx_nsclc/config.vsh.yaml +++ b/src/datasets/workflows/process_bruker_cosmx_nsclc/config.vsh.yaml @@ -5,9 +5,13 @@ argument_groups: - name: Inputs arguments: - type: string - name: --sample - example: "Lung9_Rep1" - description: "Sample id" + name: --input_raw + example: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+SMI+Flat+data.tar.gz" + description: "URL/path to the flat-files + cell-labels archive. A string (not a staged file) so the loader can stream it and extract only what sopa needs — see the loader's --input_raw." + - type: string + name: --input_morphology + example: "s3://openproblems-data/resources/raw_data/bruker_cosmx/Lung9_Rep1+RawMorphologyImages.tar.gz" + description: "URL/path to the raw morphology-images archive. Streamed in place; only the Z004 plane per FOV is extracted — see the loader's --input_morphology." - type: string name: --segmentation_id default: ["cell"] diff --git a/src/datasets/workflows/process_bruker_cosmx_nsclc/main.nf b/src/datasets/workflows/process_bruker_cosmx_nsclc/main.nf index 84dd78b09..1a32d3ca4 100644 --- a/src/datasets/workflows/process_bruker_cosmx_nsclc/main.nf +++ b/src/datasets/workflows/process_bruker_cosmx_nsclc/main.nf @@ -20,7 +20,8 @@ workflow run_wf { | bruker_cosmx_nsclc.run( fromState: [ - "sample", + "input_raw", + "input_morphology", "segmentation_id", "dataset_id", "dataset_name", diff --git a/src/methods_transcript_assignment/fastreseg/input.py b/src/methods_transcript_assignment/fastreseg/input.py index 09e59ab37..8279a4d65 100644 --- a/src/methods_transcript_assignment/fastreseg/input.py +++ b/src/methods_transcript_assignment/fastreseg/input.py @@ -129,9 +129,17 @@ def parse_arguments(): y_coords = transcripts.y.compute().to_numpy(dtype=np.int64) x_coords = transcripts.x.compute().to_numpy(dtype=np.int64) if isinstance(sdata_segm["segmentation"], xr.DataTree): - label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() + label_image = sdata_segm["segmentation"]["scale0"].image.to_numpy() else: label_image = sdata_segm["segmentation"].to_numpy() +# 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 +# raise IndexError. Matches basic_transcript_assignment; edge/background at the border. +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) cell_id_dask_series = dask.dataframe.from_dask_array( dask.array.from_array( label_image[y_coords, x_coords], chunks=tuple(sdata['transcripts'].map_partitions(len).compute())