Add EPU XML decoding and a cross-source acquisition check

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

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`

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

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

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"</" + _NS + rb"Key>\s*<" + _NS +
+           rb"Value[^>]*>([^<]*)</")
+    m = re.search(b"<" + pat, data)
+    return m.group(1).decode("utf-8", "replace").strip() if m else None
+
+
+def _num(text: Optional[str]) -> 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"<MicroscopeImage" not in data:
+        return None
+
+    volts = _first_float(data, rb"<AccelerationVoltage>([^<]*)</AccelerationVoltage>")
+    px_m = _first_float(data, rb"<pixelSize>\s*<x>\s*<numericValue>([^<]*)</numericValue>")
+    mag = _first_float(data, rb"<NominalMagnification>([^<]*)</NominalMagnification>")
+    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"<ExposureTime>([^<]*)</ExposureTime>", data))
+                     if v and v > 0), None)
+
+    bin_m = re.search(rb"<Binning\b[^>]*>\s*<" + _NS + rb"x>(\d+)</" + _NS + rb"x>\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"<ReadoutArea\b[^>]*>\s*<" + _NS + rb"height>(\d+)</" + _NS + rb"height>\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"<EnergySelectionSlitWidth>([^<]*)</EnergySelectionSlitWidth>")
+    sw = re.search(rb"<ApplicationSoftware>([^<]*)</ApplicationSoftware>", data)
+    sw_v = re.search(rb"<ApplicationSoftwareVersion>([^<]*)</ApplicationSoftwareVersion>", 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