diff --git a/CHANGELOG.md b/CHANGELOG.md index fb1f011..e5c5e54 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,21 @@ # Changelog +## Unreleased + +- Added `decode_epu_xml`: reads the acquisition settings out of a Thermo Fisher + EPU image XML (voltage, pixel size as recorded, binning, super-resolution + factor, fractions, exposure, dose, applied defocus, detector, energy-filter + slit). Not registered for automatic dispatch, because `.xml` is too generic an + extension to claim. `epu_pixel_size_candidates` lists the plausible pixel sizes + for the saved movie, since EPU records the pixel size of the binned grid. +- Added `compare_acquisition` and `Finding`: cross-check voltage, pixel size, dose, + frame count, Cs and magnification across any number of sources (a collection + sheet, EPU, a movie header, a RELION or CryoSPARC optics table, a deposit) and + report disagreements with a severity. With `target_resolution_a` it also flags a + pixel size whose Nyquist limit cannot reach the target. Motivated by a processing + run that used 300 kV for data recorded at 200 kV and silently stalled several + angstrom short of the published resolution. + ## 0.4.1 - Added `[project.urls]` (Homepage, Repository, Issues) — the prior release had diff --git a/README.md b/README.md index 93c457c..ba5500a 100644 --- a/README.md +++ b/README.md @@ -98,6 +98,29 @@ a format. NPY and NIfTI landed with zero plumbing changes. ## Pixel size for raw movies +### Cross-check the acquisition settings + +Voltage, pixel size, dose and frame count usually arrive from several places typed +or exported by different people. A wrong one does not fail; it quietly caps the +resolution. `compare_acquisition` lines the sources up and reports disagreements: + +```python +from scigantic_headers import compare_acquisition, decode_epu_xml, read_star_optics + +epu = decode_epu_xml(open("FoilHole_..._Data_....xml", "rb").read()).fields +for f in compare_acquisition({ + "sheet": {"voltageKv": 300, "pixelSizeA": 0.84, "doseEA2": 42.4, "fractions": 40}, + "epu": epu, + "relion": read_star_optics("Import/movies.star"), +}, target_resolution_a=3.3): + print(f.severity, f.field, f.message) +``` + +It does not decide which source is right. EPU records the pixel size of the binned +readout grid, so `epu_pixel_size_candidates(epu)` lists the likely values for the +saved movie; compare with the movie's own dimensions. `decode_epu_xml` is not +auto-dispatched (`.xml` is too generic), so call it on a file you know is EPU's. + A raw MRC movie header has no pixel size: CELLA is 0, so `decode_*` returns `pixelSizeA` as `None`. The value is in the RELION STAR or CryoSPARC `.cs` file the workflow writes in the session, and the library reads it. `read_session_optics` diff --git a/src/scigantic_headers/__init__.py b/src/scigantic_headers/__init__.py index 5658884..12ec96c 100644 --- a/src/scigantic_headers/__init__.py +++ b/src/scigantic_headers/__init__.py @@ -11,9 +11,11 @@ from importlib.metadata import PackageNotFoundError, version as _pkg_version +from .acquisition import Finding, compare_acquisition from .bed import decode_bed from .cryosparc import parse_cryosparc_optics, read_cryosparc_optics from .dicom import decode_dicom +from .epu import decode_epu_xml, epu_pixel_size_candidates from .fastq import decode_fastq from .fcs import decode_fcs from .genbank import decode_genbank @@ -59,6 +61,7 @@ except PackageNotFoundError: # a source checkout with no install metadata __version__ = "0.0.0+unknown" __all__ = [ + "Finding", "compare_acquisition", "decode_epu_xml", "epu_pixel_size_candidates", "HEADER_BYTES", "DecodedHeader", "Read", "decode_bytes", "decode_mrc_header", "decode_nifti_header", "decode_npy_header", "decode_cryosparc_header", "parse_cryosparc_optics", "read_cryosparc_optics", "read_session_optics", diff --git a/src/scigantic_headers/acquisition.py b/src/scigantic_headers/acquisition.py new file mode 100644 index 0000000..bac3c36 --- /dev/null +++ b/src/scigantic_headers/acquisition.py @@ -0,0 +1,157 @@ +"""Cross-check acquisition settings that come from different places. + +A cryo-EM processing run needs voltage, pixel size, dose, frame count and Cs. They +usually arrive from several sources that were each typed or exported by someone +else: a collection sheet, the EPU or SerialEM file, the movie header, the RELION +or CryoSPARC optics table, the EMDB deposit. When two of them disagree, one is +wrong, and the wrong one silently caps the resolution rather than failing: a +processing run that used 300 kV for data recorded at 200 kV ran to the end and +simply stalled several angstrom short of the published map. + +`compare_acquisition` takes the values from any number of sources, finds the fields +where they disagree beyond a tolerance, and says how serious it is. It does not +decide which source is right; it makes the disagreement visible before compute is +spent. Sources are plain dicts using the key names this library already returns +(`pixelSizeA`, `voltageKv` from `read_star_optics` and the MRC decoder; the rest from +`decode_epu_xml`), so they can be passed straight in. + + from scigantic_headers import compare_acquisition, read_star_optics + findings = compare_acquisition({ + "sheet": {"voltageKv": 300, "pixelSizeA": 0.84, "doseEA2": 42.4, "fractions": 40}, + "relion": read_star_optics("Import/movies.star"), + "deposit": {"voltageKv": 200, "pixelSizeA": 0.88}, + }, target_resolution_a=3.3) + for f in findings: + print(f.severity, f.field, f.message) +""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass, field +from typing import Dict, List, Mapping, Optional + +from .epu import epu_pixel_size_candidates + +_SEVERITY_ORDER = {"error": 0, "warn": 1, "info": 2} + + +@dataclass(frozen=True) +class Finding: + """One disagreement (or risk) in the acquisition settings.""" + + field: str + severity: str # "error" | "warn" | "info" + message: str + values: Dict[str, object] = field(default_factory=dict) + + def to_dict(self) -> dict: + return asdict(self) + + +def _present(sources: Mapping[str, Mapping[str, object]], key: str) -> Dict[str, float]: + out = {} + for name, src in sources.items(): + v = src.get(key) if src else None + if isinstance(v, (int, float)) and not isinstance(v, bool) and v == v and v > 0: + out[name] = float(v) + return out + + +def _fmt(values: Mapping[str, float]) -> str: + return ", ".join("%s=%g" % (k, v) for k, v in values.items()) + + +def _spread(values: Mapping[str, float]) -> float: + lo, hi = min(values.values()), max(values.values()) + return hi / lo - 1.0 + + +def compare_acquisition( + sources: Mapping[str, Mapping[str, object]], + *, + target_resolution_a: Optional[float] = None, +) -> List[Finding]: + """Compare acquisition values across sources. Returns findings, worst first. + + An empty list means every field present in two or more sources agrees within + tolerance (and, if `target_resolution_a` is given, the pixel size can reach it). + A field only one source reports cannot be checked and is not flagged.""" + findings: List[Finding] = [] + + volts = _present(sources, "voltageKv") + if len(volts) >= 2 and max(volts.values()) - min(volts.values()) > 0.5: + findings.append(Finding( + "voltageKv", "error", + "voltage differs between sources (%s). A wrong voltage corrupts the CTF and " + "the dose weighting at high resolution without any error being raised." % _fmt(volts), + dict(volts))) + + pix = _present(sources, "pixelSizeA") + if len(pix) >= 2: + s = _spread(pix) + if s > 0.01: + hint = "" + for name, src in sources.items(): + for label, cand in epu_pixel_size_candidates(dict(src or {})).items(): + if label == "asRecorded": + continue + for other_name, other in pix.items(): + # Calibrated and nominal pixel sizes routinely differ by a few + # percent, so match loosely; this is a pointer, not a verdict. + if other_name != name and abs(cand / other - 1) < 0.06: + hint = (" %s=%g is within %.1f%% of %g, the %s pixel size derived from '%s' " + "(recorded %g). Check which grid the movies are saved on from " + "their own dimensions." + % (other_name, other, 100 * abs(cand / other - 1), cand, label, name, + src.get("pixelSizeA"))) + break + if hint: + break + if hint: + break + findings.append(Finding( + "pixelSizeA", "error" if s > 0.05 else "warn", + "pixel size differs by %.1f%% (%s).%s" % (100 * s, _fmt(pix), hint), dict(pix))) + + dose = _present(sources, "doseEA2") + if len(dose) >= 2: + s = _spread(dose) + if s > 0.10: + findings.append(Finding( + "doseEA2", "error" if s > 0.25 else "warn", + "total dose differs by %.0f%% (%s). Dose weighting uses the per-frame dose." % (100 * s, _fmt(dose)), + dict(dose))) + + frames = _present(sources, "fractions") + if len(frames) >= 2 and len({int(v) for v in frames.values()}) > 1: + findings.append(Finding( + "fractions", "error", + "frame count differs between sources (%s)." % _fmt(frames), dict(frames))) + + cs = _present(sources, "csMm") + if len(cs) >= 2 and max(cs.values()) - min(cs.values()) > 0.1: + findings.append(Finding("csMm", "warn", "Cs differs between sources (%s)." % _fmt(cs), dict(cs))) + + mag = _present(sources, "magnification") + if len(mag) >= 2 and _spread(mag) > 0.01: + findings.append(Finding("magnification", "warn", + "nominal magnification differs (%s)." % _fmt(mag), dict(mag))) + + if target_resolution_a: + for px in sorted({round(v, 4) for v in pix.values()}): + nyq = 2 * px + if nyq > target_resolution_a: + findings.append(Finding( + "pixelSizeA", "error", + "at %g A/px the Nyquist limit is %.2f A, coarser than the %g A target; the " + "target cannot be reached on this grid." % (px, nyq, target_resolution_a), + {"pixelSizeA": px, "nyquistA": nyq})) + elif nyq * 1.1 > target_resolution_a: + findings.append(Finding( + "pixelSizeA", "warn", + "at %g A/px the Nyquist limit is %.2f A, within 10%% of the %g A target; " + "there is almost no room above it." % (px, nyq, target_resolution_a), + {"pixelSizeA": px, "nyquistA": nyq})) + + findings.sort(key=lambda f: _SEVERITY_ORDER[f.severity]) + return findings diff --git a/src/scigantic_headers/epu.py b/src/scigantic_headers/epu.py new file mode 100644 index 0000000..4e2be4d --- /dev/null +++ b/src/scigantic_headers/epu.py @@ -0,0 +1,130 @@ +"""Decode the acquisition metadata of a Thermo Fisher EPU image XML. + +EPU writes one XML next to each movie (`FoilHole_..._Data_....xml`, root element +`MicroscopeImage`). It carries the settings the movie was recorded with: voltage, +pixel size, camera, binning, number of fractions, exposure, dose, applied defocus. +These are the values a processing run needs and the ones most often copied in by +hand from somewhere else, so they are worth reading from the source. + +This is NOT registered for automatic dispatch. `.xml` is far too generic an +extension to claim (mzML and Illumina use it too), and the EPU file name carries no +fixed marker. Call `decode_epu_xml` on bytes you already know are an EPU file; it +returns None for anything without a `MicroscopeImage` root. + +Units in the file, and what is returned here: + AccelerationVoltage volts -> voltageKv (kV) + pixelSize metres -> pixelSizeA (angstrom, AS RECORDED, see below) + AppliedDefocus metres -> appliedDefocusUm (micrometres, negative = underfocus) + Dose electrons/m^2 -> doseEA2 (electrons/angstrom^2) + +`pixelSizeA` is the pixel size of the grid EPU recorded on. A camera read out with +2x2 binning reports a pixel twice the physical one, and a super-resolution factor +means the saved movie can sit on a finer grid still. Which grid a saved movie uses +depends on the export setting, so the file alone cannot say. `epu_pixel_size_candidates` +lists the plausible values; check the movie's own dimensions to choose. +""" + +from __future__ import annotations + +import re +from typing import Dict, Optional + +from .decoders import DecodedHeader + +_NS = rb"(?:\w+:)?" # EPU XML prefixes its collection elements (a:, b:); match with or without + + +def _kv(data: bytes, key: str) -> Optional[str]: + """Value of a `CustomData` key/value pair.""" + pat = (_NS + rb"Key>" + re.escape(key.encode()) + rb"\s*<" + _NS + + rb"Value[^>]*>([^<]*) Optional[float]: + try: + return float(text) if text is not None else None + except ValueError: + return None + + +def _first_float(data: bytes, pattern: bytes) -> Optional[float]: + m = re.search(pattern, data, re.S) + return _num(m.group(1).decode()) if m else None + + +def decode_epu_xml(data: bytes) -> Optional[DecodedHeader]: + if b"([^<]*)") + px_m = _first_float(data, rb"\s*\s*([^<]*)") + mag = _first_float(data, rb"([^<]*)") + fractions = _first_float(data, rb"<" + _NS + rb"NumberOffractions\b[^>]*>(\d+)<") + + # The first positive ExposureTime is the camera's; a later plateCamera one is 0. + exposure = next((v for v in (_num(x.decode()) for x in re.findall(rb"([^<]*)", data)) + if v and v > 0), None) + + bin_m = re.search(rb"]*>\s*<" + _NS + rb"x>(\d+)\s*<" + _NS + rb"y>(\d+)<", data) + binning = (int(bin_m.group(1)), int(bin_m.group(2))) if bin_m else None + ro = re.search(rb"]*>\s*<" + _NS + rb"height>(\d+)\s*<" + + _NS + rb"width>(\d+)<", data) + readout = {"height": int(ro.group(1)), "width": int(ro.group(2))} if ro else None + + superres = _num(_kv(data, "SuperResolutionFactor")) + detector = _kv(data, "DetectorCommercialName") or _kv(data, "Detectors[EF-CCD].CommercialName") + counting = _kv(data, "ElectronCountingEnabled") or _kv(data, "Detectors[EF-CCD].ElectronCounted") + defocus_m = _num(_kv(data, "AppliedDefocus")) + dose_m2 = _num(_kv(data, "Dose")) + slit = _first_float(data, rb"([^<]*)") + sw = re.search(rb"([^<]*)", data) + sw_v = re.search(rb"([^<]*)", data) + + fields: Dict[str, object] = { + "voltageKv": volts / 1000.0 if volts else None, + "pixelSizeA": px_m * 1e10 if px_m else None, + "magnification": int(mag) if mag else None, + "detector": detector, + "electronCounting": (counting or "").lower() == "true" if counting is not None else None, + "binning": list(binning) if binning else None, + "superResolutionFactor": int(superres) if superres else None, + "readoutAreaPx": readout, + "fractions": int(fractions) if fractions else None, + "exposureS": exposure, + "doseEA2": dose_m2 / 1e20 if dose_m2 else None, + "appliedDefocusUm": defocus_m * 1e6 if defocus_m is not None else None, + "energyFilterSlitEv": slit, + "phasePlate": (_kv(data, "PhasePlateUsed") or "").lower() == "true" if _kv(data, "PhasePlateUsed") is not None else None, + "software": ((sw.group(1).decode().strip() + " " + (sw_v.group(1).decode().strip() if sw_v else "")).strip() + if sw else None), + } + summary = "EPU image: %s kV, %s A/px as recorded, %s, %s fractions" % ( + ("%g" % fields["voltageKv"]) if fields["voltageKv"] else "?", + ("%.4g" % fields["pixelSizeA"]) if fields["pixelSizeA"] else "?", + detector or "unknown detector", + fields["fractions"] if fields["fractions"] else "?", + ) + return DecodedHeader(format="epu-xml", summary=summary, fields=fields) + + +def epu_pixel_size_candidates(fields: Dict[str, object]) -> Dict[str, float]: + """Plausible pixel sizes (angstrom) for the saved movie, from decoded EPU fields. + + `asRecorded` is what EPU stored. `perBinning` divides by the binning, the + physical pixel. `perBinningAndSuperRes` also divides by the super-resolution + factor, the grid of a super-resolution movie. Only candidates that differ from + `asRecorded` are added, so a plain unbinned camera returns one entry.""" + px = fields.get("pixelSizeA") + if not px: + return {} + out = {"asRecorded": float(px)} + b = fields.get("binning") + s = fields.get("superResolutionFactor") + per_bin = px / b[0] if b and b[0] and b[0] > 1 else None + if per_bin: + out["perBinning"] = per_bin + if s and s > 1: + out["perBinningAndSuperRes"] = (per_bin or px) / s + return out diff --git a/tests/test_acquisition.py b/tests/test_acquisition.py new file mode 100644 index 0000000..2c8211f --- /dev/null +++ b/tests/test_acquisition.py @@ -0,0 +1,131 @@ +"""EPU XML decoding and cross-source acquisition checks. + +The EPU document below is synthetic: it follows the structure of a real EPU image +XML but every value is made up, so no collection data lives in this public repo. +The voltage case is modelled on a real incident with public numbers: a run that used +300 kV for EMPIAR-10581 / EMD-0731, which the deposit records at 200 kV. +""" + +from scigantic_headers import ( + compare_acquisition, + decode_epu_xml, + epu_pixel_size_candidates, +) + +EPU = b""" + +Dose4.0E+21 +DetectorCommercialNameExample Detector +AppliedDefocus-2.2E-06 +SuperResolutionFactor2 +ElectronCountingEnabledtrue +PhasePlateUsedfalse + +1.92E-101.92E-10 + + +22 +2.0 +20462880 +30 +0 +EPU3.0.0.1 +300000 +20 +130000 +""" + + +def test_epu_fields(): + h = decode_epu_xml(EPU) + assert h.format == "epu-xml" + f = h.fields + assert f["voltageKv"] == 300.0 + assert abs(f["pixelSizeA"] - 1.92) < 1e-9 + assert f["binning"] == [2, 2] and f["superResolutionFactor"] == 2 + assert f["fractions"] == 30 and f["exposureS"] == 2.0 + assert abs(f["doseEA2"] - 40.0) < 1e-9 # 4.0e21 e/m^2 + assert abs(f["appliedDefocusUm"] + 2.2) < 1e-9 + assert f["electronCounting"] is True and f["phasePlate"] is False + assert f["detector"] == "Example Detector" and f["magnification"] == 130000 + assert f["energyFilterSlitEv"] == 20.0 and f["readoutAreaPx"] == {"height": 2046, "width": 2880} + assert f["software"] == "EPU 3.0.0.1" + + +def test_epu_rejects_other_xml(): + assert decode_epu_xml(b"") is None + + +def test_epu_pixel_candidates(): + f = decode_epu_xml(EPU).fields + c = epu_pixel_size_candidates(f) + assert abs(c["asRecorded"] - 1.92) < 1e-9 + assert abs(c["perBinning"] - 0.96) < 1e-9 + assert abs(c["perBinningAndSuperRes"] - 0.48) < 1e-9 + # unbinned, no super-resolution: only the recorded value + assert list(epu_pixel_size_candidates({"pixelSizeA": 1.0})) == ["asRecorded"] + + +def test_voltage_mismatch_is_an_error(): + # The rehearsal incident: sheet/pipeline said 300 kV, the deposit says 200 kV. + out = compare_acquisition({"pipeline": {"voltageKv": 300.0, "pixelSizeA": 0.8844}, + "deposit": {"voltageKv": 200.0, "pixelSizeA": 0.88}}) + assert out[0].field == "voltageKv" and out[0].severity == "error" + assert out[0].values == {"pipeline": 300.0, "deposit": 200.0} + # 0.8844 vs 0.88 is 0.5%, under the 1% tolerance, so nothing is said about it + assert [f.field for f in out] == ["voltageKv"] + + +def test_agreeing_sources_give_no_findings(): + assert compare_acquisition({ + "sheet": {"voltageKv": 300, "pixelSizeA": 0.84, "doseEA2": 42.4, "fractions": 40}, + "star": {"voltageKv": 300.0, "pixelSizeA": 0.8401}, + "header": {"pixelSizeA": 0.84}, + }) == [] + + +def test_single_source_is_not_flagged(): + assert compare_acquisition({"only": {"voltageKv": 300, "pixelSizeA": 1.0}}) == [] + + +def test_pixel_size_mismatch_points_at_the_binning_candidate(): + src = {"epu": {"pixelSizeA": 1.92, "binning": [2, 2], "superResolutionFactor": 2}, + "sheet": {"pixelSizeA": 0.97}} + out = compare_acquisition(src) + assert out[0].field == "pixelSizeA" and out[0].severity == "error" + assert "perBinning" in out[0].message and "dimensions" in out[0].message + + +def test_dose_and_frames(): + # 42.4 vs 34.0 is a 25% gap: worth saying, but dose weighting is forgiving + out = compare_acquisition({"a": {"doseEA2": 42.4, "fractions": 40}, + "b": {"doseEA2": 34.0, "fractions": 40}}) + assert [(f.field, f.severity) for f in out] == [("doseEA2", "warn")] + out = compare_acquisition({"a": {"doseEA2": 42.4}, "b": {"doseEA2": 30.0}}) + assert [(f.field, f.severity) for f in out] == [("doseEA2", "error")] + out = compare_acquisition({"a": {"fractions": 40}, "b": {"fractions": 49}}) + assert out[0].field == "fractions" and out[0].severity == "error" + + +def test_nyquist_against_target(): + out = compare_acquisition({"x": {"pixelSizeA": 1.77}}, target_resolution_a=3.3) + assert out[0].severity == "error" and "Nyquist" in out[0].message + out = compare_acquisition({"x": {"pixelSizeA": 1.55}}, target_resolution_a=3.3) + assert out[0].severity == "warn" # 3.10 A Nyquist, within 10% of 3.3 + assert compare_acquisition({"x": {"pixelSizeA": 0.88}}, target_resolution_a=3.3) == [] + + +def test_findings_are_ordered_worst_first_and_serialisable(): + out = compare_acquisition({"a": {"voltageKv": 300, "csMm": 2.7, "doseEA2": 40.0}, + "b": {"voltageKv": 200, "csMm": 2.0, "doseEA2": 40.0}}) + assert [f.severity for f in out] == ["error", "warn"] + assert out[0].to_dict()["field"] == "voltageKv" + + +def test_hint_tolerates_a_calibrated_pixel_size(): + # recorded 1.7675 at 2x2 binning -> 0.8837 per binning, against a stated 0.84 (5.2% apart) + src = {"epu": {"pixelSizeA": 1.7675, "binning": [2, 2], "superResolutionFactor": 2}, + "sheet": {"pixelSizeA": 0.84}} + out = compare_acquisition(src) + assert out[0].field == "pixelSizeA" and "perBinning" in out[0].message + assert "5.2%" in out[0].message