|
| 1 | +# --- |
| 2 | +# jupyter: |
| 3 | +# jupytext: |
| 4 | +# formats: ipynb,py:light |
| 5 | +# text_representation: |
| 6 | +# extension: .py |
| 7 | +# format_name: light |
| 8 | +# format_version: '1.5' |
| 9 | +# jupytext_version: 1.19.1 |
| 10 | +# kernelspec: |
| 11 | +# display_name: CytoDataFrame (3.13.5.final.0) |
| 12 | +# language: python |
| 13 | +# name: python3 |
| 14 | +# --- |
| 15 | + |
| 16 | +# # NF1 3D pilot cropped voxel view |
| 17 | +# |
| 18 | +# View real 3D single-cell crops (with mask overlays) from two NF1 wells in one `CytoDataFrame`. |
| 19 | +# |
| 20 | +# Most of the setup below isn't about `CytoDataFrame` -- it's working around this |
| 21 | +# data's layout, which doesn't look like typical CellProfiler output: |
| 22 | +# |
| 23 | +# - There are no `Image_FileName_*`/`Image_PathName_*` columns, so each row's raw |
| 24 | +# image path is built from the patient/well/field metadata (`load_image_set`). |
| 25 | +# - Every well's mask file has the exact same generic name (`nuclei_mask.tiff`), |
| 26 | +# so combining two wells into one table needs `stage_mask` to tell them apart |
| 27 | +# -- a single-well notebook wouldn't need this at all. |
| 28 | +# |
| 29 | +# The actual `CytoDataFrame(...)` call at the end is just a handful of arguments. |
| 30 | + |
| 31 | +# + |
| 32 | +import os |
| 33 | +import pathlib |
| 34 | +import re |
| 35 | +import tempfile |
| 36 | + |
| 37 | +import pandas as pd |
| 38 | + |
| 39 | +from cytodataframe import CytoDataFrame |
| 40 | + |
| 41 | +# Local mount point for bandicoot. Set BANDICOOT_LOCAL_BASE if it is mounted |
| 42 | +# somewhere else on your machine. |
| 43 | +BANDICOOT_LOCAL_BASE = pathlib.Path( |
| 44 | + os.environ.get("BANDICOOT_LOCAL_BASE", "~/mnt/bandicoot") |
| 45 | +).expanduser() |
| 46 | +DATA_DIR = BANDICOOT_LOCAL_BASE / "NF1_organoid_data/data" |
| 47 | +PROFILES_DIR = ( |
| 48 | + DATA_DIR / "image_based_profiles_production_zedprofiler/ibp/sc_profiles_related" |
| 49 | +) |
| 50 | +if not PROFILES_DIR.is_dir(): |
| 51 | + raise FileNotFoundError( |
| 52 | + f"Profiles not found at {PROFILES_DIR}. Set the BANDICOOT_LOCAL_BASE " |
| 53 | + "environment variable to wherever bandicoot is mounted on this machine." |
| 54 | + ) |
| 55 | + |
| 56 | +# An "image set" here is one well + field-of-view (one 3D image). Add every |
| 57 | +# image set you want to view to this list -- e.g. all wells/fields in a |
| 58 | +# plate for per-plate QC -- not just the two shown here. |
| 59 | +IMAGE_SETS = [ |
| 60 | + "NF0055_T1__NF0055_T1__B10__F1", |
| 61 | + "NF0014_T1__NF0014_T1__C4__F2", |
| 62 | +] |
| 63 | +CHANNEL = "DNA" |
| 64 | +CHANNEL_CODE = "405" |
| 65 | +COMPARTMENT = "Nuclei" |
| 66 | + |
| 67 | +# Every well's mask file has the same generic name (e.g. "nuclei_mask.tiff"), |
| 68 | +# but CytoDataFrame's data_mask_context_dir only matches masks by filename |
| 69 | +# pattern within one shared directory -- it can't tell two identically-named |
| 70 | +# files in different wells apart. MASK_LINK_DIR is a scratch directory (not |
| 71 | +# part of the repo) that `stage_mask` below fills with per-well symlinks to |
| 72 | +# the real mask files, renamed to embed each well's own identifier, so that |
| 73 | +# matching works. |
| 74 | +MASK_LINK_DIR = pathlib.Path(tempfile.gettempdir()) / "cytodataframe_nf1_3d_mask_links" |
| 75 | +MASK_LINK_DIR.mkdir(exist_ok=True) |
| 76 | + |
| 77 | + |
| 78 | +def load_image_set(image_set: str) -> tuple[pd.DataFrame, pathlib.Path, pathlib.Path]: |
| 79 | + profiles = pd.read_parquet(PROFILES_DIR / f"{image_set}.parquet").head(1) |
| 80 | + row = profiles.iloc[0] |
| 81 | + patient = row["Metadata_Biology_PatientTumor"] |
| 82 | + # e.g. B10-1 |
| 83 | + well_field = ( |
| 84 | + f"{row['Metadata_Experiment_WellID']}-{row['Metadata_Imaging_FieldID']}" |
| 85 | + ) |
| 86 | + |
| 87 | + channel_path = ( |
| 88 | + DATA_DIR |
| 89 | + / patient |
| 90 | + / "zstack_images" |
| 91 | + / well_field |
| 92 | + / f"{well_field}_{CHANNEL_CODE}.tif" |
| 93 | + ) |
| 94 | + mask_path = ( |
| 95 | + DATA_DIR |
| 96 | + / patient |
| 97 | + / "segmentation_masks" |
| 98 | + / well_field |
| 99 | + / f"{COMPARTMENT.lower()}_mask.tiff" |
| 100 | + ) |
| 101 | + profiles[f"Image_FileName_{CHANNEL}"] = str(channel_path) |
| 102 | + return profiles, channel_path, mask_path |
| 103 | + |
| 104 | + |
| 105 | +def stage_mask(channel_path: pathlib.Path, mask_path: pathlib.Path) -> pathlib.Path: |
| 106 | + # See MASK_LINK_DIR above for why this exists. |
| 107 | + link = MASK_LINK_DIR / f"{channel_path.stem}__{mask_path.name}" |
| 108 | + if not link.exists(): |
| 109 | + link.symlink_to(mask_path) |
| 110 | + return link |
| 111 | + |
| 112 | + |
| 113 | +# - |
| 114 | + |
| 115 | +# Load one object per well into a single two-row table. |
| 116 | + |
| 117 | +# + |
| 118 | +profile_rows = [] |
| 119 | +for image_set in IMAGE_SETS: |
| 120 | + profiles, channel_path, mask_path = load_image_set(image_set) |
| 121 | + stage_mask(channel_path, mask_path) |
| 122 | + profile_rows.append(profiles) |
| 123 | + |
| 124 | +profiles = pd.concat(profile_rows, ignore_index=True) |
| 125 | +mask_name = mask_path.name |
| 126 | + |
| 127 | +# + |
| 128 | +bbox_column_map = { |
| 129 | + "x_min": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MinX", |
| 130 | + "x_max": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MaxX", |
| 131 | + "y_min": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MinY", |
| 132 | + "y_max": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MaxY", |
| 133 | + "z_min": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MinZ", |
| 134 | + "z_max": f"{COMPARTMENT}_NoChannel_VolumeSizeShape_MaxZ", |
| 135 | +} |
| 136 | + |
| 137 | +center_columns = [ |
| 138 | + f"{COMPARTMENT}_NoChannel_VolumeSizeShape_Center{axis}" for axis in "XYZ" |
| 139 | +] |
| 140 | + |
| 141 | +voxel_view = CytoDataFrame( |
| 142 | + data=profiles[[f"Image_FileName_{CHANNEL}"]], |
| 143 | + data_bounding_box=profiles[list(bbox_column_map.values())], |
| 144 | + compartment_center_xy=profiles[center_columns], |
| 145 | + data_mask_context_dir=str(MASK_LINK_DIR), |
| 146 | + segmentation_file_regex={rf"__{re.escape(mask_name)}$": r"_\d+\.tif$"}, |
| 147 | + display_options={ |
| 148 | + "width": 260, |
| 149 | + "height": 260, |
| 150 | + "table_max_height": "580px", |
| 151 | + "label_overlay_mode": "filled", |
| 152 | + # Voxel size (x, y, z) in um. These TIFFs carry no voxel-size |
| 153 | + # metadata, and z is coarser than xy, so without this the objects are |
| 154 | + # drawn too flat. xy_spacing (0.1) is the NF1 3D pipeline manifests' |
| 155 | + # recorded value; z_spacing (0.5) was chosen visually against these |
| 156 | + # two objects rather than read from the manifests (which record 1.0). |
| 157 | + "volume_spacing": (0.1, 0.1, 0.5), |
| 158 | + "volume_bbox_column_map": bbox_column_map, |
| 159 | + }, |
| 160 | +) |
| 161 | +# The rendered widget includes a "Mask" checkbox to toggle the overlay on/off, |
| 162 | +# and a red dot marking each object's center. |
| 163 | +voxel_view |
0 commit comments