diff --git a/CHANGELOG.md b/CHANGELOG.md index ce025b9..0231dd0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,17 @@ All notable changes to this project are documented here. The format follows ## [Unreleased] ### Added +- `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 + result for a DataFrame. The Python counterpart of `chromExtract()` in the Bioconductor + package Chromatograms (Louail et al., Anal Chem 2026). +- `peak_boundary()`, `peak_metrics()` / `PeakMetrics`, `beta_values()`: peak boundaries, + in-boundary area, FWHM, asymmetry, S/N, prominence and Kumler's beta shape metrics, ported + from the Chromatograms, MsQuality and MetaboCoreUtils R sources and pinned against them + (`tests/data/r_peak_reference.json`). Deviations from R are opt-in: `contains_rt=` (the + peak containing an MS2 time), `apex_window=`, `nearest_valley=` (R's left boundary runs to + the start of a zero-filled trace). - `iter_spectra(ms_level=, rt_range=)` on both readers; `SpectrumCache` filters in SQL, so an MS1 sweep never decodes the MS2 scans. `extract_xics` pushes its filters down when the reader accepts them. @@ -33,6 +44,8 @@ All notable changes to this project are documented here. The format follows - Main-only git workflow: pre-commit branch guard and `scripts/dev-setup.sh`. ### Fixed +- `XIC.area` docstring said "peak area"; it integrates the whole trace (use `peak_metrics` + for the area between peak boundaries). - XIC area under numpy ≥ 2 (the `trapz` fallback was evaluated eagerly). - `pyteomics[xml]` is declared as the dependency; a plain `pyteomics` install could not import the package. - README development-install clone URL pointed at a non-existent GitHub user (#5). diff --git a/README.md b/README.md index fcf5cf3..e079000 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,8 @@ pip install -e ".[dev]" |--------|-------------| | `reader` | mzML file I/O with indexed and sequential access | | `ions` | Ion searching, mass matching, tolerance calculations | -| `xic` | Extracted ion chromatograms (intensity vs retention time) over a run; isotope co-elution evidence for a precursor | +| `xic` | Extracted ion chromatograms (intensity vs retention time) over a run; isotope co-elution evidence for a precursor; peak boundaries and peak-quality metrics | +| `chrom_extract` | Peak table → chromatograms across many runs, every row column carried through (the Bioconductor `chromExtract` analogue) | | `isotopes` | Exact-composition isotope distributions, for envelopes the averagine cannot express (halogen tags, heavy labels) | | `pairing` | MS1 duty-cycle grouping, HCD-EThcD scan pairing | | `fragments` | Fragment ion calculator (b/y/c/z/Y/oxonium), peak matching, false match rate | @@ -99,6 +100,27 @@ dist = isotope_distribution(comp, n_peaks=3, charge=2) print(dist.mz, dist.abundance) # pass distribution=dist to isotope_trace_evidence ``` +### Extract a peak table across many runs, with peak metrics + +```python +import pandas as pd +from mzml_utils import extract_chromatograms, chromatogram_rows + +# One row per target: m/z, RT window, run name, plus any annotation columns +table = psms[["psg", "run", "mz", "rt_min", "rt_max", "ms2_rt"]].to_dict("records") +res = extract_chromatograms(table, mzml_paths, by="run", apex_rt="ms2_rt", + tolerance=10.0, nearest_valley=True) +df = pd.DataFrame(chromatogram_rows(res)) # row columns + apex, boundaries, area, FWHM, ... +``` + +Each row gets its own m/z and RT window, and one pass per run serves every row of that +run. Without `by=`, every row is extracted from every run (cross-sample overlays). +`peak_metrics` / `peak_boundary` port the Chromatograms, MsQuality and MetaboCoreUtils R +code (boundaries, S/N, prominence, Kumler beta shape) and are pinned against those R +functions. `apex_rt=` measures the peak whose boundaries contain the MS2 time, and +`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. + ### Calculate fragment ions ```python diff --git a/src/mzml_utils/__init__.py b/src/mzml_utils/__init__.py index 549521c..5fc92a8 100644 --- a/src/mzml_utils/__init__.py +++ b/src/mzml_utils/__init__.py @@ -99,6 +99,18 @@ coelution_score, isotope_trace_evidence, IsotopeTraceEvidence, + peak_boundary, + peak_metrics, + PeakMetrics, + beta_values, +) + +# chrom_extract -- peak table -> chromatograms across many runs (chromExtract analogue) +from .chrom_extract import ( + extract_chromatograms, + chromatogram_rows, + ExtractedChromatogram, + source_key, ) # canonical match layer @@ -207,6 +219,15 @@ "coelution_score", "isotope_trace_evidence", "IsotopeTraceEvidence", + "peak_boundary", + "peak_metrics", + "PeakMetrics", + "beta_values", + # chrom_extract (peak table -> chromatograms across runs) + "extract_chromatograms", + "chromatogram_rows", + "ExtractedChromatogram", + "source_key", # constants "PROTON", "H2O", diff --git a/src/mzml_utils/chrom_extract.py b/src/mzml_utils/chrom_extract.py new file mode 100644 index 0000000..4a273c8 --- /dev/null +++ b/src/mzml_utils/chrom_extract.py @@ -0,0 +1,299 @@ +""" +Table-driven chromatogram extraction across many runs. + +The Python counterpart of ``chromExtract()`` in the Bioconductor package +Chromatograms (Louail, Gatto, Gibb & Rainer, Anal Chem 2026). A peak table from any +software -- MSFragger PSMs, a published supplementary table, a hand-written target +list -- names an m/z window and a retention-time window per row. Each row's +chromatogram is pulled from the raw data, and every column of the row is carried +through unchanged, so the result stays traceable to the table that asked for it. + +What differs from :func:`~mzml_utils.xic.extract_xics`: + +* every row has its OWN retention-time window (and m/z window), so glycoforms that + elute minutes apart are extracted in one pass without a whole-run trace each; +* one pass per run serves all rows of that run: each MS1 scan is decoded once, and + all active rows are summed in one vectorised ``searchsorted`` over a cumulative + intensity sum instead of a mask per target; +* many runs in one call, with ``by=`` restricting a row to its own run (a PSM + belongs to the file it was identified in) or, without ``by``, every row extracted + from every run (cross-sample overlays: one glycoform across all fractions); +* optional :func:`~mzml_utils.xic.peak_metrics` per chromatogram. + +Typical use:: + + from mzml_utils import extract_chromatograms, chromatogram_rows + + table = [ + {"psg": "P01024_N1115_N4H5F1A1", "mz": 1204.564, "rt_min": 61.5, "rt_max": 63.5, + "run": "fraction_07", "ms2_rt": 62.4}, + ... + ] + res = extract_chromatograms(table, ["/data/fraction_07.mzML", ...], by="run", + apex_rt="ms2_rt", nearest_valley=True) + rows = chromatogram_rows(res) # list of flat dicts; pandas.DataFrame(rows) +""" + +from __future__ import annotations + +import inspect +import os +import warnings +from dataclasses import dataclass +from typing import Any, Dict, Iterable, List, Mapping, Optional, Sequence, Tuple, Union + +import numpy as np + +from .xic import XIC, PeakMetrics, peak_metrics + +SourceSpec = Union[Sequence[Union[str, os.PathLike]], Mapping[str, Any]] + +_SOURCE_SUFFIXES = (".spectra.db", ".mzml", ".mzxml", ".db") + + +def source_key(path: Union[str, os.PathLike]) -> str: + """The run name of a spectrum source: its file name without the mzML/cache suffix. + + ``/data/spectra_cache/run_01.spectra.db`` and ``/data/run_01.mzML`` both give + ``run_01``; this is what a ``by=`` column must contain. + """ + name = os.path.basename(str(path)) + low = name.lower() + for suf in _SOURCE_SUFFIXES: + if low.endswith(suf): + return name[: len(name) - len(suf)] + return name + + +@dataclass +class ExtractedChromatogram: + """One row of the peak table, extracted from one run.""" + + row: Dict[str, Any] + """The peak-table row, carried through unchanged.""" + source: str + """Run name (:func:`source_key`, or the key given in a ``sources`` mapping).""" + xic: XIC + metrics: Optional[PeakMetrics] = None + + +def _normalize_sources(sources: SourceSpec) -> List[Tuple[str, Any]]: + if isinstance(sources, Mapping): + return [(str(k), v) for k, v in sources.items()] + if isinstance(sources, (str, os.PathLike)): + sources = [sources] + out = [(source_key(p), p) for p in sources] + keys = [k for k, _ in out] + if len(set(keys)) != len(keys): + dup = sorted({k for k in keys if keys.count(k) > 1}) + raise ValueError(f"sources share a run name {dup}; pass a {{name: path}} mapping") + return out + + +def _row_windows(rows: List[Mapping[str, Any]], mz: str, mz_min: str, mz_max: str, + rt_min: str, rt_max: str, tolerance: float, unit: str): + """Arrays (mz_lo, mz_hi, rt_lo, rt_hi, target_mz), validated like chromExtract.""" + if unit not in ("ppm", "Da"): + raise ValueError(f"Unknown unit '{unit}'. Use 'ppm' or 'Da'.") + n = len(rows) + lo = np.empty(n) + hi = np.empty(n) + target = np.empty(n) + rlo = np.empty(n) + rhi = np.empty(n) + for i, r in enumerate(rows): + try: + rlo[i] = float(r[rt_min]) + rhi[i] = float(r[rt_max]) + except KeyError as e: + raise ValueError(f"peak table row {i} has no {e.args[0]!r} column") from None + if np.isnan(rlo[i]) or np.isnan(rhi[i]): + raise ValueError(f"peak table row {i}: {rt_min}/{rt_max} is missing (NaN)") + if rlo[i] > rhi[i]: + raise ValueError(f"peak table row {i}: {rt_min} > {rt_max}") + has_min, has_max = mz_min in r, mz_max in r + if has_min != has_max: + raise ValueError(f"peak table row {i}: give both {mz_min} and {mz_max}, or neither") + if has_min: + lo[i], hi[i] = float(r[mz_min]), float(r[mz_max]) + target[i] = float(r[mz]) if mz in r else (lo[i] + hi[i]) / 2.0 + else: + if mz not in r: + raise ValueError(f"peak table row {i} has neither {mz!r} nor {mz_min!r}/{mz_max!r}") + t = float(r[mz]) + tol = t * tolerance / 1e6 if unit == "ppm" else tolerance + lo[i], hi[i], target[i] = t - tol, t + tol, t + if np.isnan(lo[i]) or np.isnan(hi[i]): + raise ValueError(f"peak table row {i}: m/z is missing (NaN)") + return lo, hi, rlo, rhi, target + + +def _iter_pushed(reader, ms_level: Optional[int], rt_range: Tuple[float, float]): + """iter_spectra with the MS-level / RT filters pushed down when accepted.""" + accepted = inspect.signature(reader.iter_spectra).parameters + pushed = {k: v for k, v in (("ms_level", ms_level), ("rt_range", rt_range)) + if v is not None and k in accepted} + for spec in reader.iter_spectra(**pushed): + if ms_level is not None and spec.ms_level != ms_level: + continue + if not (rt_range[0] <= spec.rt <= rt_range[1]): + continue + yield spec + + +def _extract_one_source(reader, idx: np.ndarray, lo, hi, rlo, rhi, ms_level): + """(rt, intensity) per row index in ``idx`` from one open reader, in one pass.""" + s_lo, s_hi, s_rlo, s_rhi = lo[idx], hi[idx], rlo[idx], rhi[idx] + rt_parts: List[np.ndarray] = [] + row_parts: List[np.ndarray] = [] + val_parts: List[np.ndarray] = [] + for spec in _iter_pushed(reader, ms_level, (float(s_rlo.min()), float(s_rhi.max()))): + rt = float(spec.rt) + active = np.flatnonzero((s_rlo <= rt) & (rt <= s_rhi)) + if active.size == 0: + continue + mz_a = np.asarray(spec.mz, dtype=float) + int_a = np.asarray(spec.intensity, dtype=float) + if mz_a.size and np.any(np.diff(mz_a) < 0): + order = np.argsort(mz_a, kind="stable") + mz_a, int_a = mz_a[order], int_a[order] + csum = np.concatenate(([0.0], np.cumsum(int_a))) + a = np.searchsorted(mz_a, s_lo[active], side="left") + b = np.searchsorted(mz_a, s_hi[active], side="right") + val_parts.append(csum[b] - csum[a]) + row_parts.append(active) + rt_parts.append(np.full(active.size, rt)) + + out: Dict[int, Tuple[np.ndarray, np.ndarray]] = {} + if row_parts: + rows = np.concatenate(row_parts) + rts = np.concatenate(rt_parts) + vals = np.concatenate(val_parts) + order = np.lexsort((rts, rows)) # by row, then retention time (stable) + rows, rts, vals = rows[order], rts[order], vals[order] + cuts = np.flatnonzero(np.diff(rows)) + 1 + for r_rows, r_rt, r_val in zip(np.split(rows, cuts), np.split(rts, cuts), + np.split(vals, cuts)): + out[int(idx[r_rows[0]])] = (r_rt, r_val) + return out + + +def extract_chromatograms(peak_table: Iterable[Mapping[str, Any]], + sources: SourceSpec, *, + by: Optional[str] = None, + mz: str = "mz", + mz_min: str = "mz_min", + mz_max: str = "mz_max", + rt_min: str = "rt_min", + rt_max: str = "rt_max", + tolerance: float = 10.0, + unit: str = "ppm", + ms_level: Optional[int] = 1, + metrics: bool = True, + apex_rt: Optional[str] = None, + **boundary_kwargs) -> List[ExtractedChromatogram]: + """Extract one chromatogram per (peak-table row, run), carrying the row through. + + Args: + peak_table: Rows as mappings (``df.to_dict("records")`` for a DataFrame). Each + row needs ``rt_min``/``rt_max`` and either ``mz`` (window = +/- ``tolerance``) + or both ``mz_min``/``mz_max``. A missing or NaN window is an error, as in + chromExtract -- a row is never silently extracted over the whole run. + sources: Spectrum files (mzML or ``.spectra.db``; opened with + :func:`~mzml_utils.open_spectra`, so a co-located cache is used), or a + ``{run_name: path_or_open_reader}`` mapping. Readers passed in are not closed. + by: Row column holding the run name. When given, a row is extracted only from + the run whose name equals it; a run name matching no source is an error. + When ``None``, every row is extracted from every source. + mz, mz_min, mz_max, rt_min, rt_max: Column names. + tolerance, unit: m/z window for rows given by ``mz`` (``'ppm'`` or ``'Da'``). + ms_level: MS level extracted (1 = precursor traces). ``None`` uses every scan. + metrics: Compute :func:`~mzml_utils.xic.peak_metrics` for each chromatogram. + apex_rt: Optional column of retention times the measured peak must contain + (e.g. the PSM's MS2 scan time); passed as ``contains_rt`` to + :func:`~mzml_utils.xic.peak_boundary`. Without it the tallest peak in the + window is measured. A row whose value is missing falls back to the tallest. + **boundary_kwargs: ``threshold``, ``baseline_threshold``, ``baseline_quantile``, + ``nearest_valley`` passed to :func:`~mzml_utils.xic.peak_boundary`. + ``nearest_valley=True`` is recommended for these zero-filled traces. + + Returns: + One :class:`ExtractedChromatogram` per (row, run) pair, in row order then source + order. A row whose window holds no scan gets an empty XIC and empty metrics. + Intensities within an m/z window are summed (the ``mode='sum'`` of + ``extract_xics``); they agree with ``extract_xics`` to floating-point rounding. + """ + rows = [dict(r) for r in peak_table] + srcs = _normalize_sources(sources) + if not rows or not srcs: + return [] + lo, hi, rlo, rhi, target = _row_windows(rows, mz, mz_min, mz_max, rt_min, rt_max, + tolerance, unit) + + names = [k for k, _ in srcs] + if by is not None: + missing = [i for i, r in enumerate(rows) if by not in r] + if missing: + raise ValueError(f"peak table rows {missing[:5]} have no {by!r} column") + unknown = sorted({str(r[by]) for r in rows} - set(names)) + if unknown: + raise ValueError(f"{by!r} values match no source: {unknown[:5]} " + f"(sources: {names[:5]}{' ...' if len(names) > 5 else ''})") + + traces: Dict[Tuple[int, str], Tuple[np.ndarray, np.ndarray]] = {} + for name, src in srcs: + if by is None: + idx = np.arange(len(rows)) + else: + idx = np.array([i for i, r in enumerate(rows) if str(r[by]) == name], dtype=int) + if idx.size == 0: + continue + if hasattr(src, "iter_spectra"): + got = _extract_one_source(src, idx, lo, hi, rlo, rhi, ms_level) + else: + from .spectrum_cache import open_spectra + with open_spectra(str(src)) as reader: + got = _extract_one_source(reader, idx, lo, hi, rlo, rhi, ms_level) + for i, tr in got.items(): + traces[(i, name)] = tr + + empty = (np.array([], dtype=float), np.array([], dtype=float)) + out: List[ExtractedChromatogram] = [] + for i, row in enumerate(rows): + for name in names: + if by is not None and str(row[by]) != name: + continue + rt_arr, int_arr = traces.get((i, name), empty) + xic = XIC(target_mz=float(target[i]), rt=rt_arr, intensity=int_arr, + name=f"{name}:{i}", ms_level=ms_level if ms_level is not None else 0) + pm = None + if metrics: + t = row.get(apex_rt) if apex_rt is not None else None + t = None if t is None or np.isnan(float(t)) else float(t) + pm = peak_metrics(xic, contains_rt=t, **boundary_kwargs) + out.append(ExtractedChromatogram(row=row, source=name, xic=xic, metrics=pm)) + return out + + +def chromatogram_rows(results: Sequence[ExtractedChromatogram]) -> List[Dict[str, Any]]: + """Flatten results to one dict per chromatogram: the row's own columns, then + ``source``, ``n_scans`` and the metric columns (``pandas.DataFrame(rows)`` ready). + + A row column that shares a name with an added column is replaced, with a warning + (chromExtract does the same). + """ + flat: List[Dict[str, Any]] = [] + warned = set() + for res in results: + d = dict(res.row) + added: Dict[str, Any] = {"source": res.source, "n_scans": res.xic.n_points} + if res.metrics is not None: + added.update(res.metrics.as_dict()) + for k, v in added.items(): + if k in d and k not in warned: + warnings.warn(f"chromatogram_rows: peak-table column {k!r} replaced by the " + f"computed value", stacklevel=2) + warned.add(k) + d[k] = v + flat.append(d) + return flat diff --git a/src/mzml_utils/xic.py b/src/mzml_utils/xic.py index c441f1a..ef14d21 100644 --- a/src/mzml_utils/xic.py +++ b/src/mzml_utils/xic.py @@ -38,6 +38,7 @@ from __future__ import annotations import inspect +import math import numpy as np from dataclasses import dataclass from typing import Dict, List, Mapping, Optional, Sequence, Tuple, Union @@ -82,23 +83,31 @@ def apex_rt(self) -> float: @property def area(self) -> float: - """Trapezoidal peak area over retention time.""" + """Trapezoidal area over the WHOLE trace (every point in the extraction window). + + This is not a peak area: a second elution or background inside the window is + integrated too. For the area between the peak's boundaries use + :func:`peak_metrics`. + """ if len(self.rt) < 2: return 0.0 return float(_trapz(self.intensity, self.rt)) - def _half_max_span(self): + def _half_max_span(self, apex: Optional[int] = None): """(apex_idx, half, lo_idx, hi_idx) of the CONTIGUOUS half-maximum region. Contiguous on purpose: an XIC often contains more than one peak (a second elution, or an interfering species at the same m/z), and counting every point above half maximum anywhere in the trace would merge them into one absurdly wide peak. Only the region containing the apex is the peak. + + ``apex`` is the index of the peak to measure; ``None`` takes the tallest point. """ n = len(self.intensity) if n == 0: return None - apex = int(np.argmax(self.intensity)) + if apex is None: + apex = int(np.argmax(self.intensity)) half = float(self.intensity[apex]) / 2.0 if half <= 0.0: return None @@ -110,21 +119,15 @@ def _half_max_span(self): hi += 1 return apex, half, lo, hi - @property - def fwhm(self) -> float: - """Full width at half maximum, in the retention-time unit of `rt` (minutes). - - The crossings are linearly interpolated between the bracketing scans rather - than snapped to them, because at 9-15 points across a peak, snapping - quantises the width in steps of a whole cycle time. + def _half_max_crossings(self, apex: Optional[int] = None): + """(left_rt, right_rt, lo_idx, hi_idx) where the trace crosses half maximum. - Returns 0.0 when there is no determinable peak (empty trace, flat, or a - single point). Check `fwhm_is_truncated` before treating the value as a - measurement: a peak running off either end of the trace yields a LOWER BOUND. + The crossings are linearly interpolated between the bracketing scans. At an + end of the trace the crossing is snapped to that end point (a lower bound). """ - span = self._half_max_span() + span = self._half_max_span(apex) if span is None or len(self.rt) < 2: - return 0.0 + return None _apex, half, lo, hi = span def cross(i_in, i_out): @@ -137,6 +140,24 @@ def cross(i_in, i_out): left = cross(lo, lo - 1) if lo > 0 else float(self.rt[lo]) right = cross(hi, hi + 1) if hi < len(self.rt) - 1 else float(self.rt[hi]) + return left, right, lo, hi + + @property + def fwhm(self) -> float: + """Full width at half maximum, in the retention-time unit of `rt` (minutes). + + The crossings are linearly interpolated between the bracketing scans rather + than snapped to them, because at 9-15 points across a peak, snapping + quantises the width in steps of a whole cycle time. + + Returns 0.0 when there is no determinable peak (empty trace, flat, or a + single point). Check `fwhm_is_truncated` before treating the value as a + measurement: a peak running off either end of the trace yields a LOWER BOUND. + """ + cr = self._half_max_crossings() + if cr is None: + return 0.0 + left, right, _lo, _hi = cr return max(0.0, right - left) @property @@ -376,6 +397,329 @@ def extract_xic(source, target_mz: float, *, return cs.xics[nm] +# --------------------------------------------------------------------------- +# Peak boundaries and peak-quality metrics +# --------------------------------------------------------------------------- +# +# Ported from the R sources, not their documentation, so the numbers match: +# * boundaries -- Chromatograms (Bioconductor, Louail et al. Anal Chem 2026), +# R/helpers.R `.peak_boundary_one` +# * snr, prominence -- MsQuality R/function_Chromatograms_metrics.R +# (`signalToNoiseRatio`, `.peakProminenceSingle`) +# * beta_cor, beta_snr -- MetaboCoreUtils R/peak-shape-quality.R `betaValues` +# (Kumler, Hazelton & Ingalls, BMC Bioinformatics 2023, 24:404) +# tests/test_xic.py pins the port against values produced by running those R +# functions verbatim. + +_MAD_CONSTANT = 1.4826 # R stats::mad default: consistent with the SD of a normal + + +def _bounds_from_apex(inten: np.ndarray, apex: int, baseline: float, threshold: float, + baseline_threshold: float, nearest_valley: bool) -> Tuple[int, int]: + """Boundary indices of the peak at ``apex`` (steps 2-4 of :func:`peak_boundary`).""" + n = len(inten) + height = float(inten[apex]) - baseline + valley_ok = baseline + height * baseline_threshold + if nearest_valley: + # Leave a flat top, then descend while strictly falling: the first point whose + # outward neighbour is not lower is the valley's nearest point. + left = apex + while left > 0 and inten[left - 1] == inten[apex]: + left -= 1 + while left > 0 and inten[left - 1] < inten[left]: + left -= 1 + right = apex + while right < n - 1 and inten[right + 1] == inten[apex]: + right += 1 + while right < n - 1 and inten[right + 1] < inten[right]: + right += 1 + else: + left = apex + while left > 0 and not (inten[left] < inten[left - 1] + and (left == n - 1 or inten[left] <= inten[left + 1])): + left -= 1 + right = apex + while right < n - 1 and not ((right == 0 or inten[right] < inten[right - 1]) + and inten[right] <= inten[right + 1]): + right += 1 + + if not (inten[left] <= valley_ok and inten[right] <= valley_ok): + thresh = baseline + height * threshold + left_cand = np.flatnonzero(inten[:apex + 1] <= thresh) + right_cand = np.flatnonzero(inten[apex:] <= thresh) + left = int(left_cand[-1]) if left_cand.size else 0 + right = apex + int(right_cand[0]) if right_cand.size else n - 1 + return int(left), int(right) + + +def peak_boundary(xic: "XIC", *, + contains_rt: Optional[float] = None, + apex_window: Optional[Tuple[float, float]] = None, + threshold: float = 0.1, + baseline_threshold: float = 0.1, + baseline_quantile: float = 0.1, + nearest_valley: bool = False) -> Optional[Tuple[int, int, int]]: + """Left boundary, apex and right boundary of one chromatographic peak, as indices. + + The algorithm of ``Chromatograms::peakBoundary``: + + 1. baseline = the ``baseline_quantile`` quantile of the trace's intensities; + 2. walk outwards from the apex to the first local minimum (valley) on each side; + 3. keep the two valleys if BOTH are at or below + ``baseline + baseline_threshold * (apex - baseline)``; + 4. otherwise -- for either side failing, both sides are redone -- take the nearest + points on each side at or below ``baseline + threshold * (apex - baseline)``, + or the trace ends. + + With the defaults this reproduces R exactly: the apex is the tallest point in the + trace. Three options deviate from R, all off by default. + + ``contains_rt`` selects the peak that CONTAINS a retention time, e.g. the PSM's MS2 + scan time. A glycopeptide XIC often holds more than one peak (an isomer, a +/-1 Da + neighbour, a second elution), and the tallest need not be the assigned species. + Peaks are taken tallest first; each is bounded as in steps 2-4 on the whole trace, + and the first whose boundaries enclose ``contains_rt`` is returned, with its apex + re-taken as the tallest point inside those boundaries. A peak that does not enclose + it is set aside and the next tallest remaining point is tried. When ``contains_rt`` + lies on a tail below ``threshold`` of the main apex, R's boundaries exclude it; + the enclosing boundaries then come from a tail point's own, lower threshold and are + wider than R would draw for the apex they are reported with. Selecting by + containment rather than by a window around the MS2 time matters in DDA: the MS2 + can be triggered on a peak's tail, a minute or more from its apex, and a window + would then measure a point on the tail as if it were the apex. + + ``apex_window=(lo, hi)`` takes the apex as the tallest point inside the window. + + ``nearest_valley=True`` changes the walk on a FLAT valley. R's walk is asymmetric + there: the right side stops at the first flat point, but the left side walks + through the whole flat run (it keeps a flat minimum at its first, leftmost index). + A peak on a perfectly flat floor -- in particular a zero-filled XIC, where every + scan without the ion is exactly 0 -- therefore gets its left boundary at the START + OF THE TRACE, and every in-boundary metric is computed over the empty floor (a + perfect Gaussian on a flat floor scores beta_cor 0.35 in R). With + ``nearest_valley=True`` both sides stop at the flat run's point nearest the apex. + Traces with a noisy baseline rarely have exact ties and give the same result + either way. + + Any valley walk, R's included, stops at an internal DROPOUT: a scan inside the + peak where the ion is absent (intensity 0) is a valley at the baseline. Check the + trace before reading a boundary as the end of elution. + + Returns ``None`` when there are fewer than 3 points, the apex intensity is zero, + the apex window holds no point, or no peak above baseline contains + ``contains_rt``. The trace must not contain NaN; extraction writes 0 for an absent + ion. + """ + if contains_rt is not None and apex_window is not None: + raise ValueError("peak_boundary: give contains_rt or apex_window, not both") + inten = np.asarray(xic.intensity, dtype=float) + rt = np.asarray(xic.rt, dtype=float) + n = len(inten) + if np.isnan(inten).any(): + raise ValueError("peak_boundary: intensity contains NaN (write 0 for an absent ion)") + if n < 3: + return None + baseline = float(np.quantile(inten, baseline_quantile)) # numpy 'linear' == R type 7 + args = (baseline, threshold, baseline_threshold, nearest_valley) + + if contains_rt is not None: + t = float(contains_rt) + if not (rt[0] <= t <= rt[-1]): + return None + free = np.ones(n, dtype=bool) + while True: + cand = np.flatnonzero(free & (inten > baseline)) + if cand.size == 0: + return None + apex = int(cand[int(np.argmax(inten[cand]))]) + left, right = _bounds_from_apex(inten, apex, *args) + if rt[left] <= t <= rt[right]: + # The apex is the tallest point inside the accepted boundaries. It differs + # from `apex` when `t` lies on a tail below the main peak's threshold: the + # tail point's own (lower) threshold then draws boundaries that span the + # main peak as well. + top = left + int(np.argmax(inten[left:right + 1])) + return left, top, right + free[left:right + 1] = False + + if apex_window is None: + apex = int(np.argmax(inten)) + else: + inside = np.flatnonzero((rt >= apex_window[0]) & (rt <= apex_window[1])) + if inside.size == 0: + return None + apex = int(inside[int(np.argmax(inten[inside]))]) + if inten[apex] <= baseline: + return None + if inten[apex] == 0.0: + return None + left, right = _bounds_from_apex(inten, apex, *args) + return left, apex, right + + +def _scale_zero_one(v: np.ndarray) -> np.ndarray: + return (v - v.min()) / (v.max() - v.min()) + + +def _beta_pdf(x: np.ndarray, a: float, b: float) -> np.ndarray: + """Beta(a, b) density on [0, 1] (a, b > 1, so it is 0 at both ends).""" + log_norm = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + out = np.zeros_like(x, dtype=float) + inner = (x > 0.0) & (x < 1.0) + xi = x[inner] + out[inner] = np.exp(log_norm + (a - 1.0) * np.log(xi) + (b - 1.0) * np.log1p(-xi)) + return out + + +def beta_values(intensity, rtime=None, + skews: Sequence[float] = (3.0, 3.5, 4.0, 4.5, 5.0)) -> Tuple[float, float]: + """Kumler's beta peak-shape metrics: ``(beta_cor, beta_snr)``. + + A port of ``MetaboCoreUtils::betaValues``. Zero intensities are dropped, retention + time is rescaled to [0, 1], and the trace is correlated (Pearson) with Beta(skew, 5) + densities for each ``skew`` (values below 5 are right-skewed, i.e. tailing). + + * ``beta_cor`` -- the best correlation. Near 1 for a clean single peak; low for + noise or a trace with several maxima. + * ``beta_snr`` -- ``log10(max(intensity) / sd(diff(scaled best curve - scaled + intensity)))``. HIGHER is better. MsQuality labels this column + ``gaussian_residuals`` and documents it as a residual SD where lower is better; + its code returns this log10 ratio. Note the numerator is the RAW maximum + intensity while the noise is computed on 0-1 scaled curves, so the value grows + with absolute signal and is not a pure shape score. + + Returns ``(nan, nan)`` with fewer than 5 non-zero points (R returns NA), or when + the non-zero intensities are all equal (R's correlation is NA there). + """ + y = np.asarray(intensity, dtype=float) + x = np.arange(1, len(y) + 1, dtype=float) if rtime is None else np.asarray(rtime, dtype=float) + keep = y > 0 + y, x = y[keep], x[keep] + if len(y) < 5 or y.max() == y.min(): + return float("nan"), float("nan") + xs = _scale_zero_one(x) + curves = [_beta_pdf(xs, float(s), 5.0) for s in skews] + cors = [float(np.corrcoef(y, c)[0, 1]) for c in curves] + best = int(np.argmax(cors)) # first maximum, as R's which.max + noise = float(np.std(np.diff(_scale_zero_one(curves[best]) - _scale_zero_one(y)), ddof=1)) + with np.errstate(divide="ignore"): + beta_snr = float(np.log10(y.max() / noise)) if noise > 0 else float("inf") + return cors[best], beta_snr + + +@dataclass +class PeakMetrics: + """Boundaries and quality metrics of one chromatographic peak (see :func:`peak_metrics`). + + Retention times are in the trace's unit (minutes). NaN means "not determinable"; + a metric that needs a peak is NaN when no boundary was found. + """ + + apex_rt: float = float("nan") + apex_intensity: float = 0.0 + left_rt: float = float("nan") + right_rt: float = float("nan") + boundary_is_truncated: bool = False + """A boundary sits on an end of the trace: the extraction window cut the peak.""" + n_points: int = 0 + """Scans between the boundaries, inclusive.""" + area: float = 0.0 + """Trapezoidal area between the boundaries.""" + fwhm: float = float("nan") + fwhm_is_truncated: bool = False + points_across_peak: int = 0 + asymmetry: float = float("nan") + """(right half-max crossing - apex) / (apex - left half-max crossing). 1 is + symmetric; above 1 tails, below 1 fronts.""" + snr: float = float("nan") + """Apex intensity / MAD of the whole trace (MsQuality ``signalToNoiseRatio``).""" + prominence: float = float("nan") + """(apex - q10) / q10 over the whole trace (MsQuality ``peakProminence``).""" + beta_cor: float = float("nan") + beta_snr: float = float("nan") + + @property + def found(self) -> bool: + return self.apex_intensity > 0.0 and not math.isnan(self.left_rt) + + @property + def width(self) -> float: + return self.right_rt - self.left_rt + + def as_dict(self) -> Dict[str, object]: + d = {f: getattr(self, f) for f in self.__dataclass_fields__} + d["width"] = self.width + return d + + +def peak_metrics(xic: "XIC", *, + contains_rt: Optional[float] = None, + apex_window: Optional[Tuple[float, float]] = None, + threshold: float = 0.1, + baseline_threshold: float = 0.1, + baseline_quantile: float = 0.1, + nearest_valley: bool = False) -> PeakMetrics: + """Boundaries, in-boundary area and quality metrics of one peak in an XIC. + + The boundary arguments are those of :func:`peak_boundary`: ``contains_rt`` (or + ``apex_window``) selects which peak is measured, and ``nearest_valley=True`` is + recommended for zero-filled traces (see there). Metric definitions follow their R sources: + + * ``area``, ``n_points`` and the ``beta_*`` metrics use the points between the + boundaries (MsQuality ``gaussianSimilarity`` does the same); + * ``fwhm``, ``points_across_peak``, ``asymmetry`` use the contiguous half-maximum + region around the apex, crossings interpolated (as :attr:`XIC.fwhm`); + * ``snr`` and ``prominence`` use the whole trace, as MsQuality does by default, + with the APEX intensity as numerator. With neither selector the apex is the + trace maximum, so both equal MsQuality's values. Both are NaN when more than + half (MAD) or a tenth (q10) of the trace is zero, as MsQuality returns NA -- + common for a zero-filled XIC over a wide window. + + These are measurements, not verdicts: no threshold is applied here. + """ + b = peak_boundary(xic, contains_rt=contains_rt, apex_window=apex_window, threshold=threshold, + baseline_threshold=baseline_threshold, + baseline_quantile=baseline_quantile, + nearest_valley=nearest_valley) + if b is None: + return PeakMetrics() + left, apex, right = b + rt = np.asarray(xic.rt, dtype=float) + inten = np.asarray(xic.intensity, dtype=float) + apex_int = float(inten[apex]) + apex_rt = float(rt[apex]) + + seg_rt, seg_int = rt[left:right + 1], inten[left:right + 1] + area = float(_trapz(seg_int, seg_rt)) if right > left else 0.0 + + fwhm = asym = float("nan") + fwhm_trunc, pts = False, 0 + cr = xic._half_max_crossings(apex) + if cr is not None: + l_x, r_x, lo, hi = cr + fwhm = max(0.0, r_x - l_x) + fwhm_trunc = lo == 0 or hi == len(inten) - 1 + pts = int(hi - lo + 1) + if apex_rt - l_x > 0 and r_x - apex_rt > 0: + asym = (r_x - apex_rt) / (apex_rt - l_x) + + med = float(np.median(inten)) + mad = _MAD_CONSTANT * float(np.median(np.abs(inten - med))) + snr = apex_int / mad if mad > 0 else float("nan") + q10 = float(np.quantile(inten, baseline_quantile)) + prominence = (apex_int - q10) / q10 if q10 > 0 else float("nan") + + beta_cor, beta_snr = beta_values(seg_int, seg_rt) if len(seg_int) >= 5 else (float("nan"),) * 2 + + return PeakMetrics(apex_rt=apex_rt, apex_intensity=apex_int, + left_rt=float(rt[left]), right_rt=float(rt[right]), + boundary_is_truncated=(left == 0 or right == len(inten) - 1), + n_points=int(right - left + 1), area=area, + fwhm=fwhm, fwhm_is_truncated=fwhm_trunc, points_across_peak=pts, + asymmetry=asym, snr=snr, prominence=prominence, + beta_cor=beta_cor, beta_snr=beta_snr) + + # --------------------------------------------------------------------------- # Isotope co-elution -- trace-level precursor evidence # --------------------------------------------------------------------------- diff --git a/tests/data/r_peak_reference.json b/tests/data/r_peak_reference.json new file mode 100644 index 0000000..46a8fe9 --- /dev/null +++ b/tests/data/r_peak_reference.json @@ -0,0 +1,1449 @@ +{ + "provenance": "Reference values from R 4.5.2 running, verbatim: Chromatograms R/helpers.R `.peak_boundary_one` (main, 2026-09-08), MetaboCoreUtils R/peak-shape-quality.R `betaValues` + `.scale_zero_one` (main), and the MsQuality signalToNoiseRatio (max / mad) and peakProminence ((max - q10) / q10) formulas, whole trace. Beta 'in_boundary' uses rt >= left & rt <= right, as MsQuality gaussianSimilarity. null = R Inf (MAD or q10 is zero; MsQuality itself returns NA there). Built 2026-09-24; 'kumler_skinny' is the skinny_peak example from the betaValues docs.", + "traces": { + "gauss_floor": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0001, + 50.0003, + 50.0015, + 50.0066, + 50.027, + 50.1035, + 50.3727, + 51.2607, + 54.0065, + 61.9613, + 83.5463, + 138.3826, + 268.7491, + 558.6069, + 1160.8997, + 2329.4181, + 4443.6934, + 8005.9509, + 13583.5283, + 21676.5167, + 32515.2467, + 45833.3362, + 60703.066, + 75533.9602, + 88299.6903, + 96973.3234, + 100050.0, + 96973.3234, + 88299.6903, + 75533.9602, + 60703.066, + 45833.3362, + 32515.2467, + 21676.5167, + 13583.5283, + 8005.9509, + 4443.6934, + 2329.4181, + 1160.8997, + 558.6069, + 268.7491, + 138.3826, + 83.5463, + 61.9613, + 54.0065, + 51.2607, + 50.3727, + 50.1035, + 50.027, + 50.0066, + 50.0015, + 50.0003, + 50.0001, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0 + ], + "r": { + "left_rt": 0.0, + "right_rt": 7.7, + "beta_cor_in_boundary": 0.3516653011061723, + "beta_snr_in_boundary": 6.230102222277778, + "beta_cor_whole": 0.6307063740494046, + "beta_snr_whole": 6.329354278564223, + "snr": 224942668.2833049, + "prominence": 2000.0 + } + }, + "two_peaks": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 100.0, + 376.5071, + 117.1655, + 267.4668, + 212.1202, + 236.0442, + 287.8276, + 149.6807, + 198.7603, + 261.9909, + 391.8264, + 347.6176, + 157.7197, + 410.5128, + 378.3538, + 320.9276, + 596.8128, + 738.8656, + 1404.1555, + 2476.2112, + 4721.8066, + 8222.76, + 13831.7913, + 21852.3182, + 32605.9759, + 45982.4771, + 60901.6632, + 75584.0555, + 88398.262, + 97149.3815, + 100395.5246, + 97240.3459, + 88435.4007, + 75709.3733, + 60813.2453, + 46056.1124, + 32859.7029, + 22002.6873, + 13774.6113, + 8286.4124, + 4672.6897, + 2659.6979, + 1346.8269, + 673.2471, + 527.205, + 427.3231, + 219.0077, + 267.7344, + 151.5502, + 148.9947, + 226.3857, + 222.3704, + 547.0693, + 710.5343, + 1133.1619, + 2119.3477, + 3424.7819, + 5580.9088, + 8843.4092, + 13138.5743, + 18464.2773, + 24385.8853, + 30536.8532, + 35404.8773, + 38895.7462, + 40123.2581, + 39048.615, + 35472.3198, + 30410.0049, + 24459.6796, + 18523.6896, + 13238.308, + 8802.5868, + 5522.0779, + 3339.461, + 2035.9451, + 1058.0041, + 752.1873, + 358.6912, + 200.755, + 355.5169, + 242.8584, + 131.6442, + 184.8157, + 128.7266, + 213.0563, + 245.4843, + 207.3994, + 233.2446, + 378.4084, + 344.6245, + 142.4389, + 258.9162, + 321.209, + 334.2826, + 325.3333, + 378.3756, + 143.8745, + 116.7319, + 246.9796, + 179.2416 + ], + "r": { + "left_rt": 1.5, + "right_rt": 4.6, + "beta_cor_in_boundary": 0.9746864784421962, + "beta_snr_in_boundary": 6.490369936824226, + "beta_cor_whole": 0.5200140647352194, + "beta_snr_whole": 6.288183995375695, + "snr": 116.97531817439128, + "prominence": 635.5439738979976 + } + }, + "zero_filled": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1432.6729, + 3375.9768, + 7331.6084, + 14673.9303, + 27067.0566, + 46013.2598, + 72089.5577, + 104090.0242, + 138513.8648, + 169873.1633, + 192001.0883, + 200000.0, + 192001.0883, + 169873.1633, + 138513.8648, + 104090.0242, + 72089.5577, + 46013.2598, + 27067.0566, + 14673.9303, + 7331.6084, + 3375.9768, + 1432.6729, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "r": { + "left_rt": 0.0, + "right_rt": 6.2, + "beta_cor_in_boundary": 0.9985572643492612, + "beta_snr_in_boundary": 7.144119225855359, + "beta_cor_whole": 0.9985572643492612, + "beta_snr_whole": 7.144119225855359, + "snr": null, + "prominence": null + } + }, + "tailing": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 200.0, + 476.5071, + 217.1655, + 367.4668, + 312.1201, + 336.0439, + 387.826, + 249.674, + 298.7333, + 361.8874, + 491.4537, + 446.3569, + 253.7131, + 498.5515, + 444.8076, + 332.5449, + 478.0636, + 330.2587, + 393.2559, + 296.7931, + 428.1132, + 366.8092, + 398.2629, + 325.8015, + 240.7292, + 299.1409, + 348.5972, + 200.0953, + 248.5717, + 326.058, + 495.5246, + 417.0225, + 285.7105, + 325.4131, + 260.1794, + 372.7762, + 494.4562, + 476.1706, + 341.083, + 430.4615, + 378.9962, + 43421.0678, + 57687.0372, + 59692.1083, + 56847.0961, + 52055.443, + 46656.4429, + 41648.3736, + 36853.0929, + 32619.0889, + 28904.6781, + 25502.2065, + 22749.1304, + 20053.9439, + 17697.2508, + 15796.6704, + 13875.6144, + 12210.6513, + 10832.6602, + 9553.8952, + 8459.4294, + 7468.6286, + 6836.0525, + 5846.6139, + 5205.1231, + 4616.9513, + 4256.7063, + 3694.2555, + 3336.1591, + 2963.3629, + 2662.1297, + 2427.6431, + 2083.544, + 1825.0161, + 1683.504, + 1637.282, + 1357.1365, + 1388.193, + 1120.418, + 976.7648, + 1093.9585, + 924.0616, + 751.6115, + 746.305, + 636.8995, + 673.5636, + 663.721, + 588.268, + 581.1171, + 697.1569, + 637.6698, + 412.8008, + 509.2601, + 553.8871, + 551.3706, + 528.6631, + 569.5638, + 324.3478, + 287.7493, + 409.6522, + 334.5501 + ], + "r": { + "left_rt": 4.0, + "right_rt": 7.6, + "beta_cor_in_boundary": 0.31545214399012944, + "beta_snr_in_boundary": 5.598402916729488, + "beta_cor_whole": 0.6579881760302444, + "beta_snr_whole": 5.872167748144791, + "snr": 142.86623102085298, + "prominence": 200.12363899295502 + } + }, + "high_valley": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 100.0, + 376.5071, + 117.1655, + 267.4668, + 212.1201, + 236.0439, + 287.826, + 149.674, + 198.7333, + 261.8874, + 391.4537, + 346.3569, + 153.7131, + 398.5516, + 344.8077, + 232.5454, + 378.0652, + 230.2636, + 293.2714, + 196.8396, + 328.2476, + 267.1818, + 299.2559, + 228.3434, + 146.9814, + 213.9162, + 282.1448, + 173.28, + 301.9667, + 534.9715, + 993.2478, + 1428.2575, + 2170.7136, + 3632.4463, + 5779.2826, + 9178.2353, + 13958.1761, + 20231.9027, + 28182.856, + 38139.5558, + 49492.0639, + 62033.1549, + 74636.5269, + 86755.9214, + 97672.1778, + 106361.6044, + 112365.1326, + 116084.9207, + 117477.7508, + 117439.3355, + 116613.2667, + 115375.5926, + 114372.2385, + 112971.5968, + 111105.5944, + 108369.6205, + 103775.9298, + 97277.5398, + 88886.75, + 78731.5401, + 67488.4656, + 55823.3179, + 44748.5744, + 34191.8716, + 25303.1549, + 18007.3213, + 12493.0074, + 8190.1641, + 5274.8016, + 3265.2661, + 1997.0466, + 1252.3916, + 689.9566, + 386.6976, + 295.1391, + 344.335, + 176.43, + 321.1253, + 160.8755, + 115.5432, + 321.0575, + 229.7753, + 126.9806, + 183.255, + 128.2363, + 212.9117, + 245.4443, + 207.389, + 233.242, + 378.4078, + 344.6244, + 142.4389, + 258.9162, + 321.209, + 334.2826, + 325.3333, + 378.3756, + 143.8745, + 116.7319, + 246.9796, + 179.2416 + ], + "r": { + "left_rt": 2.7, + "right_rt": 7.4, + "beta_cor_in_boundary": 0.9959606342756425, + "beta_snr_in_boundary": 6.894747360854513, + "beta_cor_whole": 0.8360245359461937, + "beta_snr_whole": 6.573666974893433, + "snr": 340.61286921392013, + "prominence": 763.2663559579503 + } + }, + "noisy_flat": { + "rt": [ + 0.0, + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, + 1.3, + 1.4, + 1.5, + 1.6, + 1.7, + 1.8, + 1.9, + 2.0, + 2.1, + 2.2, + 2.3, + 2.4, + 2.5, + 2.6, + 2.7, + 2.8, + 2.9, + 3.0, + 3.1, + 3.2, + 3.3, + 3.4, + 3.5, + 3.6, + 3.7, + 3.8, + 3.9, + 4.0, + 4.1, + 4.2, + 4.3, + 4.4, + 4.5, + 4.6, + 4.7, + 4.8, + 4.9, + 5.0, + 5.1, + 5.2, + 5.3, + 5.4, + 5.5, + 5.6, + 5.7, + 5.8, + 5.9, + 6.0, + 6.1, + 6.2, + 6.3, + 6.4, + 6.5, + 6.6, + 6.7, + 6.8, + 6.9, + 7.0, + 7.1, + 7.2, + 7.3, + 7.4, + 7.5, + 7.6, + 7.7, + 7.8, + 7.9, + 8.0, + 8.1, + 8.2, + 8.3, + 8.4, + 8.5, + 8.6, + 8.7, + 8.8, + 8.9, + 9.0, + 9.1, + 9.2, + 9.3, + 9.4, + 9.5, + 9.6, + 9.7, + 9.8, + 9.9, + 10.0 + ], + "intensity": [ + 1000.0, + 1276.5071, + 1017.1655, + 1167.4668, + 1112.1201, + 1136.0439, + 1187.826, + 1049.674, + 1098.7333, + 1161.8874, + 1291.4537, + 1246.3569, + 1053.7131, + 1298.5515, + 1244.8076, + 1132.5449, + 1278.0636, + 1130.2587, + 1193.2559, + 1096.7931, + 1228.1132, + 1166.8092, + 1198.2629, + 1125.8015, + 1040.7292, + 1099.1409, + 1148.5972, + 1000.0953, + 1048.5717, + 1126.058, + 1295.5246, + 1217.0225, + 1085.7105, + 1125.4131, + 1060.1794, + 1172.7762, + 1294.4562, + 1276.1706, + 1141.083, + 1230.4615, + 1178.9962, + 1280.2792, + 1135.9246, + 1064.6293, + 1208.4145, + 1238.7914, + 1084.9571, + 1154.1705, + 1042.7591, + 1034.3155, + 1090.66, + 1034.7673, + 1243.5995, + 1166.1679, + 1121.3931, + 1261.87, + 1142.4015, + 1067.4974, + 1092.8025, + 1052.4756, + 1050.9428, + 1024.6589, + 1243.2691, + 1005.0012, + 1026.4168, + 1023.2581, + 1179.2856, + 1072.4437, + 1116.4208, + 1098.4532, + 1110.3551, + 1152.2093, + 1051.9801, + 1008.6666, + 1057.0806, + 1178.4678, + 1046.2368, + 1207.8275, + 1055.2485, + 1013.2554, + 1220.1638, + 1129.4399, + 1026.8597, + 1083.2131, + 1028.2223, + 1112.9072, + 1145.4429, + 1107.3886, + 1133.2419, + 1278.4078, + 1244.6243, + 1042.4389, + 1158.9162, + 1221.209, + 1234.2826, + 1225.3333, + 1278.3756, + 1043.8745, + 1016.7319, + 1146.9796, + 1079.2416 + ], + "r": { + "left_rt": 1.2, + "right_rt": 2.4, + "beta_cor_in_boundary": 0.31510082643953274, + "beta_snr_in_boundary": 3.4081140315106215, + "beta_cor_whole": 0.14648907492670468, + "beta_snr_whole": 3.536847106419052, + "snr": 11.566063632088909, + "prominence": 0.2645851229724957 + } + }, + "kumler_skinny": { + "rt": [ + 101.0, + 102.0, + 103.0, + 104.0, + 105.0, + 106.0, + 107.0, + 108.0, + 109.0, + 110.0, + 111.0, + 112.0, + 113.0, + 114.0, + 115.0, + 116.0, + 117.0, + 118.0, + 119.0, + 120.0, + 121.0, + 122.0, + 123.0, + 124.0, + 125.0, + 126.0, + 127.0, + 128.0, + 129.0, + 130.0, + 131.0, + 132.0, + 133.0, + 134.0, + 135.0, + 136.0, + 137.0, + 138.0, + 139.0, + 140.0, + 141.0, + 142.0, + 143.0, + 144.0, + 145.0, + 146.0, + 147.0, + 148.0, + 149.0, + 150.0, + 151.0, + 152.0, + 153.0, + 154.0, + 155.0, + 156.0, + 157.0, + 158.0, + 159.0, + 160.0 + ], + "intensity": [ + 9107.0, + 3326.0, + 9523.0, + 3245.0, + 3429.0, + 9394.0, + 1123.0, + 935.0, + 5128.0, + 8576.0, + 2711.0, + 3427.0, + 7294.0, + 8109.0, + 9288.0, + 6997.0, + 9756.0, + 8034.0, + 1317.0, + 8866.0, + 13877.0, + 14854.0, + 28296.0, + 57101.0, + 92209.0, + 151797.0, + 222386.0, + 299402.0, + 365045.0, + 394255.0, + 402680.0, + 363996.0, + 293985.0, + 222989.0, + 147007.0, + 94947.0, + 52924.0, + 32438.0, + 11511.0, + 10836.0, + 8046.0, + 601.0, + 889.0, + 5917.0, + 2690.0, + 5381.0, + 9901.0, + 8494.0, + 3349.0, + 8283.0, + 3410.0, + 5935.0, + 3332.0, + 7041.0, + 3284.0, + 7478.0, + 76.0, + 3739.0, + 2158.0, + 5507.0 + ], + "r": { + "left_rt": 119.0, + "right_rt": 142.0, + "beta_cor_in_boundary": 0.9885576333796767, + "beta_snr_in_boundary": 7.070372735943331, + "beta_cor_whole": 0.727436991527969, + "beta_snr_whole": 6.789728356498464, + "snr": 55.805206292589965, + "prominence": 193.1655817541829 + } + } + } +} \ No newline at end of file diff --git a/tests/test_chrom_extract.py b/tests/test_chrom_extract.py new file mode 100644 index 0000000..3d2860b --- /dev/null +++ b/tests/test_chrom_extract.py @@ -0,0 +1,157 @@ +"""Tests for chrom_extract: peak table -> chromatograms across runs.""" + +import math + +import numpy as np +import pytest + +from mzml_utils import Spectrum +from mzml_utils.chrom_extract import (extract_chromatograms, chromatogram_rows, + source_key, ExtractedChromatogram) +from mzml_utils.xic import extract_xics + + +def _ms1(scan, rt, peaks): + mz = np.array([p[0] for p in peaks], dtype=float) + inten = np.array([p[1] for p in peaks], dtype=float) + return Spectrum(scan_num=scan, mz=mz, intensity=inten, ms_level=1, rt=rt, + tic=float(inten.sum())) + + +def _ms2(scan, rt, peaks): + mz = np.array([p[0] for p in peaks], dtype=float) + inten = np.array([p[1] for p in peaks], dtype=float) + return Spectrum(scan_num=scan, mz=mz, intensity=inten, ms_level=2, rt=rt, + precursor_mz=1000.5, activation_type="HCD") + + +class FakeReader: + def __init__(self, specs): + self._specs = specs + + def iter_spectra(self): + yield from self._specs + + +def _run(shift=0.0, scale=1.0): + """20 MS1 scans 1 min apart. Species A (1000.5000) elutes at 5, B (800.2500) at 14, + plus a neighbour 0.02 Da from A that must stay outside a 10 ppm window.""" + specs = [] + for i in range(20): + rt = float(i) + shift + a = 1e5 * scale * math.exp(-0.5 * ((i - 5) / 1.2) ** 2) + b = 5e4 * scale * math.exp(-0.5 * ((i - 14) / 1.0) ** 2) + peaks = sorted([(1000.5000, a), (1000.5200, 7e3), (800.2500, b), (650.0, 100.0)]) + specs.append(_ms1(2 * i + 1, rt, peaks)) + specs.append(_ms2(2 * i + 2, rt + 0.5, [(204.087, 9e3)])) + return FakeReader(specs) + + +TABLE = [ + {"id": "A", "mz": 1000.5, "rt_min": 2.0, "rt_max": 9.0, "run": "r1", "note": "kept"}, + {"id": "B", "mz": 800.25, "rt_min": 10.0, "rt_max": 18.0, "run": "r2", "note": "kept too"}, +] + + +def test_source_key(): + assert source_key("/x/spectra_cache/run_01.spectra.db") == "run_01" + assert source_key("/x/run_01.mzML") == "run_01" + assert source_key("C:/x/run_01.mzml") == "run_01" + + +def test_per_row_windows_and_carry_through(): + res = extract_chromatograms(TABLE, {"r1": _run()}, metrics=True) + assert [r.row["id"] for r in res] == ["A", "B"] + a, b = res + assert a.xic.rt.min() >= 2.0 and a.xic.rt.max() <= 9.0 and a.xic.n_points == 8 + assert b.xic.rt.min() >= 10.0 and b.xic.n_points == 9 + assert a.metrics.apex_rt == 5.0 and b.metrics.apex_rt == 14.0 + assert a.row == TABLE[0] # carried through unchanged + assert a.xic.intensity.max() == pytest.approx(1e5) # the 0.02 Da neighbour is excluded + + +def test_by_restricts_rows_to_their_run(): + srcs = {"r1": _run(), "r2": _run(shift=0.3)} + res = extract_chromatograms(TABLE, srcs, by="run") + assert [(r.row["id"], r.source) for r in res] == [("A", "r1"), ("B", "r2")] + assert res[1].metrics.apex_rt == pytest.approx(14.3) + + +def test_without_by_every_row_from_every_run_for_overlays(): + srcs = {"r1": _run(), "r2": _run(shift=0.4, scale=0.5)} + res = extract_chromatograms(TABLE[:1], srcs) + assert [r.source for r in res] == ["r1", "r2"] + assert res[1].metrics.apex_rt - res[0].metrics.apex_rt == pytest.approx(0.4) + assert res[1].metrics.apex_intensity == pytest.approx(0.5e5) + + +def test_single_row_matches_extract_xics(): + reader = _run() + res = extract_chromatograms(TABLE[:1], {"r1": reader}, tolerance=10.0, metrics=False) + ref = extract_xics(reader, {"A": 1000.5}, tolerance=10.0, ms_level=1, + rt_range=(2.0, 9.0), collect_tic=False)["A"] + np.testing.assert_array_equal(res[0].xic.rt, ref.rt) + np.testing.assert_allclose(res[0].xic.intensity, ref.intensity, rtol=1e-12, atol=0) + assert res[0].metrics is None + + +def test_unsorted_mz_is_handled(): + specs = [_ms1(1, 1.0, [(900.0, 1.0), (500.0, 7.0), (1000.5, 3.0)]), + _ms1(2, 2.0, [(1000.5, 4.0), (100.0, 9.0)])] + res = extract_chromatograms([{"mz": 1000.5, "rt_min": 0, "rt_max": 3}], + {"x": FakeReader(specs)}, metrics=False) + np.testing.assert_allclose(res[0].xic.intensity, [3.0, 4.0]) + + +def test_mz_min_max_override_tolerance_and_da_unit(): + row = {"mz_min": 1000.49, "mz_max": 1000.53, "rt_min": 4.0, "rt_max": 6.0} + res = extract_chromatograms([row], {"r": _run()}, metrics=False) + assert res[0].xic.intensity[1] == pytest.approx(1e5 + 7e3) # both peaks inside + assert res[0].xic.target_mz == pytest.approx(1000.51) + res = extract_chromatograms([{"mz": 1000.5, "rt_min": 4.0, "rt_max": 6.0}], {"r": _run()}, + tolerance=0.03, unit="Da", metrics=False) + assert res[0].xic.intensity[1] == pytest.approx(1e5 + 7e3) + + +def test_apex_rt_column_measures_the_peak_containing_the_ms2(): + specs = [] + for i in range(30): + y = 1e5 * math.exp(-0.5 * ((i - 6) / 1.0) ** 2) + 3e4 * math.exp(-0.5 * ((i - 20) / 1.0) ** 2) + specs.append(_ms1(i + 1, float(i), [(1000.5, y + 10.0 + 0.01 * i)])) + row = {"mz": 1000.5, "rt_min": 0.0, "rt_max": 29.0, "ms2_rt": 21.6} + res = extract_chromatograms([row], {"r": FakeReader(specs)}, apex_rt="ms2_rt") + assert res[0].metrics.apex_rt == 20.0 + res = extract_chromatograms([row], {"r": FakeReader(specs)}) + assert res[0].metrics.apex_rt == 6.0 + row["ms2_rt"] = float("nan") # missing MS2 time falls back to the tallest + res = extract_chromatograms([row], {"r": FakeReader(specs)}, apex_rt="ms2_rt") + assert res[0].metrics.apex_rt == 6.0 + + +def test_row_with_no_scans_gets_empty_xic(): + res = extract_chromatograms([{"mz": 1000.5, "rt_min": 100.0, "rt_max": 101.0}], {"r": _run()}) + assert res[0].xic.n_points == 0 and not res[0].metrics.found + + +def test_validation_errors(): + with pytest.raises(ValueError, match="rt_min"): + extract_chromatograms([{"mz": 1.0, "rt_max": 2.0}], {"r": _run()}) + with pytest.raises(ValueError, match="NaN"): + extract_chromatograms([{"mz": 1.0, "rt_min": float("nan"), "rt_max": 2.0}], {"r": _run()}) + with pytest.raises(ValueError, match="match no source"): + extract_chromatograms(TABLE, {"r1": _run()}, by="run") + with pytest.raises(ValueError, match="both"): + extract_chromatograms([{"mz_min": 1.0, "rt_min": 0, "rt_max": 2}], {"r": _run()}) + with pytest.raises(ValueError, match="share a run name"): + extract_chromatograms(TABLE, ["/a/x.mzML", "/b/x.mzML"]) + + +def test_chromatogram_rows_flattens_and_warns_on_collision(): + res = extract_chromatograms(TABLE[:1], {"r1": _run()}) + rows = chromatogram_rows(res) + assert rows[0]["id"] == "A" and rows[0]["source"] == "r1" and rows[0]["apex_rt"] == 5.0 + assert {"n_scans", "area", "beta_cor", "width"} <= set(rows[0]) + clash = [dict(TABLE[0], area=-1)] + with pytest.warns(UserWarning, match="area"): + rows = chromatogram_rows(extract_chromatograms(clash, {"r1": _run()})) + assert rows[0]["area"] > 0 diff --git a/tests/test_xic.py b/tests/test_xic.py index 2106dfc..8c38460 100644 --- a/tests/test_xic.py +++ b/tests/test_xic.py @@ -338,3 +338,160 @@ def test_missing_precursor(self): assert not ev.found assert ev.pattern_cosine == 0.0 assert not ev.m_minus_1_coelutes + + +# --------------------------------------------------------------------------- +# Peak boundaries and quality metrics (ports of Chromatograms / MsQuality / +# MetaboCoreUtils; tests/data/r_peak_reference.json holds the R values) +# --------------------------------------------------------------------------- + +import json +import math +from pathlib import Path + +from mzml_utils.xic import peak_boundary, peak_metrics, beta_values, PeakMetrics + +_R_REF = json.loads((Path(__file__).parent / "data" / "r_peak_reference.json").read_text()) + + +def _gauss(rt, c, s, h): + return h * np.exp(-0.5 * ((rt - c) / s) ** 2) + + +class TestRReference: + """The port reproduces the R functions run verbatim on the same traces.""" + + @pytest.mark.parametrize("name", sorted(_R_REF["traces"])) + def test_matches_r(self, name): + tr = _R_REF["traces"][name] + rt, y, ref = np.array(tr["rt"]), np.array(tr["intensity"]), tr["r"] + xic = XIC(target_mz=0.0, rt=rt, intensity=y) + left, _apex, right = peak_boundary(xic) + assert rt[left] == ref["left_rt"] and rt[right] == ref["right_rt"] + seg = slice(left, right + 1) + got = dict(zip(("beta_cor_in_boundary", "beta_snr_in_boundary"), beta_values(y[seg], rt[seg]))) + got.update(zip(("beta_cor_whole", "beta_snr_whole"), beta_values(y, rt))) + pm = peak_metrics(xic) + got["snr"], got["prominence"] = pm.snr, pm.prominence + assert pm.beta_cor == pytest.approx(ref["beta_cor_in_boundary"], rel=1e-12) + for k, v in got.items(): + if ref[k] is None: # R Inf: MAD or q10 is zero; MsQuality returns NA -> NaN here + assert math.isnan(v), k + else: + assert v == pytest.approx(ref[k], rel=1e-12), k + + +class TestPeakBoundary: + def test_valley_between_two_peaks(self): + rt = np.linspace(0, 10, 101) + y = _gauss(rt, 3, 0.4, 1e5) + _gauss(rt, 6.5, 0.4, 4e4) + 100.0 + left, apex, right = peak_boundary(XIC(0.0, rt, y)) + assert rt[apex] == pytest.approx(3.0) + assert 4.0 < rt[right] < 6.0 # stops in the valley, not at the second peak + + def test_apex_window_selects_the_smaller_peak(self): + rt = np.linspace(0, 10, 101) + y = _gauss(rt, 3, 0.4, 1e5) + _gauss(rt, 6.5, 0.4, 4e4) + 100.0 + 0.01 * rt + left, apex, right = peak_boundary(XIC(0.0, rt, y), apex_window=(6.0, 7.0)) + assert rt[apex] == pytest.approx(6.5) + assert 4.0 < rt[left] < 6.0 + + def test_contains_rt_takes_the_peak_whose_tail_holds_the_ms2(self): + # MS2 triggered on the tail, 1 min after the apex: the containing peak is the + # big one, not the tallest point near the MS2 time. + rt = np.linspace(60, 72, 121) + y = (_gauss(rt, 64, 0.8, 3e4) + _gauss(rt, 67.1, 0.25, 1.8e8) + + np.where(rt > 67.1, 1.5e7 * np.exp(-(rt - 67.1) / 0.6), 0.0)) + xic = XIC(0.0, rt, y) + left, apex, right = peak_boundary(xic, contains_rt=68.1, nearest_valley=True) + assert rt[apex] == pytest.approx(67.1) and rt[left] <= 68.1 <= rt[right] + w_left, w_apex, _ = peak_boundary(xic, apex_window=(67.6, 68.6)) + assert rt[w_apex] > 67.5 # the window picks a tail point + + def test_contains_rt_peels_to_the_smaller_peak(self): + rt = np.linspace(0, 10, 101) + y = _gauss(rt, 3, 0.4, 1e5) + _gauss(rt, 6.5, 0.4, 4e4) + 100.0 + 0.01 * rt + xic = XIC(0.0, rt, y) + _l, apex, _r = peak_boundary(xic, contains_rt=6.8) + assert rt[apex] == pytest.approx(6.5) + assert peak_boundary(xic, contains_rt=50.0) is None # outside the trace + with pytest.raises(ValueError): + peak_boundary(xic, contains_rt=6.8, apex_window=(6, 7)) + + def test_high_valley_falls_back_to_threshold_on_both_sides(self): + rt = np.linspace(0, 10, 101) + y = _gauss(rt, 4.6, 0.5, 1e5) + _gauss(rt, 5.6, 0.5, 9e4) + 100.0 + 0.01 * rt + left, apex, right = peak_boundary(XIC(0.0, rt, y)) + # The valley between the two peaks is far above baseline + 10 %, so both + # boundaries come from the 10 % threshold and enclose both peaks. + assert rt[left] < 4.0 and rt[right] > 6.0 + + def test_zero_filled_left_boundary_runs_to_trace_start_in_r_mode(self): + rt = np.linspace(0, 10, 101) + y = np.where(np.abs(rt - 5) < 1.2, _gauss(rt, 5, 0.35, 2e5), 0.0) + xic = XIC(0.0, rt, y) + left, _apex, right = peak_boundary(xic) + assert left == 0 # R behaviour, documented + l2, _a2, r2 = peak_boundary(xic, nearest_valley=True) + # both boundaries are the first zero next to the peak + assert y[l2] == 0 and y[l2 + 1] > 0 and y[r2] == 0 and y[r2 - 1] > 0 + assert r2 == right # the right side is unchanged + + def test_nearest_valley_on_noisy_baseline_matches_r(self): + tr = _R_REF["traces"]["two_peaks"] + xic = XIC(0.0, np.array(tr["rt"]), np.array(tr["intensity"])) + assert peak_boundary(xic) == peak_boundary(xic, nearest_valley=True) + + def test_flat_top_is_not_mistaken_for_a_valley(self): + rt = np.arange(9, dtype=float) + y = np.array([0, 1, 5, 9, 9, 9, 5, 1, 0], dtype=float) + assert peak_boundary(XIC(0.0, rt, y), nearest_valley=True) == (0, 3, 8) + + def test_degenerate_inputs(self): + assert peak_boundary(XIC(0.0, np.array([1.0, 2.0]), np.array([5.0, 3.0]))) is None + assert peak_boundary(XIC(0.0, np.arange(5.0), np.zeros(5))) is None + xic = XIC(0.0, np.arange(5.0), np.array([0, 1, 5, 1, 0.0])) + assert peak_boundary(xic, apex_window=(10.0, 11.0)) is None + with pytest.raises(ValueError): + peak_boundary(XIC(0.0, np.arange(5.0), np.array([0, 1, np.nan, 1, 0.0]))) + + +class TestPeakMetrics: + def test_area_is_within_boundaries_not_whole_trace(self): + rt = np.linspace(0, 20, 201) + y = _gauss(rt, 5, 0.3, 1e5) + _gauss(rt, 15, 0.3, 1e5) + xic = XIC(0.0, rt, y) + pm = peak_metrics(xic, nearest_valley=True) + one_peak = 1e5 * 0.3 * math.sqrt(2 * math.pi) + assert pm.area == pytest.approx(one_peak, rel=1e-3) + assert xic.area == pytest.approx(2 * one_peak, rel=1e-3) + + def test_asymmetry_symmetric_and_tailing(self): + rt = np.linspace(0, 10, 1001) + sym = peak_metrics(XIC(0.0, rt, _gauss(rt, 5, 0.3, 1e5) + 10.0)) + assert sym.asymmetry == pytest.approx(1.0, abs=0.02) + # exponentially modified Gaussian-like tail on the right + tail = np.where(rt >= 4, np.exp(-(rt - 4) / 0.8) * (1 - np.exp(-(rt - 4) / 0.15)), 0.0) + pm = peak_metrics(XIC(0.0, rt, 1e5 * tail + 10.0)) + assert pm.asymmetry > 1.5 + + def test_fwhm_matches_xic_property_for_the_tallest_peak(self): + rt = np.linspace(0, 10, 101) + xic = XIC(0.0, rt, _gauss(rt, 5, 0.4, 1e5) + 10.0) + pm = peak_metrics(xic) + assert pm.fwhm == pytest.approx(xic.fwhm) + assert pm.fwhm == pytest.approx(2.3548 * 0.4, rel=0.02) + assert pm.points_across_peak == xic.points_across_peak + + def test_beta_cor_high_for_a_peak_low_for_noise(self): + rt = np.linspace(0, 1, 41) + peak = np.where((rt > 0) & (rt < 1), rt ** 2 * (1 - rt) ** 4, 0.0) * 1e6 + 1.0 + assert beta_values(peak, rt)[0] > 0.99 + noise = _R_REF["traces"]["noisy_flat"] + assert beta_values(np.array(noise["intensity"]))[0] < 0.5 + + def test_no_peak_gives_empty_metrics(self): + pm = peak_metrics(XIC(0.0, np.array([]), np.array([]))) + assert isinstance(pm, PeakMetrics) and not pm.found and math.isnan(pm.apex_rt) + d = pm.as_dict() + assert "width" in d and "beta_cor" in d