diff --git a/CHANGELOG.md b/CHANGELOG.md index 0231dd0..14efde2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,25 @@ All notable changes to this project are documented here. The format follows ## [Unreleased] ### Added +- **Uncapped MS1 from the RAW file.** `spectrum_cache.patch_ms1()` / `python -m + mzml_utils.spectrum_cache patch-ms1` replaces every MS1 scan's peaks in a spectrum cache with the + full centroid list read from the Thermo RAW (`mzml_utils.thermo_raw`, RawFileReader under Mono; + the DLLs are found, not shipped). It undoes a per-scan "N most intense" cap from the conversion + preset and refuses anything else: same MS1 scans and RTs, every cached peak a RAW peak, the + cached ones the RAW's most intense, re-capping the patched scan reproduces the cache bit for bit, + and MS2 rows byte-identical (digest before and after). `mz_mode="calibrated"` handles caches + built from MSFragger's `*_calibrated.mzML`: kept peaks stay, the removed ones get the scan's own + calibration curve (held-out error p99 under 0.2 ppm). Provenance goes into the cache meta + (`ms1_patch`: RAW path, size, sha256, cap undone, peaks before/after). +- `open_spectra(db, ms1_cap=N)` / `SpectrumCache(ms1_cap=N)` serve MS1 re-capped with + `cap_peaks()` (msconvert's rule, ties at the threshold excluded), so one tool can run on capped + and uncapped MS1 side by side. `ms1_cap_status()` / `warn_if_ms1_capped()` tell a tool when its + MS1 is capped and unpatched. +- `mzml_utils.ms1quant`: precursor MS1 intensity, IonQuant-like (sum of the top-3 isotope-cluster + scans within rttol; nearest peak per isotope within mztol; minscans, minisotopes), parameters + from `fragpipe.workflow`. Ported from NGlycoMM's `ms1_xic_extractor.py` and identical to it on + 9,215 of 9,215 precursors of its validation fraction. Reports top-3, apex, area, apex RT and, + on capped MS1, `apex_floor_ratio` / `cap_censored`. - `mzml_utils.chrom_extract`: `extract_chromatograms()` takes a peak table (rows with an m/z and RT window plus any annotation columns) and extracts one chromatogram per row and run in one pass per run, carrying every row column through; `chromatogram_rows()` flattens the diff --git a/README.md b/README.md index e079000..d08ed3a 100644 --- a/README.md +++ b/README.md @@ -34,6 +34,9 @@ pip install -e ".[dev]" | `fragments` | Fragment ion calculator (b/y/c/z/Y/oxonium), peak matching, false match rate | | `constants` | Physical constants, amino acid masses, common ion lists | | `utils` | Spectrum ID parsing, file finding, modification string parsing | +| `spectrum_cache` | Build-once local cache of an mzML; `open_spectra()`; `patch_ms1()` replaces capped MS1 with the RAW's full centroid list | +| `ms1quant` | Precursor MS1 intensity, IonQuant-like (top-3 isotope-cluster scans), with a cap-censoring flag | +| `thermo_raw` | Thermo RAW MS1 centroids through RawFileReader under Mono (DLLs found, not shipped) | ## Quick Start @@ -121,6 +124,27 @@ functions. `apex_rt=` measures the peak whose boundaries contain the MS2 time, a `nearest_valley=True` stops R's left boundary from running to the start of a zero-filled trace. These are measurements, not verdicts: no threshold is applied. +### Undo a per-scan peak cap and quantify precursors from MS1 + +A conversion preset that keeps only the N most intense peaks per scan truncates busy MS1 surveys. +If the Thermo RAW is at hand, patch the spectrum cache with every MS1 centroid from it (needs +`mono`, `mcs` and Thermo's RawFileReader DLLs; set `THERMO_PWIZ_DIR` if they are not found): + +```bash +python -m mzml_utils.spectrum_cache patch-ms1 --dry-run run.mzML # check only +python -m mzml_utils.spectrum_cache patch-ms1 run.mzML # RAW beside the mzML +python -m mzml_utils.spectrum_cache ms1-status run.mzML +``` + +```python +from mzml_utils import open_spectra, Ms1QuantParams, quantify_precursors + +params = Ms1QuantParams.from_workflow("fragpipe.workflow") # ionquant.* settings +rows = [{"ms2_scan": 12345, "charge": 3}] # anchor: the selected ion +uncapped = quantify_precursors(open_spectra("run.mzML"), rows, params) +capped = quantify_precursors(open_spectra("run.mzML", ms1_cap=1000), rows, params) # as searched +``` + ### Calculate fragment ions ```python diff --git a/src/mzml_utils/__init__.py b/src/mzml_utils/__init__.py index 3126017..28b4694 100644 --- a/src/mzml_utils/__init__.py +++ b/src/mzml_utils/__init__.py @@ -17,6 +17,9 @@ utils - Spectrum ID parsing, file finding, modification parsing structure - Protein structure databases (AlphaFold, PDBe, UniProt, InterPro) behind one resolver + ms1quant - Precursor MS1 intensity (IonQuant-like top-3), cap-censoring flag + thermo_raw - Thermo RAW MS1 centroids via RawFileReader under Mono, for + spectrum_cache.patch_ms1 (undoing a per-scan peak cap) ``structure`` is imported lazily: it pulls in ``requests``, and the spectrum side of this package has no use for HTTP. ``mzml_utils.structure`` and @@ -38,8 +41,17 @@ cache_path_for, build_cache, verify_cache, + patch_ms1, + ms1_cap_status, + warn_if_ms1_capped, + cap_peaks, ) +# ms1quant -- IonQuant-like precursor MS1 intensity (top-3 isotope-cluster scans), with a flag +# for results the per-scan peak cap may have censored. thermo_raw (Thermo RAW via Mono) is +# imported on use by spectrum_cache.patch_ms1. +from .ms1quant import Ms1QuantParams, Ms1Quant, Ms1Index, precursor_ms1, quantify_precursors + # ions from .ions import ( ppm_error, @@ -163,6 +175,8 @@ from .utils import parse_spectrum_id, find_mzml_file, parse_modifications __all__ = [ + "patch_ms1", "ms1_cap_status", "warn_if_ms1_capped", "cap_peaks", + "Ms1QuantParams", "Ms1Quant", "Ms1Index", "precursor_ms1", "quantify_precursors", # version "__version__", # reader diff --git a/src/mzml_utils/ms1quant.py b/src/mzml_utils/ms1quant.py new file mode 100644 index 0000000..101e922 --- /dev/null +++ b/src/mzml_utils/ms1quant.py @@ -0,0 +1,227 @@ +"""Precursor MS1 intensity, IonQuant-like, for every PSM of a run -- and whether the cap censored it. + +Ported from NGlycoMM ``data_analysis/tools/ms1_xic_extractor.py`` (2026-06-16), which was tuned +against IonQuant 1.11.20 ``Intensity`` on 4,887 precursors: Pearson(log) 0.924, Spearman 0.926, +median ratio 0.94, uniform across glycan classes. The rule, unchanged: + +* For every MS1 scan within ``rttol`` minutes of the MS2 scan's retention time, sum the isotope + cluster at the precursor m/z: isotopes ``iso_lo..iso_hi`` (relative to the anchor m/z) at + spacing 1.0033548378 / z, each taken as the *nearest* peak within ``mztol`` ppm (on an exact + tie, the higher-m/z peak). A scan counts when at least ``min_isotopes`` isotopes are found. +* A feature needs at least ``min_scans`` such scans; its intensity is the sum of the ``topn`` + largest scan sums ("top-3"). The apex scan sum, the trace area and the apex RT are reported too. + +Parameters come from the search: ``Ms1QuantParams.from_workflow`` reads ``ionquant.mztol``, +``rttol``, ``minscans`` and ``minisotopes`` from a FragPipe ``fragpipe.workflow``. + +Why not ``xic.extract_xics``: it makes one pass over the run per target and sums every peak in the +window. Here MS1 is indexed once per run and sliced by retention time, and each isotope takes the +nearest peak, the rule that was validated against IonQuant. + +**Capped MS1.** On an mzML capped at N peaks per scan, a weak isotope can be missing because it fell +below the cap, not because it is absent. ``apex_floor_ratio`` is the strongest isotope of the apex +scan over that scan's lowest kept peak (inf when the scan is not at the cap). ``cap_censored`` is +True when that ratio is at most ``censor_factor`` (3), i.e. the feature itself sits near the cap +floor, or, for a missing feature, when any MS1 scan in its window sits at the cap. Weak isotopes +of a strong feature routinely sit near the floor and do not set the flag: they change the value +little, while a feature near the floor can appear, vanish or change scans when uncapped. +On a cache patched from the RAW (``spectrum_cache.patch_ms1``) nothing is capped and the flag is +always False. Run the same estimator on ``open_spectra(db, ms1_cap=N)`` to see the capped result. +""" +from __future__ import annotations + +import configparser +import re +from dataclasses import dataclass, asdict +from typing import Dict, Iterable, List, Optional + +import numpy as np + +ISOTOPE_SPACING = 1.0033548378 # 13C - 12C + + +@dataclass(frozen=True) +class Ms1QuantParams: + mztol_ppm: float = 10.0 + rttol_min: float = 0.4 + min_scans: int = 3 + min_isotopes: int = 2 + iso_lo: int = -1 + iso_hi: int = 4 + topn: int = 3 + censor_factor: float = 3.0 + + @classmethod + def from_workflow(cls, path: str, **overrides) -> "Ms1QuantParams": + """IonQuant's mztol / rttol / minscans / minisotopes from a fragpipe.workflow.""" + text = open(path, encoding="utf-8", errors="replace").read() + cp = configparser.ConfigParser(interpolation=None, strict=False) + cp.read_string("[w]\n" + re.sub(r"(?m)^#.*$", "", text)) + w = cp["w"] + kw = {} + for key, name, typ in (("ionquant.mztol", "mztol_ppm", float), ("ionquant.rttol", "rttol_min", float), + ("ionquant.minscans", "min_scans", int), ("ionquant.minisotopes", "min_isotopes", int)): + if key in w and w[key].strip(): + kw[name] = typ(float(w[key])) + kw.update(overrides) + return cls(**kw) + + +@dataclass +class Ms1Quant: + top: float # sum of the topn largest scan sums (0 when no feature) + apex: float # largest scan sum + area: float # trapezoid over the window's scan sums (scans without a cluster count 0) + n_scans: int # scans with >= min_isotopes isotopes + apex_rt: float # RT of the apex scan (nan when no feature) + n_isotopes_apex: int # isotopes found in the apex scan + apex_floor_ratio: float # strongest apex isotope / that scan's lowest kept peak (inf if not capped) + cap_censored: bool # see the module docstring + + def as_dict(self) -> dict: + return asdict(self) + + +class Ms1Index: + """Every MS1 scan of a run in memory, sorted by RT, plus every scan's RT. + + ``reader`` is anything from ``open_spectra`` (a cache, possibly re-capped with ``ms1_cap``, + or an MzMLReader). The cap is read with ``spectrum_cache.ms1_cap_status``. + """ + + def __init__(self, reader): + from .spectrum_cache import ms1_cap_status + self.cap_status = ms1_cap_status(reader) + cap = self.cap_status["cap"] if self.cap_status["capped"] else None + self.scan_rt: Dict[int, float] = {} + conn = getattr(reader, "_conn", None) + if conn is not None: # a SpectrumCache: every RT from SQL, MS1 peaks decoded only + self.scan_rt = {int(s): float(rt) for s, rt in conn.execute("SELECT scan_num, rt FROM scans")} + spectra = reader.iter_spectra(ms_level=1) + else: + spectra = self._ms1_and_rts(reader) + rows = [] + for sp in spectra: + mz = np.asarray(sp.mz, dtype=np.float64) + it = np.asarray(sp.intensity, dtype=np.float64) + if len(mz) and np.any(np.diff(mz) < 0): + o = np.argsort(mz, kind="stable") + mz, it = mz[o], it[o] + at_cap = bool(cap) and len(mz) >= cap - 1 # N-1 when ties at the threshold were dropped + floor = float(it.min()) if (at_cap and len(it)) else 0.0 + rows.append((float(sp.rt), mz, it, floor)) + # RT in seconds (some converters) -> minutes, as the NGlycoMM extractor did. + rt_max = max([r[0] for r in rows] + list(self.scan_rt.values()) + [0.0]) + if rt_max > 600: + rows = [(r[0] / 60.0,) + r[1:] for r in rows] + self.scan_rt = {k: v / 60.0 for k, v in self.scan_rt.items()} + rows.sort(key=lambda r: r[0]) + self.rt = np.array([r[0] for r in rows], dtype=np.float64) + self.mz = [r[1] for r in rows] + self.intensity = [r[2] for r in rows] + self.floor = np.array([r[3] for r in rows], dtype=np.float64) + self.capped = cap is not None + + def _ms1_and_rts(self, reader): + for sp in reader.iter_spectra(): + self.scan_rt[int(sp.scan_num)] = float(sp.rt) + if sp.ms_level == 1: + yield sp + + def __len__(self): + return len(self.rt) + + +def _nearest_peak(mz: np.ndarray, target: float, tol_da: float) -> int: + """Index of the peak nearest ``target`` within ``tol_da`` (ties: the higher m/z), else -1.""" + k = int(np.searchsorted(mz, target)) + best, bestd = -1, tol_da + for kk in (k - 1, k): + if 0 <= kk < mz.shape[0]: + d = abs(float(mz[kk]) - target) + if d <= bestd: + bestd, best = d, kk + return best + + +def isotope_cluster(mz: np.ndarray, intensity: np.ndarray, anchor_mz: float, charge: int, + params: Ms1QuantParams = Ms1QuantParams()): + """(sum, n_found, strongest_found) of the isotope cluster in one scan.""" + s, n, strongest = 0.0, 0, 0.0 + for i in range(params.iso_lo, params.iso_hi + 1): + target = anchor_mz + i * ISOTOPE_SPACING / charge + j = _nearest_peak(mz, target, target * params.mztol_ppm * 1e-6) + if j >= 0: + v = float(intensity[j]) + s += v + n += 1 + strongest = max(strongest, v) + return s, n, strongest + + +def precursor_ms1(index: Ms1Index, anchor_mz: float, charge: int, rt_center: float, + params: Ms1QuantParams = Ms1QuantParams()) -> Ms1Quant: + """The feature at ``anchor_mz`` / ``charge`` around ``rt_center`` (minutes).""" + nan = float("nan") + if not charge or charge < 1 or not (anchor_mz > 0) or not (rt_center == rt_center): + return Ms1Quant(0.0, 0.0, 0.0, 0, nan, 0, nan, False) + i0 = int(np.searchsorted(index.rt, rt_center - params.rttol_min)) + i1 = int(np.searchsorted(index.rt, rt_center + params.rttol_min)) + sums, rts, kept = [], [], [] + window_at_cap = False + for idx in range(i0, i1): + s, n, strongest = isotope_cluster(index.mz[idx], index.intensity[idx], anchor_mz, charge, params) + window_at_cap |= index.floor[idx] > 0 + ok = n >= params.min_isotopes and s > 0 + rts.append(index.rt[idx]) + sums.append(s if ok else 0.0) + if ok: + kept.append((s, idx, n, strongest)) + if len(kept) < params.min_scans: + return Ms1Quant(0.0, 0.0, 0.0, len(kept), nan, 0, nan, bool(window_at_cap)) + vals = np.array([k[0] for k in kept]) + top = float(np.sort(vals)[-params.topn:].sum()) + a = int(np.argmax(vals)) + _, aidx, an, astrong = kept[a] + floor = index.floor[aidx] + ratio = float(astrong / floor) if floor > 0 else float("inf") + trapz = getattr(np, "trapezoid", None) or np.trapz + area = float(trapz(sums, rts)) if len(rts) > 1 else 0.0 + return Ms1Quant(top, float(vals[a]), area, len(kept), float(index.rt[aidx]), int(an), ratio, + bool(ratio <= params.censor_factor)) + + +def quantify_precursors(source, rows: Iterable[dict], params: Ms1QuantParams = Ms1QuantParams(), + ms1_cap: Optional[int] = None) -> List[dict]: + """MS1 for many precursors of one run. + + ``source`` is a path (opened with ``open_spectra``, honouring ``ms1_cap``) or an open reader. + Each row needs ``ms2_scan`` and ``charge``; ``anchor_mz`` defaults to the MS2 scan's selected + ion m/z from the reader (instrument coordinates) and ``rt_center`` to the MS2 scan's RT. + Returns one dict per row: the row's keys plus the ``Ms1Quant`` fields and ``anchor_source``. + """ + from .spectrum_cache import open_spectra + own = not hasattr(source, "iter_spectra") + reader = open_spectra(str(source), **({"ms1_cap": ms1_cap} if ms1_cap else {})) if own else source + try: + index = Ms1Index(reader) + out = [] + for r in rows: + r = dict(r) + scan = int(r["ms2_scan"]) + anchor = r.get("anchor_mz") + src = "given" + if anchor is None or not (anchor == anchor): + sp = reader.get_spectrum(scan) + anchor = float(sp.precursor_mz) if sp is not None else float("nan") + src = "selected_ion" + rtc = r.get("rt_center") + if rtc is None or not (rtc == rtc): + rtc = index.scan_rt.get(scan, float("nan")) + q = precursor_ms1(index, float(anchor), int(r["charge"]), float(rtc), params) + r.update(q.as_dict(), anchor_mz=float(anchor), anchor_source=src, rt_center=float(rtc)) + out.append(r) + return out + finally: + if own: + reader.close() diff --git a/src/mzml_utils/spectrum_cache.py b/src/mzml_utils/spectrum_cache.py index be32488..b9f8bb9 100644 --- a/src/mzml_utils/spectrum_cache.py +++ b/src/mzml_utils/spectrum_cache.py @@ -482,7 +482,11 @@ class SpectrumCache: """ def __init__(self, db_path: str, mirror_local: bool = True, in_memory: bool = False, - local_dir: Optional[str] = None): + local_dir: Optional[str] = None, ms1_cap: Optional[int] = None): + """``ms1_cap``: serve MS1 scans re-capped to their N most intense peaks (``cap_peaks``). + On a cache patched from the RAW this reproduces the capped mzML exactly, so the same + tool can be run on capped and uncapped MS1 side by side.""" + self.ms1_cap = ms1_cap self.source_path = db_path path = db_path if mirror_local and self._is_network(db_path): @@ -513,9 +517,14 @@ def _mirror(src: str, local_dir: str) -> str: return dst # -- MzMLReader-compatible API ---------------------------------------- # + def _view(self, sp: Spectrum) -> Spectrum: + if self.ms1_cap and sp.ms_level == 1 and len(sp.mz) > self.ms1_cap: + sp.mz, sp.intensity = cap_peaks(sp.mz, sp.intensity, self.ms1_cap) + return sp + def get_spectrum(self, scan_num: int) -> Optional[Spectrum]: row = self._conn.execute("SELECT * FROM scans WHERE scan_num=?", (scan_num,)).fetchone() - return _spectrum_from_row(row) if row else None + return self._view(_spectrum_from_row(row)) if row else None def iter_spectra(self, ms_level: Optional[int] = None, rt_range: Optional[Tuple[float, float]] = None) -> Iterator[Spectrum]: @@ -537,7 +546,7 @@ def iter_spectra(self, ms_level: Optional[int] = None, if where: sql += " WHERE " + " AND ".join(where) for row in self._conn.execute(sql + " ORDER BY scan_num", params): - yield _spectrum_from_row(row) + yield self._view(_spectrum_from_row(row)) def find_best_ms1(self, ms2_scan: int, precursor_mz: float, n_ms1: int = 3, tol_da: float = 0.02) -> Optional[int]: @@ -657,7 +666,9 @@ def open_spectra(path: str, prefer_cache: bool = True, in_memory=None, """Drop-in for MzMLReader(path). Given a .db/.sqlite/.cache path, opens the cache. Given an mzML path and prefer_cache=True, auto-uses the co-located cache (/spectra_cache/.spectra.db) when it exists, else MzMLReader. - Network caches load into RAM by default so nothing is left on local disk.""" + Network caches load into RAM by default so nothing is left on local disk. + ``ms1_cap=N`` (caches only) serves MS1 re-capped to its N most intense peaks (``cap_peaks``), + which on a cache patched from the RAW reproduces the capped mzML exactly.""" p = str(path) if p.lower().endswith((".db", ".sqlite", ".cache")): return _open_cache(p, in_memory, mirror_local, **kw) @@ -669,6 +680,393 @@ def open_spectra(path: str, prefer_cache: bool = True, in_memory=None, return MzMLReader(p) +# --------------------------------------------------------------------------- # +# MS1 from the RAW file: undo a per-scan peak cap (e.g. msconvert "1000 most intense") +# --------------------------------------------------------------------------- # +class PatchError(RuntimeError): + """The RAW export does not match the cache; nothing was written.""" + + +def _ms2_digest(conn: sqlite3.Connection) -> str: + """sha256 over every non-MS1 row's scan number and peak blobs, in scan order.""" + import hashlib + h = hashlib.sha256() + for r in conn.execute("SELECT scan_num, mz_blob, int_blob FROM scans WHERE ms_level != 1 ORDER BY scan_num"): + h.update(str(r[0]).encode()) + h.update(r[1] or b"") + h.update(r[2] or b"") + return h.hexdigest() + + +def cap_peaks(mz, intensity, n: int): + """The ``n`` most intense peaks, excluding ties at the threshold (msconvert's rule), in m/z order. + + msconvert's "threshold count N most-intense" keeps the N most intense peaks and drops every + peak tied at the threshold when keeping them all would exceed N, so a scan can end with N - 1. + Applying this to a patched scan returns the capped scan exactly; ``patch_ms1`` checks that. + """ + mz = np.asarray(mz) + it = np.asarray(intensity) + if n is None or len(it) <= n: + return mz, it + desc = np.sort(it)[::-1] + thr = desc[n - 1] + keep = it > thr if desc[n] == thr else it >= thr + return mz[keep], it[keep] + + +def check_ms1_against_raw(cached_mz, cached_int, raw_mz, raw_int, mz_ppm: float = 0.5) -> Optional[str]: + """None when the cached peaks are the RAW peaks after a most-intense cap; else the reason. + + The cached peaks must each match a distinct RAW peak (m/z within ``mz_ppm``, intensity equal + at float32), and when the RAW scan has more peaks, the cached ones must be its most intense: + no RAW peak left out may be more intense than the weakest one kept. + """ + cmz = np.asarray(cached_mz, dtype=np.float64) + cit = np.asarray(cached_int, dtype=np.float32) + rmz = np.asarray(raw_mz, dtype=np.float64) + rit = np.asarray(raw_int, dtype=np.float64) + if len(cmz) == 0: + return None if len(rmz) == 0 else f"cache scan is empty, the RAW scan has {len(rmz)} peaks" + if len(rmz) == 0: + return f"{len(cmz)} cached peaks, none in the RAW scan" + if len(cmz) > len(rmz): + return f"cache has more peaks ({len(cmz)}) than the RAW scan ({len(rmz)})" + idx = np.clip(np.searchsorted(rmz, cmz), 1, len(rmz) - 1) if len(rmz) > 1 else np.zeros(len(cmz), int) + if len(rmz) > 1: + left, right = rmz[idx - 1], rmz[idx] + j = np.where(np.abs(cmz - left) <= np.abs(cmz - right), idx - 1, idx) + else: + j = idx + ppm = np.abs(rmz[j] - cmz) / cmz * 1e6 + if (ppm > mz_ppm).any(): + k = int(np.argmax(ppm)) + return f"cached peak m/z {cmz[k]:.5f} has no RAW peak within {mz_ppm} ppm (nearest {ppm[k]:.2f} ppm)" + if len(np.unique(j)) != len(j): + return "two cached peaks map to one RAW peak" + rel = np.abs(rit[j].astype(np.float32) - cit) / np.maximum(np.abs(cit), 1.0) + if (rel > 1e-6).any(): + k = int(np.argmax(rel)) + return f"intensity differs at m/z {cmz[k]:.5f}: cache {cit[k]:.6g}, RAW {rit[j][k]:.6g}" + if len(j) < len(rmz): + left_out = np.ones(len(rmz), bool) + left_out[j] = False + if rit[left_out].max() > float(cit.min()) * (1 + 1e-6): + return (f"not a most-intense cap: a RAW peak left out ({rit[left_out].max():.6g}) is more " + f"intense than the weakest kept ({float(cit.min()):.6g})") + return None + + +def merge_calibrated_ms1(cached_mz, cached_int, raw_mz, raw_int, max_shift_ppm: float = 50.0, + max_heldout_ppm: float = 1.0): + """Uncap a scan whose kept peaks were m/z-recalibrated (e.g. MSFragger's *_calibrated.mzML). + + Each kept peak is paired with the RAW peak of identical intensity (float32) nearest in m/z, + within ``max_shift_ppm``; every kept peak must pair, uniquely, and the kept peaks must be the + RAW scan's most intense. The calibration shift (ppm) is then a function of m/z inside the scan + and is interpolated from the pairs onto the peaks the cap removed (constant beyond the ends). + Kept peaks keep their calibrated values exactly; added peaks get RAW intensity and + interpolated calibrated m/z. Accuracy is measured on the scan itself: shifts at odd-indexed + pairs are predicted from the even-indexed ones, and the 99th percentile must stay under + ``max_heldout_ppm``. + + Returns (mz, intensity, info) or raises ValueError with the reason. + """ + cmz = np.asarray(cached_mz, dtype=np.float64) + cit = np.asarray(cached_int, dtype=np.float32) + rmz = np.asarray(raw_mz, dtype=np.float64) + rit = np.asarray(raw_int, dtype=np.float64) + info = dict(n_kept=len(cmz), n_added=0, max_abs_shift_ppm=0.0, heldout_p99_ppm=0.0, heldout_max_ppm=0.0) + if len(cmz) == 0: + if len(rmz): + raise ValueError(f"cache scan is empty, the RAW scan has {len(rmz)} peaks") + return cmz, cit.astype(np.float64), info + r32 = rit.astype(np.float32) + by_int: Dict[float, List[int]] = {} + for k, v in enumerate(r32): + by_int.setdefault(float(v), []).append(k) + used = set() + pair = np.empty(len(cmz), dtype=int) + for k, (m, v) in enumerate(zip(cmz, cit)): + cand = [q for q in by_int.get(float(v), ()) if q not in used and abs(rmz[q] - m) / m * 1e6 <= max_shift_ppm] + if not cand: + raise ValueError(f"kept peak m/z {m:.5f} (intensity {v:.6g}) has no RAW peak of that intensity " + f"within {max_shift_ppm} ppm") + q = min(cand, key=lambda q: abs(rmz[q] - m)) + used.add(q) + pair[k] = q + if len(pair) < len(rmz): + left_out = np.ones(len(rmz), bool) + left_out[pair] = False + if rit[left_out].max() > float(cit.min()) * (1 + 1e-6): + raise ValueError("not a most-intense cap: a RAW peak left out is more intense than the weakest kept") + o = np.argsort(rmz[pair]) + a, b = rmz[pair][o], cmz[o] + if (np.diff(b) < 0).any(): + raise ValueError("recalibration reorders peaks; not a smooth m/z function") + shift = (b - a) / a * 1e6 + info["max_abs_shift_ppm"] = float(np.abs(shift).max()) + if len(a) >= 10: + tr_i, te_i = np.arange(0, len(a), 2), np.arange(1, len(a), 2) + res = np.abs(np.interp(a[te_i], a[tr_i], shift[tr_i]) - shift[te_i]) + info["heldout_p99_ppm"] = float(np.percentile(res, 99)) + info["heldout_max_ppm"] = float(res.max()) + if info["heldout_p99_ppm"] > max_heldout_ppm: + raise ValueError(f"calibration shift is not smooth in m/z (held-out p99 {info['heldout_p99_ppm']:.3f} ppm)") + if len(pair) == len(rmz): + return cmz, cit.astype(np.float64), info + add = np.setdiff1d(np.arange(len(rmz)), pair) + add_mz = rmz[add] * (1 + np.interp(rmz[add], a, shift) / 1e6) + mz = np.concatenate([cmz, add_mz]) + it = np.concatenate([cit.astype(np.float64), rit[add]]) + srt = np.argsort(mz, kind="stable") + info["n_added"] = int(len(add)) + return mz[srt], it[srt], info + + +def _default_raw_for(mzml_or_db: str) -> Optional[str]: + """/.raw for an mzML (also stripping _calibrated / _uncalibrated).""" + p = os.path.abspath(mzml_or_db) + if p.lower().endswith(".db"): + d = os.path.dirname(os.path.dirname(p)) + stem = os.path.basename(p)[: -len(".spectra.db")] if p.endswith(".spectra.db") else os.path.splitext(os.path.basename(p))[0] + else: + d, stem = os.path.dirname(p), os.path.splitext(os.path.basename(p))[0] + for s in (stem, re.sub(r"_(un)?calibrated$", "", stem)): + for ext in (".raw", ".RAW"): + cand = os.path.join(d, s + ext) + if os.path.isfile(cand): + return cand + return None + + +def patch_ms1(cache: str, raw_path: Optional[str] = None, *, export: Optional[dict] = None, + pwiz_dir: Optional[str] = None, dry_run: bool = False, force: bool = False, + mz_ppm: float = 0.5, rt_tol_min: float = 1e-4, mz_mode: str = "auto", + log=None) -> dict: + """Replace every MS1 scan's peaks in a spectrum cache with the full centroid list from the RAW. + + ``cache`` is the .spectra.db or the mzML it sits beside. ``raw_path`` defaults to the RAW next + to the mzML. ``export`` ({scan: thermo_raw.RawMs1Scan}) bypasses the RAW read (tests). + + Checked before anything is written; any failure raises PatchError and leaves the cache as it + was: + * the cache and the RAW hold the same MS1 scan numbers, at the same retention times; + * every cached MS1 peak is a RAW peak (m/z within ``mz_ppm``, same intensity), and where the + RAW scan has more peaks the cached ones are its most intense, i.e. the only difference is + a most-intense cap. A different centroider would fail here, so the patch never mixes two; + * every RAW MS1 scan is centroid data (exporter source 1 or 2), not profile. + ``mz_mode``: "exact" requires the cached m/z to be the RAW m/z and writes the RAW peaks. + "calibrated" is for a cache built from a recalibrated mzML (MSFragger's *_calibrated.mzML): + kept peaks stay as they are and the removed ones are added with m/z recalibrated by the + scan's own shift curve (``merge_calibrated_ms1``). "auto" picks "calibrated" when the cache's + source mzML name ends in _calibrated, else "exact". + MS2 and higher rows are not touched: a digest of their peak blobs is compared before and after. + TIC and base peak stay as recorded (the mzML takes them from the scan header). + + The cache records ``ms1_patch`` in its meta: RAW path, size and sha256, exporter version, the + cap that was undone, peaks before and after, and the MS2 digest. A cache already patched from + the same RAW is left alone; from a different RAW it raises unless ``force``. + """ + from . import thermo_raw + + def _log(m): + if log: + log(m) + + db = cache if cache.lower().endswith((".db", ".sqlite", ".cache")) else cache_path_for(cache) + if not os.path.isfile(db): + raise FileNotFoundError(f"no spectrum cache at {db}; build it first") + if export is None: + raw_path = raw_path or _default_raw_for(cache) + if not raw_path: + raise FileNotFoundError(f"no RAW file found beside {cache}; pass raw_path / --raw") + t0 = time.time() + + # Work on a local copy: SQLite writes over SMB are unsafe, and the original stays intact + # until the patched copy is complete and checked. + work_dir = tempfile.mkdtemp(prefix="ms1patch_", dir=_build_tmpdir()) + work = os.path.join(work_dir, os.path.basename(db)) + try: + shutil.copy2(db, work) + conn = sqlite3.connect(work) + conn.row_factory = sqlite3.Row + meta = {r["key"]: json.loads(r["value"]) for r in conn.execute("SELECT key, value FROM meta")} + if mz_mode == "auto": + src = os.path.splitext(os.path.basename(str(meta.get("source_mzml") or db)))[0] + src = src[: -len(".spectra")] if src.endswith(".spectra") else src + mz_mode = "calibrated" if src.endswith("_calibrated") else "exact" + if mz_mode not in ("exact", "calibrated"): + raise ValueError(f"mz_mode must be auto, exact or calibrated, not {mz_mode!r}") + raw_id = None + if raw_path: + st = os.stat(raw_path) + raw_id = dict(raw_path=os.path.abspath(raw_path), raw_size=st.st_size) + prev = meta.get("ms1_patch") + if prev: + same = raw_id and prev.get("raw_size") == raw_id["raw_size"] and \ + os.path.basename(prev.get("raw_path", "")) == os.path.basename(raw_id["raw_path"]) + if same and not force: + conn.close() + _log(f"already patched from {prev.get('raw_path')}; nothing to do") + return dict(prev, status="already_patched", db_path=db) + if not force: + conn.close() + raise PatchError(f"{db} was patched from {prev.get('raw_path')}; pass force=True to redo it") + + if export is None: + _log(f"exporting MS1 from {raw_path}") + export = thermo_raw.export_ms1(raw_path, pwiz_dir=pwiz_dir) + rows = conn.execute("SELECT scan_num, rt, n_peaks, mz_blob, int_blob FROM scans " + "WHERE ms_level = 1 ORDER BY scan_num").fetchall() + cache_scans = [r["scan_num"] for r in rows] + if set(cache_scans) != set(export): + only_c = sorted(set(cache_scans) - set(export))[:5] + only_r = sorted(set(export) - set(cache_scans))[:5] + raise PatchError(f"MS1 scan numbers differ: cache-only {only_c}, RAW-only {only_r}") + profile = [s for s, x in export.items() if x.source not in (1, 2)] + if profile: + raise PatchError(f"{len(profile)} RAW MS1 scans have no centroids (e.g. scan {profile[0]})") + + before_digest = _ms2_digest(conn) + counts: Dict[int, int] = {} + for r in rows: + counts[r["n_peaks"]] = counts.get(r["n_peaks"], 0) + 1 + cap_hint = max(counts) if counts else None + n_before = n_after = 0 + truncated, cap_values = 0, {} + updates = [] + heldout = [] + max_shift = 0.0 + for r in rows: + s = r["scan_num"] + x = export[s] + if r["rt"] is not None and abs(float(r["rt"]) - x.rt) > rt_tol_min: + raise PatchError(f"scan {s}: retention time {r['rt']} in the cache, {x.rt} in the RAW") + cmz, cit = decode_peaks(r["mz_blob"], r["int_blob"]) + if mz_mode == "exact": + why = check_ms1_against_raw(cmz, cit, x.mz, x.intensity, mz_ppm=mz_ppm) + if why: + raise PatchError(f"scan {s}: {why}" + ( + " (if the cache was built from a recalibrated mzML, use mz_mode='calibrated')" + if "no RAW peak" in why else "")) + new_mz, new_it = x.mz, x.intensity + else: + try: + new_mz, new_it, info = merge_calibrated_ms1(cmz, cit, x.mz, x.intensity) + except ValueError as e: + raise PatchError(f"scan {s}: {e}") + heldout.append(info["heldout_p99_ppm"]) + max_shift = max(max_shift, info["max_abs_shift_ppm"]) + if len(new_mz) > len(cmz): + # Re-capping the patched scan must give back the cached scan, bit for bit at float32. + rmz_c, rit_c = cap_peaks(np.asarray(new_mz, np.float32), np.asarray(new_it, np.float32), cap_hint) + if len(rmz_c) != len(cmz) or not (np.array_equal(rmz_c, cmz.astype(np.float32)) + and np.array_equal(rit_c, cit.astype(np.float32))): + raise PatchError(f"scan {s}: re-capping the patched scan at {cap_hint} does not " + f"reproduce the cached scan ({len(rmz_c)} vs {len(cmz)} peaks)") + n_before += len(cmz) + n_after += len(new_mz) + if len(x.mz) > len(cmz): + truncated += 1 + cap_values[len(cmz)] = cap_values.get(len(cmz), 0) + 1 + mzb, itb = encode_peaks(new_mz, new_it) + updates.append((mzb, itb, int(len(new_mz)), s)) + cap = max(cap_values, key=cap_values.get) if cap_values else None + report = dict( + db_path=db, raw_path=os.path.abspath(raw_path) if raw_path else None, + raw_size=raw_id["raw_size"] if raw_id else None, + exporter_version=thermo_raw.EXPORTER_VERSION, n_ms1=len(rows), + ms1_truncated_before=truncated, cap_undone=cap, recap_reproduces_cache=True, + truncated_peak_counts={str(k): v for k, v in sorted(cap_values.items())}, + ms1_peaks_before=n_before, ms1_peaks_after=n_after, ms2_digest=before_digest, + mz_ppm=mz_ppm, rt_tol_min=rt_tol_min, mz_mode=mz_mode, + ) + if mz_mode == "calibrated": + report["max_abs_calibration_shift_ppm"] = round(max_shift, 3) + report["heldout_p99_ppm_worst_scan"] = round(max(heldout), 4) if heldout else 0.0 + report["heldout_p99_ppm_median_scan"] = round(float(np.median(heldout)), 4) if heldout else 0.0 + if dry_run: + conn.close() + report["status"] = "checked (dry run, nothing written)" + return report + + if raw_path: + _log("hashing the RAW file") + report["raw_sha256"] = thermo_raw.sha256_file(raw_path) + conn.executemany("UPDATE scans SET mz_blob = ?, int_blob = ?, n_peaks = ? WHERE scan_num = ?", updates) + if _ms2_digest(conn) != before_digest: + raise PatchError("MS2 rows changed during the patch") + report["patched_at"] = time.strftime("%Y-%m-%dT%H:%M:%S") + report["seconds"] = round(time.time() - t0, 1) + set_meta(conn, ms1_patch={k: v for k, v in report.items() if k != "db_path"}) + conn.commit() + conn.execute("VACUUM") + conn.close() + # Publish: copy beside the original, then an atomic rename over it. + staged = db + ".ms1patch.tmp" + shutil.copy2(work, staged) + if os.path.getsize(staged) != os.path.getsize(work): + os.remove(staged) + raise PatchError("staged copy size mismatch; original left in place") + os.replace(staged, db) + report["status"] = "patched" + report["db_bytes"] = os.path.getsize(db) + _log(f"patched {len(rows)} MS1 scans ({truncated} were truncated at {cap}); " + f"{n_before} -> {n_after} peaks") + return report + finally: + shutil.rmtree(work_dir, ignore_errors=True) + + +def ms1_cap_status(reader) -> dict: + """Whether a reader's MS1 peaks look capped (a per-scan most-intense limit) or were patched. + + For a SpectrumCache: patched when its meta records ``ms1_patch``; otherwise capped when at + least 5% of MS1 scans share exactly the maximum peak count (a natural peak-count distribution + has no such pile-up). For an MzMLReader: read from the mzML header's dataProcessing + ("most intense count"). Returns {capped, cap, patched, frac_at_cap, source}. + """ + if isinstance(reader, SpectrumCache): + m = reader.meta() + if reader.ms1_cap: + return dict(capped=True, cap=reader.ms1_cap, patched=bool(m.get("ms1_patch")), + frac_at_cap=float("nan"), source="re-capped view") + if m.get("ms1_patch"): + return dict(capped=False, cap=m["ms1_patch"].get("cap_undone"), patched=True, + frac_at_cap=0.0, source="cache-meta") + tot = reader._conn.execute("SELECT COUNT(*) FROM scans WHERE ms_level = 1").fetchone()[0] + top = reader._conn.execute("SELECT n_peaks, COUNT(*) FROM scans WHERE ms_level = 1 " + "GROUP BY n_peaks ORDER BY n_peaks DESC LIMIT 1").fetchone() + if not tot or not top: + return dict(capped=False, cap=None, patched=False, frac_at_cap=0.0, source="peak-counts") + frac = top[1] / tot + capped = frac >= 0.05 and top[0] >= 50 + return dict(capped=capped, cap=int(top[0]) if capped else None, patched=False, + frac_at_cap=round(frac, 4), source="peak-counts") + path = getattr(reader, "path", None) or getattr(reader, "filepath", None) or getattr(reader, "filename", None) + cap = None + if path and os.path.isfile(str(path)): + with open(str(path), "rb") as fh: + head = fh.read(200_000).decode("utf-8", "replace") + m = re.search(r'most intense count[^"]*"\s+value="(\d+)"', head) + if m: + cap = int(m.group(1)) + return dict(capped=cap is not None, cap=cap, patched=False, frac_at_cap=float("nan"), + source="mzml-header" if path else "unknown") + + +def warn_if_ms1_capped(reader, what: str = "MS1 peaks") -> dict: + """Warn once per call site when MS1 is capped and not patched; returns ms1_cap_status.""" + import warnings + st = ms1_cap_status(reader) + if st["capped"]: + warnings.warn(f"{what}: MS1 scans are capped at {st['cap']} peaks and the cache is not " + f"patched from the RAW (spectrum_cache patch-ms1); absent MS1 peaks may be " + f"below the cap, not missing.", stacklevel=2) + return st + + # --------------------------------------------------------------------------- # # verify (run against the real mzML when the user chooses) # --------------------------------------------------------------------------- # @@ -740,7 +1138,45 @@ def _cli(): v.add_argument("--allow-network", action="store_true", help="allow verifying against a /Volumes source (can hang on SMB)") + pm = sub.add_parser("patch-ms1", + help="replace capped MS1 peaks in a cache with the full centroids from the RAW") + pm.add_argument("target", nargs="+", help="spectrum cache .db or the mzML it sits beside") + pm.add_argument("--raw", help="RAW file (default: /.raw); only with one target") + pm.add_argument("--pwiz-dir", help="directory with the Thermo RawFileReader DLLs") + pm.add_argument("--dry-run", action="store_true", help="check only, write nothing") + pm.add_argument("--force", action="store_true", help="redo a cache patched from another RAW") + pm.add_argument("--mz-ppm", type=float, default=0.5) + pm.add_argument("--mz-mode", choices=["auto", "exact", "calibrated"], default="auto", + help="calibrated: the cache was built from a recalibrated mzML (default auto: by name)") + + ms = sub.add_parser("ms1-status", help="is MS1 capped, and was the cache patched?") + ms.add_argument("target", nargs="+") + args = ap.parse_args() + if args.cmd == "patch-ms1": + if args.raw and len(args.target) > 1: + ap.error("--raw takes one target") + bad = 0 + for tgt in args.target: + try: + res = patch_ms1(tgt, args.raw, pwiz_dir=args.pwiz_dir, dry_run=args.dry_run, + force=args.force, mz_ppm=args.mz_ppm, mz_mode=args.mz_mode, + log=lambda m: print(f" {m}", flush=True)) + print(f"[{res['status']}] {tgt}: {res.get('n_ms1')} MS1, truncated " + f"{res.get('ms1_truncated_before')} at {res.get('cap_undone')}, peaks " + f"{res.get('ms1_peaks_before')} -> {res.get('ms1_peaks_after')}", flush=True) + except (PatchError, FileNotFoundError, RuntimeError) as e: + bad += 1 + print(f"[FAILED] {tgt}: {e}", flush=True) + sys.exit(1 if bad else 0) + if args.cmd == "ms1-status": + for tgt in args.target: + r = open_spectra(tgt) + try: + print(tgt, json.dumps(ms1_cap_status(r))) + finally: + r.close() + return if args.cmd == "build": if args.beside: out = cache_path_for(args.mzml) diff --git a/src/mzml_utils/thermo_raw.py b/src/mzml_utils/thermo_raw.py new file mode 100644 index 0000000..e8bf1d5 --- /dev/null +++ b/src/mzml_utils/thermo_raw.py @@ -0,0 +1,271 @@ +"""Thermo RAW access through Thermo's RawFileReader, run under Mono. + +Only what the spectrum cache needs: every MS1 scan's centroids, read from the RAW file itself. + +Why it exists. A conversion preset that keeps only the N most intense peaks per scan (msconvert's +``threshold count N most-intense``) truncates every busy MS1 survey, while MS2 scans rarely reach +the cap. Everything that reads MS1 peaks from such an mzML -- precursor isotope checks, +co-isolation, XICs, MS1 quantification -- then reads censored data. The RAW file still holds every +centroid, so the fix is to read MS1 from it and patch the spectrum cache +(``spectrum_cache.patch_ms1``) instead of re-converting. + +The centroids are the instrument's own label peaks, *including* reference and exception peaks +(``GetCentroidStream(scan, true)``). That is what ProteoWizard's "Thermo/Xcalibur peak picking" +writes: on an Orbitrap Eclipse run, every capped mzML MS1 peak was found in this export at 0 ppm +and identical intensity, and every capped scan held exactly the export's N most intense peaks. +Leaving reference/exception peaks out drops about 4% of them, mostly shoulders beside intense ions. + +The Thermo DLLs are not redistributed here. ``find_pwiz_dir`` looks for +``ThermoFisher.CommonCore.{Data,RawFileReader}.dll`` in, in order: an explicit directory, +``$THERMO_PWIZ_DIR`` (one directory, or several separated by ``os.pathsep``), a local copy made on +first use (``~/.cache/mzml_utils/thermo/dlls``), and a ThermoRawFileParser install under +``~/.local/opt/ThermoRawFileParser`` (it bundles the same two DLLs). A ProteoWizard install has +them too. The first hit off the local copy is copied into it, so later runs never read the DLLs +from a network share. ``mono`` and ``mcs`` (or ``csc``) must be on PATH. +""" +from __future__ import annotations + +import hashlib +import os +import shutil +import struct +import subprocess +import tempfile +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, Iterable, Optional + +import numpy as np + +DLLS = ("ThermoFisher.CommonCore.Data.dll", "ThermoFisher.CommonCore.RawFileReader.dll") + +#: Local installs that bundle the DLLs, searched last. +LOCAL_DLL_DIRS = (str(Path.home() / ".local" / "opt" / "ThermoRawFileParser"),) + + +def _cache_root() -> Path: + return Path(os.environ.get("MZML_UTILS_CACHE", Path.home() / ".cache" / "mzml_utils")) + + +def _has_dlls(d: Path) -> bool: + return all((d / name).is_file() for name in DLLS) + + +def find_pwiz_dir(user_value: Optional[str] = None, required: bool = True) -> Optional[Path]: + """Directory holding the two RawFileReader DLLs, or None (``required=False``). + + Search order: ``user_value``, ``$THERMO_PWIZ_DIR``, the local copy, ``LOCAL_DLL_DIRS``. A hit + outside the local copy is copied into it so the next call does not touch the network. + """ + local = _cache_root() / "thermo" / "dlls" + env = [d for d in os.environ.get("THERMO_PWIZ_DIR", "").split(os.pathsep) if d] + candidates = [user_value, *env, str(local), *LOCAL_DLL_DIRS] + for c in candidates: + if not c: + continue + d = Path(c) + try: + if not _has_dlls(d): + continue + except OSError: # an unreachable share + continue + if d != local: + try: + local.mkdir(parents=True, exist_ok=True) + for name in DLLS: + if not (local / name).is_file() or (local / name).stat().st_size != (d / name).stat().st_size: + shutil.copy2(d / name, local / name) + except OSError: + return d + return local if _has_dlls(local) else d + if required: + raise FileNotFoundError( + "Thermo RawFileReader DLLs not found (ThermoFisher.CommonCore.Data.dll, " + "ThermoFisher.CommonCore.RawFileReader.dll). Pass pwiz_dir= / --pwiz-dir or set " + "THERMO_PWIZ_DIR.") + return None + + +def raw_reader_available(pwiz_dir: Optional[str] = None) -> bool: + """True when mono, a C# compiler and the DLLs are all present.""" + return bool(shutil.which("mono") and (shutil.which("mcs") or shutil.which("csc")) + and find_pwiz_dir(pwiz_dir, required=False)) + + +# The exporter. Binary output, little-endian: +# b"MS1X", int32 version +# per MS1 scan: int32 scan, float64 rt_min, int32 source, int32 n, float64[n] mz, float64[n] intensity +# int32 -1 +# source: 1 = centroid stream (label peaks incl. reference/exception), 2 = centroid scan data, +# 3 = profile data (no centroid stream; NOT what msconvert's vendor peak picking writes). +_MS1_EXPORT_CS = r''' +using System; +using System.IO; +using System.Reflection; +using System.Runtime.CompilerServices; +using ThermoFisher.CommonCore.Data.Business; +using ThermoFisher.CommonCore.Data.FilterEnums; +using ThermoFisher.CommonCore.Data.Interfaces; +using ThermoFisher.CommonCore.RawFileReader; + +public static class Ms1Export +{ + public static int Main(string[] args) + { + if (args.Length < 3) { Console.Error.WriteLine("usage: Ms1Export "); return 2; } + var dllDir = args[0]; + AppDomain.CurrentDomain.AssemblyResolve += (s, e) => + { + var p = Path.Combine(dllDir, new AssemblyName(e.Name).Name + ".dll"); + return File.Exists(p) ? Assembly.LoadFrom(p) : null; + }; + return Run(args[1], args[2]); + } + + [MethodImpl(MethodImplOptions.NoInlining)] + static int Run(string raw, string outPath) + { + using (var rf = RawFileReaderAdapter.FileFactory(raw)) + { + if (rf == null || !rf.IsOpen || rf.IsError) { Console.Error.WriteLine("cannot open " + raw); return 3; } + rf.SelectInstrument(Device.MS, 1); + int first = rf.RunHeaderEx.FirstSpectrum, last = rf.RunHeaderEx.LastSpectrum; + int n = 0; + using (var bw = new BinaryWriter(File.Create(outPath))) + { + bw.Write(new byte[] { (byte)'M', (byte)'S', (byte)'1', (byte)'X' }); + bw.Write(1); + for (int s = first; s <= last; s++) + { + var f = rf.GetFilterForScanNumber(s); + if (f.MSOrder != MSOrderType.Ms) continue; + double rt = rf.RetentionTimeFromScanNumber(s); + double[] mz, it; int src; + var cs = rf.GetCentroidStream(s, true); + if (cs != null && cs.Length > 0) { mz = cs.Masses; it = cs.Intensities; src = 1; } + else + { + var st = rf.GetScanStatsForScanNumber(s); + var seg = rf.GetSegmentedScanFromScanNumber(s, st); + mz = seg.Positions; it = seg.Intensities; src = st.IsCentroidScan ? 2 : 3; + } + bw.Write(s); bw.Write(rt); bw.Write(src); bw.Write(mz.Length); + foreach (var v in mz) bw.Write(v); + foreach (var v in it) bw.Write(v); + n++; + } + bw.Write(-1); + } + Console.WriteLine("ms1_scans\t" + n + "\tfirst\t" + first + "\tlast\t" + last); + } + return 0; + } +} +''' + +EXPORTER_VERSION = hashlib.sha256(_MS1_EXPORT_CS.encode()).hexdigest()[:12] + + +@dataclass +class RawMs1Scan: + scan: int + rt: float # minutes + source: int # see the exporter comment above + mz: np.ndarray # float64 + intensity: np.ndarray # float64 + + +def compile_ms1_exporter(pwiz_dir: Optional[str] = None) -> Path: + """Compile the exporter once per source version; returns the .exe (cached).""" + dlls = find_pwiz_dir(pwiz_dir) + compiler = shutil.which("mcs") or shutil.which("csc") + if not shutil.which("mono") or not compiler: + raise RuntimeError("mono and mcs (or csc) are required to read Thermo RAW files") + out_dir = _cache_root() / "thermo" / f"ms1export_{EXPORTER_VERSION}" + exe = out_dir / "Ms1Export.exe" + if exe.is_file(): + return exe + out_dir.mkdir(parents=True, exist_ok=True) + src = out_dir / "Ms1Export.cs" + src.write_text(_MS1_EXPORT_CS) + refs = [f"-r:{dlls / name}" for name in DLLS] + tmp = out_dir / "Ms1Export.exe.tmp" + proc = subprocess.run([compiler, "-nologo", f"-out:{tmp}", *refs, str(src)], + capture_output=True, text=True) + if proc.returncode != 0: + raise RuntimeError(f"compiling the MS1 exporter failed:\n{proc.stdout}\n{proc.stderr}") + os.replace(tmp, exe) + return exe + + +def read_ms1x(path: str) -> Dict[int, RawMs1Scan]: + """Read the exporter's binary file into {scan: RawMs1Scan}.""" + out: Dict[int, RawMs1Scan] = {} + with open(path, "rb") as fh: + if fh.read(4) != b"MS1X": + raise ValueError(f"{path}: not an MS1 export") + (version,) = struct.unpack(" None: + """Write scans in the exporter's format (tests and hand-built fixtures).""" + with open(path, "wb") as fh: + fh.write(b"MS1X" + struct.pack(" str: + h = hashlib.sha256() + with open(path, "rb") as fh: + for block in iter(lambda: fh.read(chunk), b""): + h.update(block) + return h.hexdigest() + + +def export_ms1(raw_path: str, pwiz_dir: Optional[str] = None, out_path: Optional[str] = None, + timeout: Optional[float] = None) -> Dict[int, RawMs1Scan]: + """Every MS1 scan's centroids from a Thermo RAW file, as {scan: RawMs1Scan}. + + The RAW path is resolved first: RawFileReader can open a symlink yet fail to select the MS + device. ``out_path`` keeps the binary export; by default a temp file is used and removed. + """ + raw = os.path.realpath(raw_path) + if not os.path.isfile(raw): + raise FileNotFoundError(raw) + exe = compile_ms1_exporter(pwiz_dir) + dlls = find_pwiz_dir(pwiz_dir) + keep = out_path is not None + if out_path is None: + fd, out_path = tempfile.mkstemp(suffix=".ms1x") + os.close(fd) + try: + proc = subprocess.run(["mono", str(exe), str(dlls), raw, out_path], + capture_output=True, text=True, timeout=timeout) + if proc.returncode != 0: + raise RuntimeError(f"MS1 export failed for {raw} (exit {proc.returncode}):\n{proc.stderr}") + return read_ms1x(out_path) + finally: + if not keep and os.path.exists(out_path): + os.remove(out_path) diff --git a/tests/test_ms1_patch.py b/tests/test_ms1_patch.py new file mode 100644 index 0000000..8579594 --- /dev/null +++ b/tests/test_ms1_patch.py @@ -0,0 +1,295 @@ +"""patch_ms1 / ms1_cap_status / the MS1 export format -- synthetic, no RAW needed. + +A last test runs against a real Thermo RAW when MZML_UTILS_TEST_RAW names one whose spectrum +cache exists beside its mzML (it never writes: dry run only), and skips otherwise. +""" +import hashlib +import os +import shutil +import sqlite3 +import tempfile + +import numpy as np +import pytest + +from mzml_utils import spectrum_cache as sc +from mzml_utils import thermo_raw as tr +from mzml_utils.reader import Spectrum + +CAP = 60 + + +def _raw_scan(scan, rt, n, seed): + rng = np.random.default_rng(seed) + mz = np.sort(rng.uniform(400, 1800, n)) + it = rng.uniform(1e3, 1e7, n).astype(np.float32).astype(np.float64) # float32-exact + return tr.RawMs1Scan(scan, rt, 1, mz, it) + + +def _capped(x, cap=CAP): + keep = np.sort(np.argsort(x.intensity)[::-1][:cap]) + return x.mz[keep], x.intensity[keep] + + +def _spec(scan, level, mz, it, rt): + return Spectrum(scan_num=scan, mz=np.asarray(mz, float), intensity=np.asarray(it, float), + ms_level=level, rt=rt, filter_string="FTMS + p NSI Full ms" if level == 1 else "FTMS + c NSI d Full ms2 800.0@hcd30.00", + activation_type="MS1" if level == 1 else "HCD", precursor_mz=0.0 if level == 1 else 800.0, + precursor_charge=0 if level == 1 else 3) + + +def _make_cache(tmp, raw, n_ms2=3, tamper=None): + """A cache whose MS1 rows hold the RAW scans capped at CAP; MS2 rows in between.""" + db = os.path.join(tmp, "spectra_cache", "run.spectra.db") + os.makedirs(os.path.dirname(db), exist_ok=True) + conn = sc.open_db(db) + scan = 1 + for x in raw.values(): + mz, it = _capped(x) + if tamper: + mz, it = tamper(x.scan, mz, it) + sc.insert_scan(conn, _spec(x.scan, 1, mz, it, x.rt)) + for k in range(1, n_ms2 + 1): + sc.insert_scan(conn, _spec(x.scan + k, 2, [200.0 + k, 300.0], [10.0, 20.0], x.rt + 0.001 * k)) + sc.set_meta(conn, schema_version=1) + conn.commit() + conn.close() + return db + + +def _md5(p): + return hashlib.md5(open(p, "rb").read()).hexdigest() + + +@pytest.fixture +def raw(): + return {s: _raw_scan(s, s / 100.0, n, s) for s, n in ((1, 150), (5, 40), (9, 300), (13, CAP))} + + +def test_ms1x_round_trip(tmp_path, raw): + p = tmp_path / "x.ms1x" + tr.write_ms1x(str(p), raw.values()) + back = tr.read_ms1x(str(p)) + assert sorted(back) == sorted(raw) + for s in raw: + assert np.array_equal(back[s].mz, raw[s].mz) and np.array_equal(back[s].intensity, raw[s].intensity) + assert back[s].rt == raw[s].rt and back[s].source == 1 + + +def test_ms1x_truncated_raises(tmp_path, raw): + p = tmp_path / "x.ms1x" + tr.write_ms1x(str(p), raw.values()) + data = open(p, "rb").read() + open(p, "wb").write(data[:-20]) + with pytest.raises(ValueError): + tr.read_ms1x(str(p)) + + +def test_check_accepts_a_most_intense_cap(raw): + x = raw[1] + mz, it = _capped(x) + assert sc.check_ms1_against_raw(mz.astype(np.float32), it.astype(np.float32), x.mz, x.intensity) is None + # an uncapped scan (RAW <= cap) is identical + y = raw[5] + assert sc.check_ms1_against_raw(y.mz, y.intensity, y.mz, y.intensity) is None + + +def test_check_rejects_other_differences(raw): + x = raw[1] + mz, it = _capped(x) + shifted = mz.copy(); shifted[3] += 0.01 + assert "no RAW peak" in sc.check_ms1_against_raw(shifted, it, x.mz, x.intensity) + scaled = it.copy(); scaled[0] *= 1.01 + assert "intensity differs" in sc.check_ms1_against_raw(mz, scaled, x.mz, x.intensity) + # keep the weakest peaks instead of the strongest: not a most-intense cap + weak = np.sort(np.argsort(x.intensity)[:CAP]) + assert "not a most-intense cap" in sc.check_ms1_against_raw(x.mz[weak], x.intensity[weak], x.mz, x.intensity) + assert "empty" in sc.check_ms1_against_raw([], [], x.mz, x.intensity) + + +def test_patch_replaces_ms1_and_leaves_ms2(tmp_path, raw): + db = _make_cache(str(tmp_path), raw) + with sc.SpectrumCache(db, mirror_local=False) as c: + before_ms2 = {s: c.get_spectrum(s) for s in (2, 3, 4)} + assert sc.ms1_cap_status(c)["capped"] is True and sc.ms1_cap_status(c)["cap"] == CAP + rep = sc.patch_ms1(db, export=raw) + assert rep["status"] == "patched" + assert rep["ms1_truncated_before"] == 2 and rep["cap_undone"] == CAP + assert rep["ms1_peaks_after"] == sum(len(x.mz) for x in raw.values()) + with sc.SpectrumCache(db, mirror_local=False) as c: + for s, x in raw.items(): + sp = c.get_spectrum(s) + assert sp.n_peaks == len(x.mz) + assert np.allclose(sp.mz, x.mz, rtol=1e-7) and np.allclose(sp.intensity, x.intensity, rtol=1e-6) + for s, old in before_ms2.items(): + new = c.get_spectrum(s) + assert np.array_equal(new.mz, old.mz) and np.array_equal(new.intensity, old.intensity) + st = sc.ms1_cap_status(c) + assert st["patched"] is True and st["capped"] is False + assert c.meta()["ms1_patch"]["ms2_digest"] == rep["ms2_digest"] + + +def test_patch_is_refused_twice_without_force(tmp_path, raw): + db = _make_cache(str(tmp_path), raw) + sc.patch_ms1(db, export=raw) + with pytest.raises(sc.PatchError): + sc.patch_ms1(db, export=raw) + assert sc.patch_ms1(db, export=raw, force=True)["status"] == "patched" + + +@pytest.mark.parametrize("case", ["rt", "scan_set", "foreign_peak", "dry_run"]) +def test_failed_or_dry_patch_writes_nothing(tmp_path, raw, case): + tamper = None + exp = dict(raw) + if case == "rt": + exp[9] = tr.RawMs1Scan(9, raw[9].rt + 0.01, 1, raw[9].mz, raw[9].intensity) + elif case == "scan_set": + exp.pop(13) + elif case == "foreign_peak": + def tamper(scan, mz, it): + if scan == 9: + mz = mz.copy(); mz[0] = 401.23456 + order = np.argsort(mz); mz, it = mz[order], it[order] + return mz, it + db = _make_cache(str(tmp_path), raw, tamper=tamper) + h = _md5(db) + if case == "dry_run": + rep = sc.patch_ms1(db, export=exp, dry_run=True) + assert rep["status"].startswith("checked") and rep["cap_undone"] == CAP + else: + with pytest.raises(sc.PatchError): + sc.patch_ms1(db, export=exp) + assert _md5(db) == h + assert not [f for f in os.listdir(os.path.dirname(db)) if f.endswith(".tmp")] + + +def test_profile_scan_is_refused(tmp_path, raw): + db = _make_cache(str(tmp_path), raw) + exp = dict(raw) + exp[5] = tr.RawMs1Scan(5, raw[5].rt, 3, raw[5].mz, raw[5].intensity) + with pytest.raises(sc.PatchError, match="no centroids"): + sc.patch_ms1(db, export=exp) + + +def test_default_raw_is_found_beside_the_mzml(tmp_path): + (tmp_path / "run_1.raw").write_bytes(b"x") + assert sc._default_raw_for(str(tmp_path / "run_1_calibrated.mzML")) == str(tmp_path / "run_1.raw") + db = tmp_path / "spectra_cache" / "run_1.spectra.db" + assert sc._default_raw_for(str(db)) == str(tmp_path / "run_1.raw") + + +def test_mzml_header_cap_is_read(tmp_path): + p = tmp_path / "x.mzML" + p.write_text('' + '' + '') + + class FakeReader: + path = str(p) + st = sc.ms1_cap_status(FakeReader()) + assert st["capped"] is True and st["cap"] == 1000 and st["source"] == "mzml-header" + + +RAW = os.environ.get("MZML_UTILS_TEST_RAW") + + +@pytest.mark.skipif(not RAW or not tr.raw_reader_available(), reason="set MZML_UTILS_TEST_RAW to a Thermo RAW with a cache beside its mzML") +def test_real_raw_export_matches_its_capped_cache(): + db = sc.cache_path_for(os.path.splitext(RAW)[0] + ".mzML") + assert os.path.isfile(db), db + rep = sc.patch_ms1(db, RAW, dry_run=True) + assert rep["n_ms1"] > 0 and rep["ms1_peaks_after"] >= rep["ms1_peaks_before"] + + +def _shift_ppm(mz): + return -4.0 + 4.5 * np.tanh((mz - 1100.0) / 150.0) # smooth, like MSFragger's curve + + +def _calibrated_tamper(rng_seed=None): + def tamper(scan, mz, it): + if rng_seed is None: + return mz * (1 + _shift_ppm(mz) / 1e6), it + rng = np.random.default_rng(rng_seed + scan) + return mz * (1 + rng.uniform(-8, 8, len(mz)) / 1e6), it # not a function of m/z + return tamper + + +def _set_source(db, name): + conn = sqlite3.connect(db) + sc.set_meta(conn, source_mzml=f"/data/{name}") + conn.commit(); conn.close() + + +def test_calibrated_cache_is_uncapped_on_its_own_scale(tmp_path, raw): + db = _make_cache(str(tmp_path), raw, tamper=_calibrated_tamper()) + _set_source(db, "run_calibrated.mzML") + with sc.SpectrumCache(db, mirror_local=False) as c: + kept = {s: c.get_spectrum(s) for s in raw} + with pytest.raises(sc.PatchError, match="calibrated"): + sc.patch_ms1(db, export=raw, mz_mode="exact") + rep = sc.patch_ms1(db, export=raw) # auto -> calibrated by the source name + assert rep["mz_mode"] == "calibrated" and rep["heldout_p99_ppm_worst_scan"] < 0.5 + with sc.SpectrumCache(db, mirror_local=False) as c: + for s, x in raw.items(): + sp = c.get_spectrum(s) + assert sp.n_peaks == len(x.mz) + truth = x.mz * (1 + _shift_ppm(x.mz) / 1e6) + err = np.abs(np.sort(sp.mz.astype(float)) - np.sort(truth)) / truth * 1e6 + assert err.max() < 0.5 # float32 storage + interpolation + assert set(kept[s].mz.tolist()) <= set(sp.mz.tolist()) # kept peaks unchanged + + +def test_calibrated_mode_refuses_a_shift_that_is_not_smooth(tmp_path, raw): + db = _make_cache(str(tmp_path), raw, tamper=_calibrated_tamper(rng_seed=7)) + h = _md5(db) + with pytest.raises(sc.PatchError, match="not smooth|reorders"): + sc.patch_ms1(db, export=raw, mz_mode="calibrated") + assert _md5(db) == h + + +def test_cap_peaks_excludes_ties_at_the_threshold(): + mz = np.arange(6, dtype=float) + it = np.array([5.0, 9.0, 3.0, 3.0, 7.0, 1.0]) + assert cap_list(mz, it, 3) == [1.0, 4.0, 0.0] # 9, 7, 5 + assert cap_list(mz, it, 4) == [1.0, 4.0, 0.0] # the two 3s tie at the 4th place: both out + assert cap_list(mz, it, 5) == [1.0, 4.0, 0.0, 2.0, 3.0] + assert len(cap_list(mz, it, 10)) == 6 # fewer peaks than the cap: all kept + + +def cap_list(mz, it, n): + m, i = sc.cap_peaks(mz, it, n) + return [float(x) for x in m[np.argsort(-i, kind="stable")]] + + +def test_recapped_view_reproduces_the_capped_cache(tmp_path, raw): + db = _make_cache(str(tmp_path), raw) + with sc.SpectrumCache(db, mirror_local=False) as c: + before = {s: c.get_spectrum(s) for s in raw} + rep = sc.patch_ms1(db, export=raw) + assert rep["recap_reproduces_cache"] is True + with sc.SpectrumCache(db, mirror_local=False, ms1_cap=CAP) as c: + for s, old in before.items(): + new = c.get_spectrum(s) + assert np.array_equal(new.mz, old.mz) and np.array_equal(new.intensity, old.intensity) + assert sc.ms1_cap_status(c)["capped"] is True + with sc.open_spectra(db, ms1_cap=CAP) as c: + assert all(sp.n_peaks <= CAP for sp in c.iter_spectra(ms_level=1)) + + +def test_find_pwiz_dir_copies_the_dlls_locally(tmp_path, monkeypatch): + src = tmp_path / "bundle" + src.mkdir() + for name in tr.DLLS: + (src / name).write_bytes(b"dll") + monkeypatch.setenv("MZML_UTILS_CACHE", str(tmp_path / "cache")) + monkeypatch.setenv("THERMO_PWIZ_DIR", os.pathsep.join([str(tmp_path / "missing"), str(src)])) + monkeypatch.setattr(tr, "LOCAL_DLL_DIRS", ()) + found = tr.find_pwiz_dir() + assert found == tmp_path / "cache" / "thermo" / "dlls" + assert all((found / n).read_bytes() == b"dll" for n in tr.DLLS) + monkeypatch.delenv("THERMO_PWIZ_DIR") + assert tr.find_pwiz_dir() == found # the local copy now suffices + shutil.rmtree(found) + assert tr.find_pwiz_dir(required=False) is None + with pytest.raises(FileNotFoundError): + tr.find_pwiz_dir() diff --git a/tests/test_ms1quant.py b/tests/test_ms1quant.py new file mode 100644 index 0000000..43d25f1 --- /dev/null +++ b/tests/test_ms1quant.py @@ -0,0 +1,109 @@ +"""ms1quant: the IonQuant-like top-3 isotope-cluster estimator, on synthetic runs (no files).""" +import os + +import numpy as np +import pytest + +from mzml_utils import ms1quant as mq +from mzml_utils import spectrum_cache as sc +from mzml_utils.reader import Spectrum + +D = mq.ISOTOPE_SPACING +ANCHOR, Z = 900.0, 3 + + +def _cluster(scale, isotopes=(0, 1, 2, 3)): + rel = {-1: 0.0, 0: 1.0, 1: 0.9, 2: 0.5, 3: 0.2, 4: 0.07} + return [(ANCHOR + i * D / Z, scale * rel[i]) for i in isotopes] + + +def _run(tmp_path, profile, filler=None, ms2_rt=10.0): + """MS1 every 0.1 min from 9.0 to 11.0 with the cluster scaled by profile(rt); one MS2 at ms2_rt.""" + db = os.path.join(str(tmp_path), "spectra_cache", "run.spectra.db") + os.makedirs(os.path.dirname(db), exist_ok=True) + conn = sc.open_db(db) + scan = 1 + for k in range(21): + rt = round(9.0 + 0.1 * k, 3) + peaks = [(400.0 + j, 1e4 * (j + 1)) for j in range(filler or 5)] + s = profile(rt) + if s: + peaks += _cluster(s) + peaks.sort() + sc.insert_scan(conn, Spectrum(scan_num=scan, mz=np.array([p[0] for p in peaks]), + intensity=np.array([p[1] for p in peaks]), ms_level=1, rt=rt, + activation_type="MS1", filter_string="FTMS + p NSI Full ms")) + scan += 1 + if abs(rt - ms2_rt) < 1e-9: + sc.insert_scan(conn, Spectrum(scan_num=scan, mz=np.array([200.0]), intensity=np.array([1.0]), + ms_level=2, rt=rt + 0.001, activation_type="HCD", + precursor_mz=ANCHOR + D / Z, precursor_charge=Z, + filter_string="FTMS + c NSI d Full ms2 900.33@hcd30.00")) + ms2 = scan + scan += 1 + sc.set_meta(conn, schema_version=1) + conn.commit() + conn.close() + return db, ms2 + + +def gauss(rt): + return 1e6 * np.exp(-((rt - 10.0) / 0.12) ** 2) if abs(rt - 10.0) <= 0.3 else 0.0 + + +def test_top3_matches_a_hand_computation(tmp_path): + db, ms2 = _run(tmp_path, gauss) + with sc.open_spectra(db) as r: + q = mq.quantify_precursors(r, [dict(ms2_scan=ms2, charge=Z, anchor_mz=ANCHOR)])[0] + per_scan = sorted(gauss(round(9.0 + 0.1 * k, 3)) * (1.0 + 0.9 + 0.5 + 0.2) for k in range(21)) + assert q["top"] == pytest.approx(sum(per_scan[-3:]), rel=1e-6) + assert q["apex"] == pytest.approx(per_scan[-1], rel=1e-6) + assert q["n_scans"] == sum(v > 0 for v in per_scan) == 5 and q["apex_rt"] == pytest.approx(10.0) + assert q["cap_censored"] is False and q["apex_floor_ratio"] == float("inf") + + +def test_anchor_defaults_to_the_selected_ion(tmp_path): + db, ms2 = _run(tmp_path, gauss) + with sc.open_spectra(db) as r: + q = mq.quantify_precursors(r, [dict(ms2_scan=ms2, charge=Z)])[0] + assert q["anchor_source"] == "selected_ion" and q["anchor_mz"] == pytest.approx(ANCHOR + D / Z) + assert q["top"] > 0 # iso_lo = -1 still reaches the monoisotope + + +def test_min_scans_and_min_isotopes(tmp_path): + narrow = lambda rt: 1e6 if abs(rt - 10.0) < 0.15 else 0.0 # 3 scans: 9.9, 10.0, 10.1 + db, ms2 = _run(tmp_path, narrow) + with sc.open_spectra(db) as r: + idx = mq.Ms1Index(r) + rtc = idx.scan_rt[ms2] + assert mq.precursor_ms1(idx, ANCHOR, Z, rtc).n_scans == 3 + assert mq.precursor_ms1(idx, ANCHOR, Z, rtc, mq.Ms1QuantParams(min_scans=4)).top == 0.0 + # one isotope only is not a cluster + assert mq.precursor_ms1(idx, ANCHOR, Z, rtc, mq.Ms1QuantParams(min_isotopes=5)).top == 0.0 + # a wrong charge finds nothing (z = 5: no spacing D/5 multiple below +4 lands on a z = 3 isotope) + assert mq.precursor_ms1(idx, ANCHOR, 5, rtc).top == 0.0 + + +def test_nearest_peak_rule(): + mz = np.array([899.9950, 900.0000, 900.0060]) + assert mq._nearest_peak(mz, 900.0010, 900 * 10e-6) == 1 + assert mq._nearest_peak(mz, 900.0030, 900 * 10e-6) == 2 # tie -> the higher m/z + assert mq._nearest_peak(mz, 900.0200, 900 * 10e-6) == -1 + + +def test_censoring_on_a_capped_view(tmp_path): + weak = lambda rt: 3e4 if abs(rt - 10.0) <= 0.3 else 0.0 # near the filler peaks + db, ms2 = _run(tmp_path, weak, filler=40) + with sc.open_spectra(db, ms1_cap=30) as r: + assert sc.ms1_cap_status(r)["capped"] + q = mq.quantify_precursors(r, [dict(ms2_scan=ms2, charge=Z, anchor_mz=ANCHOR)])[0] + assert q["cap_censored"] is True + + +def test_params_from_workflow(tmp_path): + p = tmp_path / "fragpipe.workflow" + p.write_text("# comment\nionquant.mztol=12\nionquant.rttol=0.5\nionquant.minscans=2\n" + "ionquant.minisotopes=1\nmsfragger.precursor_true_tolerance=20\n") + prm = mq.Ms1QuantParams.from_workflow(str(p)) + assert (prm.mztol_ppm, prm.rttol_min, prm.min_scans, prm.min_isotopes) == (12.0, 0.5, 2, 1) + assert prm.topn == 3