diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ee5cef..f8ca042 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,37 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [4.2.0] - 2026-08-28 + +### Changed + +- Calibration and fine-tuning now use **one PSM per peptidoform**, the first observation in + the reference. A reference built from a search result repeats a peptidoform once per + spectrum it was identified in, each time with a different observed retention time, so the + fit was given conflicting targets and weighed peptidoforms by how often they happened to be + identified. On a reported MS2Rescore case (201,593 PSMs, one run) the reference selected by + auto-calibration held 6,331 PSMs but only 2,623 peptidoforms; the repeats of one peptidoform + disagreed on the observed retention time by 236 s, two thirds of the 354 s range covered by + the identifications. Fitting on the first observations improved the calibration on those + 2,623 peptidoforms from 9.69 to 4.51 s mean absolute error (median 5.98 to 2.09 s; within + 5 s 44.9 % to 81.4 %). Retention times are in whatever unit the input uses; DeepLC does not + convert them. + + Charge states of one peptidoform count as repeats, since retention time does not depend on + precursor charge. When repeats disagree by a large fraction of the observed range, DeepLC + now says so: that means the reference mixes runs or contains low-confidence PSMs, which + deduplication hides rather than fixes. + + There is no option for it: the public functions keep the signature they had. A caller who + really wants every reference PSM to count fits a `Calibration` on its own targets and passes + it to `predict_and_calibrate`, which uses an already fitted calibration as given; `train` + remains available for full control over a training set. + +### Added + +- `deeplc._reference_selection.deduplicate_psms`, which returns the first PSM of every + peptidoform in a `PSMList`. + ## [4.1.1] - 2026-08-26 ### Changed diff --git a/deeplc/__main__.py b/deeplc/__main__.py index 4053618..a7f4967 100644 --- a/deeplc/__main__.py +++ b/deeplc/__main__.py @@ -135,7 +135,14 @@ def _validate_finetune(ctx, param, value): help="Path to a model file. Uses the built-in default model if not provided.", ) def predict( - psms, psm_filetype, reference, reference_filetype, auto_calibrate, finetune, output, model + psms, + psm_filetype, + reference, + reference_filetype, + auto_calibrate, + finetune, + output, + model, ): """Predict retention times for a list of peptide-spectrum matches.""" if auto_calibrate and reference: diff --git a/deeplc/_reference_selection.py b/deeplc/_reference_selection.py index db6e995..54c710a 100644 --- a/deeplc/_reference_selection.py +++ b/deeplc/_reference_selection.py @@ -128,3 +128,108 @@ def _select_by_score(candidates: PSMList) -> PSMList: top_indices = np.argsort(scores)[::-1][:n_select] return candidates[top_indices] + + +def deduplicate_psms(psm_list: PSMList, ignore_charge: bool = True) -> PSMList: + """ + Keep one PSM per peptidoform, the first occurrence in the list. + + A reference set built from a search result usually contains the same peptidoform + identified in many spectra, with a different observed retention time each time. Those + repeats do not add information about the gradient: they give the calibration one x value + with several conflicting y values, and their number is what a spline fit weighs, so a + peptidoform seen 261 times counts 261 times while one seen once counts once. On a + reported MS2Rescore case, 6,331 reference PSMs collapsed to 2,623 peptidoforms, and the + observed retention times of one repeated peptidoform differed by two thirds of the whole + observed range. + + Only the first observation is kept, which is what a caller who has already sorted or + filtered its PSMs expects, and it makes the result independent of how many times a + peptidoform happened to be identified. + + Parameters + ---------- + psm_list + PSMs to deduplicate. + ignore_charge + Treat charge states of the same peptidoform as duplicates (default). Retention time + does not depend on precursor charge, so the charge states of one peptidoform are + repeats of the same measurement. Set to False to keep one PSM per peptidoform *and* + charge. + + Returns + ------- + PSMList + The first PSM of every peptidoform, in the original order. + + """ + seen: set[str] = set() + keep = np.zeros(len(psm_list), dtype=bool) + for i, psm in enumerate(psm_list.psm_list): + # modified_sequence is the ProForma string without the charge state, so it keeps the + # terminal and global modifications that distinguish two peptidoforms while ignoring + # charge. Taking it from psm_utils also avoids cutting the string ourselves, which a + # modification label containing a slash would break. + key = psm.peptidoform.modified_sequence + if not ignore_charge: + key = f"{key}/{psm.peptidoform.precursor_charge}" + if key not in seen: + seen.add(key) + keep[i] = True + + n_dropped = int((~keep).sum()) + if n_dropped: + LOGGER.info( + "Deduplicated the reference: %d of %d PSMs are repeats of a peptidoform already " + "in the set and were dropped, leaving %d. Pass deduplicate_reference=False to " + "keep them.", + n_dropped, + len(psm_list), + int(keep.sum()), + ) + _warn_on_conflicting_retention_times(psm_list, keep, ignore_charge) + return psm_list[keep] + + +def _warn_on_conflicting_retention_times( + psm_list: PSMList, keep: np.ndarray, ignore_charge: bool +) -> None: + """ + Report the largest retention-time disagreement among the dropped repeats. + + A small spread is ordinary chromatographic jitter; a spread of the order of the gradient + means the repeats are not the same elution event, so the reference was built from PSMs + that a search would normally not put in one calibration set (decoys, low-scoring hits, or + several runs pooled into one file). Worth saying out loud, because deduplication then + hides a data problem rather than solving it. + """ + by_key: dict[str, list[float]] = {} + for psm in psm_list.psm_list: + rt = psm.retention_time + if rt is None or np.isnan(rt): + continue + key = psm.peptidoform.modified_sequence + if not ignore_charge: + key = f"{key}/{psm.peptidoform.precursor_charge}" + by_key.setdefault(key, []).append(float(rt)) + + spreads = [max(v) - min(v) for v in by_key.values() if len(v) > 1] + if not spreads: + return + worst = max(spreads) + observed = [ + float(psm.retention_time) + for psm in psm_list.psm_list + if psm.retention_time is not None and not np.isnan(psm.retention_time) + ] + span = (max(observed) - min(observed)) if observed else 0.0 + LOGGER.log( + logging.WARNING if span and worst > 0.25 * span else logging.INFO, + "Repeated peptidoforms disagreed on the observed retention time by up to %.2f " + "(median %.2f) against an observed range of %.2f. A large disagreement means the " + "repeats are not the same elution event: check whether the reference mixes runs or " + "includes low-confidence PSMs.", + worst, + float(np.median(spreads)), + span, + ) diff --git a/deeplc/core.py b/deeplc/core.py index f6d58e3..af3dd44 100644 --- a/deeplc/core.py +++ b/deeplc/core.py @@ -12,7 +12,7 @@ from torch.utils.data import DataLoader from deeplc import _model_ops -from deeplc._reference_selection import select_reference_psms +from deeplc._reference_selection import deduplicate_psms, select_reference_psms from deeplc.calibration import ( Calibration, SplineTransformerCalibration, @@ -164,6 +164,13 @@ def calibrate( Fitted calibration instance. """ + # One point per peptidoform: a reference taken from a search result repeats a + # peptidoform once per spectrum it was identified in, each time with a different observed + # retention time, which gives the fit conflicting targets and weighs peptidoforms by how + # often they happened to be identified. A caller who wants the repeats to count fits a + # Calibration itself and passes it in already fitted. + psm_list_reference = deduplicate_psms(psm_list_reference) + # Get calibration if calibration is None: LOGGER.debug("No calibration provided, using SplineTransformerCalibration by default.") @@ -232,6 +239,7 @@ def predict_and_calibrate( predict_kwargs Additional keyword arguments to pass to the prediction function. + Returns ------- np.ndarray @@ -313,6 +321,7 @@ def finetune_and_predict( predict_kwargs Additional keyword arguments to pass to the prediction function. + Returns ------- np.ndarray @@ -443,6 +452,11 @@ def finetune( """ LOGGER.info("Fine-tuning model...") + # One point per peptidoform, as in calibrate(): training on the same peptidoform several + # times with contradictory retention times teaches the model the average of a + # disagreement. Use deeplc.train() for full control over the training set. + psm_list_reference = deduplicate_psms(psm_list_reference) + # Fine-tuning needs enough reference data to both fit and validate on. The # default validation split leaves too few PSMs to early-stop against on a small # reference set, which is how a fit ends up worse than the model it started diff --git a/docs/source/usage.rst b/docs/source/usage.rst index 3f3a43a..2149d97 100644 --- a/docs/source/usage.rst +++ b/docs/source/usage.rst @@ -56,6 +56,19 @@ For calibration, pass a reference file with observed retention times using the ``--psm-file-reference`` option. If no reference file is provided, DeepLC attempts automatic calibration from high-confidence PSMs in the input file. +The reference is deduplicated before it is used: only the first PSM of each +peptidoform is kept, charge states included. A search result usually contains +the same peptidoform from many spectra with a different observed retention time +each time, and those repeats give the calibration conflicting targets while +weighing peptidoforms by how often they were identified. This is unconditional; +a caller that really wants every reference PSM to count fits a +:class:`~deeplc.calibration.Calibration` on its own targets and passes it to +:func:`~deeplc.core.predict_and_calibrate`, which uses an already fitted +calibration as given. When the repeats of a peptidoform disagree by a large +fraction of the observed retention-time range, DeepLC reports it: the reference +then mixes runs or contains low-confidence PSMs, and deduplication only hides +that. + For a full list of options: .. code-block:: sh diff --git a/pyproject.toml b/pyproject.toml index 1267661..b96964c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "deeplc" -version = "4.1.1" +version = "4.2.0" description = "DeepLC: Retention time prediction for (modified) peptides using Deep Learning." readme = "README.md" license = { file = "LICENSE" } diff --git a/tests/test_deduplication.py b/tests/test_deduplication.py new file mode 100644 index 0000000..2c1e2fe --- /dev/null +++ b/tests/test_deduplication.py @@ -0,0 +1,211 @@ +"""Reference deduplication: one PSM per peptidoform, the first observation.""" + +from __future__ import annotations + +import logging + +import numpy as np +from psm_utils import PSM, PSMList + +from deeplc import core +from deeplc._reference_selection import deduplicate_psms +from deeplc.calibration import SplineTransformerCalibration + +_PEPTIDES = [ + "AAGPSLSHTSGGTQSK", + "AGFAGDDAPR", + "AIQEYNQDK", + "AAYFGILEK", + "ADTQLDESSEQIDEEELTSK", + "AHQVVEDGYEFFAK", + "ALDQFVNFSEQK", + "AAPFSPAEK", + "VGAHAGEYGAEALER", + "LNLSPLGEEMR", +] + + +def _psms(pairs: list[tuple[str, float | None]], charge: int = 2) -> PSMList: + """Build a PSMList from (peptide, retention time) pairs.""" + return PSMList( + psm_list=[ + PSM(spectrum_id=str(i), peptidoform=f"{seq}/{charge}", retention_time=rt) + for i, (seq, rt) in enumerate(pairs) + ] + ) + + +def test_keeps_the_first_observation_of_each_peptidoform(): + """Repeats are dropped and the retention time of the first one survives.""" + psm_list = _psms([("PEPTIDEK", 10.0), ("PEPTIDEK", 40.0), ("ACDEFGHIK", 20.0)]) + + deduplicated = deduplicate_psms(psm_list) + + assert [str(p.peptidoform) for p in deduplicated] == ["PEPTIDEK/2", "ACDEFGHIK/2"] + assert [p.retention_time for p in deduplicated] == [10.0, 20.0] + + +def test_order_is_preserved(): + """The kept PSMs stay in the order they were given in.""" + psm_list = _psms([(s, float(i)) for i, s in enumerate(_PEPTIDES)]) + assert [str(p.peptidoform) for p in deduplicate_psms(psm_list)] == [ + str(p.peptidoform) for p in psm_list + ] + + +def test_idempotent(): + """Deduplicating an already deduplicated list changes nothing.""" + psm_list = _psms([("PEPTIDEK", 10.0), ("PEPTIDEK", 40.0), ("ACDEFGHIK", 20.0)]) + once = deduplicate_psms(psm_list) + assert len(deduplicate_psms(once)) == len(once) + + +def test_charge_states_are_duplicates_by_default(): + """Retention time does not depend on charge, so charge states are repeats.""" + psm_list = PSMList( + psm_list=[ + PSM(spectrum_id="1", peptidoform="PEPTIDEK/2", retention_time=10.0), + PSM(spectrum_id="2", peptidoform="PEPTIDEK/3", retention_time=11.0), + ] + ) + + assert len(deduplicate_psms(psm_list)) == 1 + assert len(deduplicate_psms(psm_list, ignore_charge=False)) == 2 + + +def test_modified_peptidoforms_are_not_duplicates(): + """A modification makes a different peptidoform, which elutes at a different time.""" + psm_list = PSMList( + psm_list=[ + PSM(spectrum_id="1", peptidoform="PEPTM[Oxidation]IDEK/2", retention_time=10.0), + PSM(spectrum_id="2", peptidoform="PEPTMIDEK/2", retention_time=12.0), + ] + ) + assert len(deduplicate_psms(psm_list)) == 2 + + +def test_exotic_peptidoforms_keep_their_identity(): + """ + The key is ``Peptidoform.modified_sequence``, so nothing but charge is ignored. + + Terminal, global and labile modifications all distinguish two peptidoforms, and a + modification label may itself contain a slash, which is why the charge is not cut off the + ProForma string by hand. + """ + distinct = [ + "PEPTIDEK/2", + "[Acetyl]-PEPTIDEK/2", + "PEPTIDEK-[Amidated]/2", + "PEPTM[Oxidation]IDEK/2", + "PEPT[Phospho]IDEK/2", + "<[Carbamidomethyl]@C>PEPCTIDEK/2", + "PROT[Phospho|+79.966]EIN/2", + ] + psm_list = PSMList( + psm_list=[ + PSM(spectrum_id=str(i), peptidoform=pf, retention_time=10.0 + i) + for i, pf in enumerate(distinct) + ] + ) + + assert len(deduplicate_psms(psm_list)) == len(distinct) + + +def test_charge_adducts_do_not_create_a_second_peptidoform(): + """``/2`` and ``/2[+2H]`` are the same peptidoform measured twice.""" + psm_list = PSMList( + psm_list=[ + PSM(spectrum_id="1", peptidoform="PEPTIDEK/2", retention_time=10.0), + PSM(spectrum_id="2", peptidoform="PEPTIDEK/2[+2H]", retention_time=11.0), + ] + ) + assert len(deduplicate_psms(psm_list)) == 1 + + +def test_warns_when_repeats_disagree_on_the_retention_time(caplog): + """ + A disagreement of the order of the gradient is a data problem, not jitter. + + Deduplication silently fixes the fit, so the case a user needs to hear about is when the + repeated peptidoform was not the same elution event at all. + """ + pairs = [(s, float(i)) for i, s in enumerate(_PEPTIDES)] + pairs.append((_PEPTIDES[0], 500.0)) # same peptidoform, 500 minutes later + + with caplog.at_level(logging.WARNING, logger="deeplc._reference_selection"): + deduplicate_psms(_psms(pairs)) + + assert any("disagreed on the observed retention time" in r.message for r in caplog.records) + + +def test_missing_retention_times_do_not_break_the_report(caplog): + """PSMs without an observed retention time are deduplicated like any other.""" + psm_list = _psms([("PEPTIDEK", None), ("PEPTIDEK", None), ("ACDEFGHIK", 20.0)]) + with caplog.at_level(logging.INFO, logger="deeplc._reference_selection"): + assert len(deduplicate_psms(psm_list)) == 2 + + +def _reference_with_duplicates() -> PSMList: + """Ten peptidoforms on a clean gradient, each repeated once at a wrong retention time.""" + pairs = [(s, 5.0 + 3.0 * i) for i, s in enumerate(_PEPTIDES)] + pairs += [(s, 100.0 - 2.0 * i) for i, s in enumerate(_PEPTIDES)] + return _psms(pairs) + + +def test_calibrate_always_uses_the_first_observations(): + """ + ``calibrate`` fits one point per peptidoform; there is no switch to turn that off. + + The reference holds each peptidoform twice: once on a clean 5 to 32 minute gradient and + once at a contradictory 82 to 100 minutes. The fit must reproduce the clean gradient. + """ + reference = _reference_with_duplicates() + targets = _psms([(s, None) for s in _PEPTIDES]) + + calibration = core.calibrate(reference, predict_kwargs={"device": "cpu"}) + predicted = core.predict(targets, return_matrix=True) + calibrated = calibration.transform(predicted[:, calibration.selected_model_head or 0]) + + assert np.isfinite(calibrated).all() + clean_low, clean_high = 5.0, 5.0 + 3.0 * (len(_PEPTIDES) - 1) + mean = float(calibrated.mean()) + assert clean_low <= mean <= clean_high, f"fit at {mean:.1f} left the clean gradient" + + +def test_a_prefitted_calibration_is_the_way_to_keep_the_repeats(): + """ + The escape hatch for the rare caller who wants every reference PSM to count. + + ``calibrate`` deduplicates unconditionally, so a caller who wants the repeats weighed fits + a ``Calibration`` on its own targets and passes it in; ``predict_and_calibrate`` then uses + it as given instead of fitting one. + """ + reference = _reference_with_duplicates() + psm_list = _psms([(s, None) for s in _PEPTIDES]) + + source = core.predict(reference, predict_kwargs={"device": "cpu"}, return_matrix=True) + own = SplineTransformerCalibration() + own.selected_model_head = 0 + own.fit(target=np.array(reference["retention_time"], dtype=np.float32), source=source[:, 0]) + + kept = core.predict_and_calibrate( + psm_list, psm_list_reference=reference, calibration=own, predict_kwargs={"device": "cpu"} + ) + deduplicated = core.predict_and_calibrate( + psm_list, psm_list_reference=reference, predict_kwargs={"device": "cpu"} + ) + + assert kept.shape == deduplicated.shape == (len(_PEPTIDES),) + assert not np.allclose(kept, deduplicated) + + +def test_predict_and_calibrate_runs_on_a_duplicated_reference(): + """One prediction per input PSM, in the input order, whatever the reference looks like.""" + psm_list = _psms([(s, None) for s in _PEPTIDES]) + result = core.predict_and_calibrate( + psm_list, + psm_list_reference=_reference_with_duplicates(), + predict_kwargs={"device": "cpu"}, + ) + assert result.shape == (len(_PEPTIDES),) + assert np.isfinite(result).all()