Skills › Scientific Agent Skills › cellprofiler › nuclei_assay.py
nuclei_assay.py
216 lines · skills/cellprofiler/scripts/nuclei_assay.py at 5f2bc03
1#!/usr/bin/env python3
2"""Validate a microscopy manifest, run CellProfiler, and check nuclei measurements."""
3from __future__ import annotations
5import argparse
6import csv
7from decimal import Decimal, InvalidOperation
8import hashlib
9import json
10import math
11import re
12from pathlib import Path
13import subprocess
16def _sha256(path: Path) -> str:
17 digest = hashlib.sha256()
18 with path.open("rb") as handle:
19 for chunk in iter(lambda: handle.read(1024 * 1024), b""):
20 digest.update(chunk)
21 return digest.hexdigest()
24def _csv_rows(path: Path, required: set[str]) -> list[dict]:
25 with path.open(newline="", encoding="utf-8-sig") as handle:
26 reader = csv.DictReader(handle)
27 names = reader.fieldnames or []
28 if len(names) != len(set(names)) or not required.issubset(names):
29 raise ValueError(f"{path.name}: requires unique columns including {sorted(required)}")
30 rows = list(reader)
31 for row in rows:
32 if None in row or any(value is None for value in row.values()):
33 raise ValueError(f"{path.name}: row length differs from header")
34 if any(not row[key].strip() for key in required):
35 raise ValueError(f"{path.name}: empty required value")
36 return rows
39def prepare(manifest: Path, destination: Path) -> dict:
40 """Make LoadData CSV without altering image pixels."""
41 import numpy as np
42 import tifffile
44 rows = _csv_rows(manifest, {"sample_id", "image_path", "plate", "well", "site"})
45 if not rows:
46 raise ValueError("Manifest contains no images")
47 output, qc, seen, fields_seen, paths_seen = [], [], set(), set(), set()
48 for row in rows:
49 if any(row[key] != row[key].strip() for key in ("sample_id", "plate", "well", "site")):
50 raise ValueError("Manifest identifiers must not contain surrounding whitespace")
51 sample = row["sample_id"]
52 if not re.fullmatch(r"[A-Za-z0-9][A-Za-z0-9_.-]*", sample):
53 raise ValueError("Sample IDs must be filename-safe letters, digits, dots, dashes or underscores")
54 try:
55 float(sample)
56 except ValueError:
57 pass
58 else:
59 raise ValueError("Sample IDs must not be numeric; use a prefix such as sample_001")
60 field = (row["plate"], row["well"], row["site"])
61 if sample in seen or field in fields_seen:
62 raise ValueError(f"Duplicate sample or plate/well/site: {sample}")
63 seen.add(sample)
64 fields_seen.add(field)
65 path = Path(row["image_path"])
66 path = (manifest.parent / path).resolve() if not path.is_absolute() else path.resolve()
67 if path in paths_seen:
68 raise ValueError(f"Repeated image file: {path}")
69 paths_seen.add(path)
70 with tifffile.TiffFile(path) as tif:
71 if len(tif.series) != 1 or len(tif.pages) != 1 or len(tif.series[0].levels) != 1:
72 raise ValueError(f"{sample}: requires a single TIFF plane, series and resolution level")
73 if tif.pages[0].photometric != 1 or tif.pages[0].samplesperpixel != 1:
74 raise ValueError(f"{sample}: requires black-is-zero grayscale TIFF (MINISBLACK)")
75 image = tif.asarray()
76 if image.ndim != 2 or image.dtype not in (np.dtype("uint8"), np.dtype("uint16")):
77 raise ValueError(f"{sample}: requires single-plane grayscale uint8/uint16 TIFF")
78 if min(image.shape) < 16 or np.ptp(image) == 0:
79 raise ValueError(f"{sample}: image too small or constant")
80 maximum = np.iinfo(image.dtype).max
81 qc.append({"sample_id": sample, "shape": list(image.shape), "dtype": str(image.dtype),
82 "saturated_fraction": float(np.mean(image == maximum)),
83 "sha256": _sha256(path)})
84 output.append({"URL_DNA": path.as_uri(), "Metadata_Sample": sample,
85 "Metadata_Plate": row["plate"], "Metadata_Well": row["well"],
86 "Metadata_Site": row["site"]})
87 if destination.resolve() in paths_seen | {manifest.resolve()}:
88 raise ValueError("LoadData destination must not overwrite the manifest or an input image")
89 destination.parent.mkdir(parents=True, exist_ok=True)
90 with destination.open("w", newline="", encoding="utf-8") as handle:
91 writer = csv.DictWriter(handle, fieldnames=list(output[0]))
92 writer.writeheader()
93 writer.writerows(output)
94 return {"images": qc, "intensity_scaling": "integer dtype maximum (255 or 65535)",
95 "warnings": [f"{x['sample_id']}: pixels at storage maximum" for x in qc if x["saturated_fraction"] > 0.01]}
98def _integer(value: str, name: str, minimum: int = 1) -> int:
99 try:
100 number = Decimal(value)
101 except InvalidOperation as exc:
102 raise ValueError(f"{name} must be an integer") from exc
103 if not number.is_finite() or number != number.to_integral_value() or number < minimum:
104 raise ValueError(f"{name} must be a finite integer >= {minimum}")
105 return int(number)
108def summarize(directory: Path, expected_samples: list[str] | None = None) -> dict:
109 images = _csv_rows(directory / "Image.csv", {"ImageNumber", "Metadata_Sample", "Count_Nuclei"})
110 objects = _csv_rows(directory / "Nuclei.csv", {
111 "ImageNumber", "ObjectNumber", "Intensity_MeanIntensity_DNA",
112 "Intensity_IntegratedIntensity_DNA", "AreaShape_Area",
113 })
114 if not images:
115 raise ValueError("CellProfiler produced no image measurements")
116 actual_samples = [r["Metadata_Sample"] for r in images]
117 if len(actual_samples) != len(set(actual_samples)):
118 raise ValueError("Duplicate image metadata in output")
119 if expected_samples is not None and set(actual_samples) != set(expected_samples):
120 raise ValueError("Output sample identities differ from the input manifest")
121 image_numbers = [_integer(row["ImageNumber"], "ImageNumber") for row in images]
122 if len(set(image_numbers)) != len(image_numbers):
123 raise ValueError("Duplicate ImageNumber in output")
124 by_image = {number: {} for number in image_numbers}
125 for obj in objects:
126 image_number = _integer(obj["ImageNumber"], "ImageNumber")
127 object_number = _integer(obj["ObjectNumber"], "ObjectNumber")
128 if image_number not in by_image:
129 raise ValueError("Orphan object row")
130 if object_number in by_image[image_number]:
131 raise ValueError("Duplicate ImageNumber/ObjectNumber pair")
132 by_image[image_number][object_number] = obj
133 measurements, warnings = [], []
134 for row, image_number in zip(images, image_numbers):
135 if any(float(value) != 0 for key, value in row.items() if key.startswith("ModuleError_")):
136 raise ValueError("CellProfiler reported a failed module; inspect cellprofiler.log")
137 subset = list(by_image[image_number].values())
138 count = _integer(row["Count_Nuclei"], "Count_Nuclei", minimum=0)
139 if count != len(subset) or set(by_image[image_number]) != set(range(1, count + 1)):
140 raise ValueError("Image count and per-object rows disagree")
141 intensities = [float(obj["Intensity_MeanIntensity_DNA"]) for obj in subset]
142 if any(not math.isfinite(v) or not 0 <= v <= 1 for v in intensities):
143 raise ValueError("Nonfinite or out-of-range normalized intensities")
144 areas = [float(obj["AreaShape_Area"]) for obj in subset]
145 integrated = [float(obj["Intensity_IntegratedIntensity_DNA"]) for obj in subset]
146 if any(not math.isfinite(v) or v <= 0 for v in areas):
147 raise ValueError("Nonfinite or nonpositive nuclear areas")
148 if any(not math.isfinite(v) or v < 0 for v in integrated):
149 raise ValueError("Nonfinite or negative integrated intensities")
150 if count == 0:
151 warnings.append(f"{row['Metadata_Sample']}: zero nuclei; inspect overlay and threshold")
152 measurements.append({"sample_id": row["Metadata_Sample"], "nuclei": count,
153 "mean_nuclear_intensity": sum(intensities) / count if count else None,
154 "mean_nuclear_area_pixels": sum(areas) / count if count else None})
155 return {"measurements": measurements, "warnings": warnings}
158def run(manifest: Path, output: Path, executable: str, pipeline: Path) -> dict:
159 if output.exists() and any(output.iterdir()):
160 raise ValueError("Use an empty output directory to keep runs distinct")
161 output.mkdir(parents=True, exist_ok=True)
162 qc = prepare(manifest.resolve(), output / "load_data.csv")
163 snapshot = output / "pipeline.cppipe"
164 snapshot.write_bytes(pipeline.read_bytes())
165 done_file = output / "cellprofiler.done"
166 command = [executable, "-c", "-r", "-p", str(snapshot.resolve()), "--data-file",
167 str((output / "load_data.csv").resolve()), "-o", str(output.resolve()),
168 "--done-file", str(done_file.resolve())]
169 result = {"status": "running", "input_qc": qc, "command": command,
170 "pipeline_sha256": _sha256(snapshot)}
171 report = output / "assay_qc.json"
172 report.write_text(json.dumps(result, indent=2) + "\n")
173 try:
174 with (output / "cellprofiler.log").open("w") as log:
175 subprocess.run(command, stdout=log, stderr=subprocess.STDOUT, check=True)
176 if not done_file.is_file() or done_file.read_text().strip() != "Complete":
177 raise ValueError("CellProfiler completion marker is missing or not Complete")
178 result.update(summarize(output, [r["sample_id"] for r in qc["images"]]))
179 result["status"] = "complete"
180 except (ValueError, OSError, KeyError, csv.Error, subprocess.CalledProcessError) as exc:
181 result.update({"status": "failed", "error": str(exc)})
182 raise
183 finally:
184 report.write_text(json.dumps(result, indent=2) + "\n")
185 return result
188def main() -> None:
189 parser = argparse.ArgumentParser(description=__doc__)
190 sub = parser.add_subparsers(dest="action", required=True)
191 prep = sub.add_parser("prepare", help="Validate CSV manifest and create LoadData CSV")
192 prep.add_argument("manifest", type=Path)
193 prep.add_argument("output_csv", type=Path)
194 analysis = sub.add_parser("run", help="Execute CellProfiler and check its exported measurements")
195 analysis.add_argument("manifest", type=Path)
196 analysis.add_argument("output", type=Path)
197 analysis.add_argument("--executable", default="cellprofiler")
198 analysis.add_argument("--pipeline", type=Path,
199 default=Path(__file__).resolve().parents[1] / "assets" / "nuclei.cppipe")
200 summary = sub.add_parser("summarize", help="Check existing Image.csv and Nuclei.csv")
201 summary.add_argument("output", type=Path)
202 args = parser.parse_args()
203 try:
204 if args.action == "prepare":
205 result = prepare(args.manifest, args.output_csv)
206 elif args.action == "run":
207 result = run(args.manifest, args.output, args.executable, args.pipeline)
208 else:
209 result = summarize(args.output)
210 except (ValueError, OSError, KeyError, csv.Error, subprocess.CalledProcessError) as exc:
211 parser.exit(1, f"Error: {exc}\n")
212 print(json.dumps(result, indent=2))
215if __name__ == "__main__":
216 main()