diff --git a/CHANGELOG.md b/CHANGELOG.md index 37a0f13..5539559 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,12 +32,29 @@ than a number. Both are recorded in every run's `manifest.json`. in the top 5,000; fragment charges 79% / 21%), `HCD2021` with charge-2 fragments only from precursor charge 3 14,412, the DIA-NN library 21,856. +- `peptides.tsv` and `proteins.tsv` (single-run and experiment-wide) carry + `is_transferred` and `transfer_q`, the acceptance basis of a match-between-runs row: + a transferred row keeps its grouped q next to the transfer q it was accepted at, a + tighter report threshold does not revoke a transfer that passed `mbr.q_transfer`, and a + protein group admitted through a transferred row carries the flag. The MBR worker's + augmented scored table gains `transfer_q` for it (docs/29 #19). +- `experiment_manifest.json` records the resolved `config_json` next to its hash, the + `model_identities` that produced the artifacts (RT source, fragment predictor, the + classifier that actually ran, feature schema, MBR strategy), the configured and + effective `quant.q_filter`, and input hashes taken at the start of the run rather than + at its end (docs/29 #15). - Dependabot covers the desktop application's Cargo dependencies (`/desktop`), and CI audits `desktop/Cargo.lock` with `cargo audit` next to the engine's lockfile, with one documented ignore (RUSTSEC-2024-0429: glib 0.18 through Tauri 2's gtk 0.18). ### Changed +- The candidate-audit rejection code `NO_PEAK_GROUP` is `DID_NOT_SURVIVE_EXTRACTION` + (`RejectionReason::DidNotSurviveExtraction`). The audit assigns it to every candidate + with no extracted row, and `extract` does not write the per-candidate table that would + separate presence, matched-fraction and gate failures, so the old name claimed a cause + the audit cannot see (docs/29 #16). The audit table has no versioned schema; the + metrics JSON gains `q_unit`. - MS2PIP 4.2.0 in every shipped environment (`env/docker-rescore.yml`, `env/console-ms2pip-requirements.txt`), and `env/mumdia-deeplc.yml` now carries `ms2pip==4.2.0` too, so one host environment serves DeepLC, MS2PIP and the `nn_torch` @@ -154,6 +171,30 @@ than a number. Both are recorded in every run's `manifest.json`. cancellation flag was written and never read. The waiter is now the only writer: it reads the intent after reaping the engine and publishes `cancelled`, `done` when the engine had already finished, or `failed`; until then the run shows "Stopping" (#14). +- Code review D, calibration, provenance, reporting (`docs/29`, findings 10, 15, 16, 19): + - LOESS retention-time calibration switched to the global least-squares line the + moment a query left the anchor range, while the grid just inside used the local fit, + and the two need not agree: on `y = 200 + 10x^2` (span 0.3) the prediction jumped + from 193.4 at `x = 1e-6` to 38.3 at `x = 0`, and from 1173.5 to 1018.4 at the top, + about 155 s discontinuities that misplaced gradient-edge peptides relative to their + extraction window. The map now continues the boundary local fit (its value and + slope) outside the range, and is continuous at both ends; the global line remains + only the degenerate fallback (#10). Measured on HYE B01 with the imported DIA-NN iRT + as the RT source and `native_tda`: 45,946 stripped peptides at `peptide_q_value` 1% + before, 45,957 after, at an unchanged 1.0% PSM-level decoy fraction and 6,410 protein + groups in both arms. 208,130 of 10.88 M candidates (1.9%) received a different window, + 194,698 of them with iRT above the anchor range, which the global line had placed + past the end of the 9,000 s run; the local fit places them at 8,578 to 9,100 s, and + extract accepted 454 more rows from them. A second pair on the DeepLC 4.1.1 + re-predicted precursor table (`w_rt` 414 s against 691 s): 48,533 stripped peptides in + both arms, PSM-q 1% targets 53,124 against 53,127 at the same decoy fraction, 6,519 + protein groups in both, 0.2% of candidates with a different window. Neutral on both RT + sources, which is what a boundary correction should be. + - The candidate audit's `passed_precursor_fdr` gate and `FAILED_PRECURSOR_FDR` reason + read the PSM `q_value`; they read `precursor_q`, the unit the label names, with the + PSM q as a recorded fallback on tables without it. A pooled scored table (several + `source` values) is refused, because the audit keys on `candidate_id` and would + attribute the last run's fate to every run (#16). - Desktop: the digest fields on the Search screen (missed cleavages, peptide length, charge range, carbamidomethyl, oxidation) now reach the engine on the built-in library path. They were read only by the DIA-NN library build, so with the built-in diff --git a/CLAUDE.md b/CLAUDE.md index 269396c..d01520a 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -137,8 +137,12 @@ Key semantics: eliminate a target against its decoy. Peptide-level q estimation subsequently performs picked target-decoy competition through the shared `base_peptide_id`; keep that pairing intact. -- `retain_top_peaks > 1` currently writes diagnostic peak alternatives only. - Those alternatives do not yet become feature/rescore rows. +- `extract.retain_top_peaks > 1` (default 1) writes the alternative peaks as + additional `psms_extracted` rows with `peak_rank >= 1` (plus a diagnostic + `.peaks.parquet`), `features` carries `peak_rank`, `compete` keys on it, and + `rescore` keeps one row per candidate and records `selected_peak_rank`. The + plumbing exists; what the default still lacks is entrapment validation on two + acquisitions. ## Validated sensitivity workflow @@ -203,9 +207,9 @@ The mechanism is peak-group formation rather than scoring: with most peaks gone, group. That reading comes from the cap dose-response, not from the audit ladder's own -label. `NO_PEAK_GROUP` cannot be used as evidence for it: `audit.rs` reads a +label. `DID_NOT_SURVIVE_EXTRACTION` cannot be used as evidence for it: `audit.rs` reads a per-candidate audit table that `extract` does not write (`emit_candidate_audit` -is unwired), so the reason map is always empty and the `_ => NoPeakGroup` +is unwired), so the reason map is always empty and the `_ => DidNotSurviveExtraction` catch-all absorbs presence failures, matched-fraction failures AND every extraction-gate rejection alike. Treat the label as "did not survive extract", and do not decompose it further until the audit table is actually produced. @@ -361,14 +365,15 @@ fine-tuning also is not guaranteed deterministic. single pooled run do not produce identical `q_value` columns. Batch to fit RAM, and compare per-run counts on `run_psm_q`. - Pooled rescore scales linearly, measured 0.834 ms/PSM on the streaming - backend. Two feature matrices, two widths: the Python worker's is - `n_psms x n_features x 4` bytes (f32), while the Rust `feats` that `rescore` - builds is `Vec>`, so `n_psms x n_features x 8` plus a heap allocation - and 24 bytes of spine per PSM. `native_tda` additionally runs all folds in + backend. Two feature matrices, one width: the Python worker's is + `n_psms x n_features x 4` bytes (f32), and the Rust `FeatureMatrix` that + `rescore` builds (`rescoring.rs`) is flat f32 as well, so the same + `n_psms x n_features x 4`. `native_tda` additionally runs all folds in parallel, each holding an owned standardised copy of its training slice, so its - peak is roughly `(1 + folds)x` the matrix. The stage logs the figure before it - allocates, and `rescore.max_feature_matrix_gib` turns exceeding a ceiling into - an error at startup rather than an OS kill hours in. + peak is roughly `(1 + folds)x` the matrix. `rescore.max_feature_matrix_gib` is + checked against that layout, from the parquet footers and the selected feature + count, before the allocation (docs/29 #11), so exceeding the ceiling is an error + at startup rather than an OS kill hours in. ### Rescore cost: handoff, feature selection, training-set reduction @@ -554,7 +559,8 @@ sensitivity result for it. Do not enable these by default from a single AIF count: -- model-visible top-K peaks (currently diagnostic sidecar only); +- model-visible top-K peaks (`extract.retain_top_peaks > 1`; implemented through + features, compete and rescore, default 1); - adaptive RT windows; - held-out RT window sizing (`rt_im_train.window_holdout_frac`). Implemented and measured on the AIF benchmark: +1.1% peptides with DeepLC 4.1.0 at unchanged @@ -581,8 +587,9 @@ and FDR population: the entrapment pool, HYE and AIF of docs/28 all ran under it The selected apex was historically correct/strongest only about 48-52% of the time while the correct peak appeared in the top five about 86-88%. Promoting top-K alternatives through features/rescore is therefore the best plausible -sensitivity project, but it needs a `candidate_id + peak_rank` contract and -entrapment validation before default activation. +sensitivity project. The `candidate_id + peak_rank` contract exists (`peak_rank` +on every extracted row, `selected_peak_rank` on the scored row, and the MBR worker +joins it); what default activation still needs is the entrapment validation. ## Coding conventions diff --git a/ci/check_smoke.py b/ci/check_smoke.py index 340ff4e..eb73301 100644 --- a/ci/check_smoke.py +++ b/ci/check_smoke.py @@ -425,10 +425,13 @@ def hashes(path): "run-experiment writes the experiment-wide peptides.tsv and proteins.tsv") if rep_pep.is_file(): hdr = rep_pep.read_text(encoding="utf-8").splitlines()[0].split("\t") + # The two trailing columns are the acceptance basis of a match-between-runs + # row (docs/29 #19); without MBR they read `false` and empty. c.ok(hdr == ["precursor", "stripped_sequence", "charge", "protein", "q_value", - "score", "n_runs", "quantity_a", "quantity_b"], - "experiment peptides.tsv has the experiment-wide columns and one quantity " - "column per run", "\t".join(hdr)) + "score", "n_runs", "quantity_a", "quantity_b", "is_transferred", + "transfer_q"], + "experiment peptides.tsv has the experiment-wide columns, one quantity " + "column per run and the transfer basis", "\t".join(hdr)) # Rows are asserted on the standalone rewrite at q 0.05 (smoke.sh): the pooled # peptide-level q of this fixture cannot reach 1 percent, so the root pair is # legitimately header-only at the default threshold. @@ -444,8 +447,10 @@ def hashes(path): "identical inputs: every precursor is quantified in both runs") if rep_prot.is_file(): hdr = rep_prot.read_text(encoding="utf-8").splitlines()[0].split("\t") - c.ok(hdr == ["protein_group", "q_value", "n_runs", "lfq_a", "lfq_b"], - "experiment proteins.tsv has one LFQ column per run", "\t".join(hdr)) + c.ok(hdr == ["protein_group", "q_value", "n_runs", "lfq_a", "lfq_b", + "is_transferred", "transfer_q"], + "experiment proteins.tsv has one LFQ column per run and the transfer " + "basis", "\t".join(hdr)) for r in ("a", "b"): c.ok(not (exp / r / "peptides.tsv").exists(), f"no per-run peptides.tsv under {r}: the grouped q is experiment-wide") diff --git a/desktop/src-tauri/src/run.rs b/desktop/src-tauri/src/run.rs index e5735d3..b785d5d 100644 --- a/desktop/src-tauri/src/run.rs +++ b/desktop/src-tauri/src/run.rs @@ -185,7 +185,8 @@ pub struct Snapshot { /// /// The interface needs this to label the result counts: an experiment-wide /// rescore groups the q columns experiment-wide, so those counts are NOT per - /// file, and `run-experiment` writes no report at all. + /// file, and the `peptides.tsv` / `proteins.tsv` at the experiment root are the + /// experiment-wide report, not a per-run one. pub experiment: bool, /// Stop was requested. The status stays `running` until the engine has actually /// been reaped, and the interface shows "Stopping" meanwhile. @@ -1361,10 +1362,10 @@ mod multifile_tests { #[test] fn a_pooled_experiment_reports_its_combined_table_and_says_so() { - // `run-experiment` writes `scored_combined.parquet` and never calls the report - // stage, so reading only `psms_scored.parquet.report.json` left the results - // screen blank after every experiment. And the counts it does yield are - // experiment-wide: the grouped q columns are grouped across the whole + // `run-experiment` writes `scored_combined.parquet` (and an experiment-wide + // TSV pair at the root), so reading only `psms_scored.parquet.report.json` left + // the results screen blank after every experiment. And the counts it does + // yield are experiment-wide: the grouped q columns are grouped across the whole // experiment, so a per-file reading of them is diluted by ~1/n_runs. let dir = std::env::temp_dir().join("mumdia-results-experiment"); let _ = std::fs::remove_dir_all(&dir); @@ -1383,7 +1384,7 @@ mod multifile_tests { assert!(r.experiment_wide, "a combined table is experiment-wide"); assert_eq!(r.peptides_1pct, 7); assert_eq!(r.precursors_1pct, 8); - // And it writes no report, so neither TSV exists. + // This fixture wrote no TSV, so neither is reported present. assert!(!r.has_peptides_tsv && !r.has_proteins_tsv); // A single run's own report wins, and is not labelled experiment-wide. diff --git a/docs/02_config_and_data_model.md b/docs/02_config_and_data_model.md index 966423d..a2e1406 100644 --- a/docs/02_config_and_data_model.md +++ b/docs/02_config_and_data_model.md @@ -385,7 +385,7 @@ the smaller `stage_order` and `Reported` never overrides a real rejection; | 6 | `RtPruned` / `RT_PRUNED` | candidate generation (B) | | 7 | `CandidateCapReached` / `CANDIDATE_CAP_REACHED` | candidate generation (B) | | 8 | `NoFragmentTraces` / `NO_FRAGMENT_TRACES` | extraction (C, D) | -| 9 | `NoPeakGroup` / `NO_PEAK_GROUP` | extraction (C, D) | +| 9 | `DidNotSurviveExtraction` / `DID_NOT_SURVIVE_EXTRACTION` (was `NoPeakGroup` / `NO_PEAK_GROUP` before docs/29 #16) | extraction (C, D) | | 10 | `PeakNotSelected` / `PEAK_NOT_SELECTED` | peak/peptide ranking (E) | | 11 | `OutcompetedByTarget` / `OUTCOMPETED_BY_TARGET` | competition (G) | | 12 | `OutcompetedByDecoy` / `OUTCOMPETED_BY_DECOY` | competition (G) | diff --git a/docs/04_convert.md b/docs/04_convert.md index e68009f..9617003 100644 --- a/docs/04_convert.md +++ b/docs/04_convert.md @@ -321,10 +321,10 @@ and 80.3% of what a library-free DIA-NN 2.2.0 search reports on the same file. The mechanism is downstream: with most peaks gone, candidates cannot assemble enough distinct matched fragments to satisfy `extract.presence_min_fragments` (`config.rs:523`, default 3 at `config.rs:690`), so they fail peak-group -formation and are recorded with rejection code `NO_PEAK_GROUP` +formation and are recorded with rejection code `DID_NOT_SURVIVE_EXTRACTION` (`rejection.rs:62`). This was confirmed with `mumdia audit` on the capped arm, restricted to the peptides that same DIA-NN search reports as present (78,782 -distinct `Stripped.Sequence` at DIA-NN `Q.Value` <= 0.01): 49,105 of 78,782 (62.3%) stopped at `candidate_generated` with `NO_PEAK_GROUP`, against only +distinct `Stripped.Sequence` at DIA-NN `Q.Value` <= 0.01): 49,105 of 78,782 (62.3%) stopped at `candidate_generated` with `DID_NOT_SURVIVE_EXTRACTION`, against only 5,380 lost to FDR and 355 to competition, and a counterfactual replay on the uncapped artifact recovered 41,948 (85.4%) of them. The loss is therefore extraction-side, not a scoring or competition effect. See docs/09_extract.md for diff --git a/docs/08_rt_im_train.md b/docs/08_rt_im_train.md index 28141cf..4fa1227 100644 --- a/docs/08_rt_im_train.md +++ b/docs/08_rt_im_train.md @@ -179,9 +179,10 @@ rt_im_train.rs:136-140). `predict` is a closure (rt_im_train.rs:142-150): it returns `NaN` when calibration is unavailable, otherwise `loess.predict(irt)` when a LOESS model exists, else -`slope * irt + intercept`. The linear coefficients are therefore both the primary -map (Linear method) and the extrapolation fallback (LOESS outside its training -range; see below). +`slope * irt + intercept`. The linear coefficients are the primary map under the +Linear method and, under LOESS, only the degenerate fallback (fewer than four anchors, +or a local window without spread). Outside its training range LOESS continues the +local fit at the nearest boundary, not the global line (see below). The math: @@ -191,7 +192,7 @@ The math: returns `(0, mean(ys))` (a constant map, calibrate.rs:8-16); a near-zero denominator (all `x` equal) returns `(0, Sy/n)` (calibrate.rs:22-24). - **`Loess::fit`** (calibrate.rs:42-83) first computes the global linear fit as - the extrapolation fallback (calibrate.rs:43), sorts the points by `x` + the degenerate fallback (calibrate.rs:43), sorts the points by `x` (calibrate.rs:45-48), and if fewer than 4 points are present just fills the grid from the linear line (calibrate.rs:50-66). Otherwise the local window size is `k = clamp(ceil(span * n), 3, n)` (calibrate.rs:67) and it evaluates a @@ -513,9 +514,9 @@ cheaper option with equal RT residuals on the one run measured here. | `candidate_window` | rt_im_train.rs:65-70 | Builds `(cal, lo, hi)`; returns `(NaN, -inf, +inf)` when calibrated RT or width is absent. | | `rt_im_train::run` | rt_im_train.rs:72-354 | The stage: join iRT, select anchors, fit, window, apply, write. | | `linear_fit` | calibrate.rs:6-28 | OLS `y = slope*x + intercept` with degenerate-case guards. | -| `Loess` | calibrate.rs:31-37 | Grid-based local-linear smoother; carries a linear fallback for extrapolation. | +| `Loess` | calibrate.rs:31-37 | Grid-based local-linear smoother; carries the boundary local slopes for extrapolation and the global line as the degenerate fallback. | | `Loess::fit` | calibrate.rs:42-83 | Sorts anchors, builds a `grid_n`-point local-linear grid, `k = clamp(ceil(span*n),3,n)`. | -| `Loess::predict` | calibrate.rs:87-102 | Grid interpolation inside range, linear extrapolation outside. | +| `Loess::predict` | calibrate.rs:87-102 | Grid interpolation inside range; outside it, the boundary grid value continued with the boundary local slope, so the map is continuous at both ends. Until docs/29 #10 it switched to the global line there, which on `y = 200 + 10x^2` (span 0.3) jumped from 193.4 to 38.3 at `x = 0` and from 1173.5 to 1018.4 at the top: about 155 s discontinuities for gradient-edge peptides. Measured on HYE B01 with the imported iRT and `native_tda`: 45,946 stripped peptides at 1% before, 45,957 after, decoy fraction unchanged; 1.9% of candidates got a different window, almost all with iRT above the anchor range, which the global line had placed past the end of the run. With the DeepLC 4.1.1 re-predicted precursors (`w_rt` 414 s): 48,533 in both arms, 0.2% of windows moved. | | `local_linear` | calibrate.rs:107-153 | Tricubic-weighted local least squares at one point. | | `percentile` | calibrate.rs:156-164 | Nearest-rank percentile: sorts a copy, `rank = round(p.clamp(0,1)*(len-1))`. Not interpolated. Empty input returns 0.0. | | `CalibrationMethod` | config.rs:56-61 | Enum `{ Loess, Linear, None }`; default `Loess`. `None` is rejected at load. | @@ -618,9 +619,11 @@ though the enum variant still exists. and therefore the whole run, non-reproducible. Treat it as an accuracy lever, not a deterministic default. - **`slope`/`intercept` are emitted in `cal.json` whenever calibration is available, - including under LOESS** (rt_im_train.rs:128, 286-287, 313-314). They are the LOESS - extrapolation fallback, not dead values; do not assume they were unused when - `method == "loess"`. They are serialized as `null` only when calibration is + including under LOESS** (rt_im_train.rs:128, 286-287, 313-314). Under LOESS they + are the degenerate fallback (fewer than four anchors, or a local window without + spread), not the extrapolation model: since docs/29 #10 the map continues the + boundary local fit outside the anchor range. Do not read them as the calibration + when `method == "loess"`. They are serialized as `null` only when calibration is unavailable (`n_train < 2`), since the fit is not computed in that case. - **`CalibrationMethod::None` still exists but is rejected** at config load (config.rs:1336-1342). The stage would otherwise fall through to the linear path diff --git a/docs/09_extract.md b/docs/09_extract.md index 5efb1e2..1cc7e5d 100644 --- a/docs/09_extract.md +++ b/docs/09_extract.md @@ -327,21 +327,22 @@ rejected before any gate or score is reached. Measured on a 50-window Orbitrap DIA run, `--top-peaks-ms2 300` discarded 78.6% of MS2 peaks; a `mumdia audit` ladder restricted to peptides that an external library-free search confirms are present in the run showed 49,105 of 78,782 (62.3%) stopping at -`candidate_generated` with rejection code `NO_PEAK_GROUP` (`rejection.rs:62`). +`candidate_generated` with rejection code `DID_NOT_SURVIVE_EXTRACTION` (`rejection.rs:62`). Only 5,380 were lost to FDR and 355 to competition. Replaying those candidates against the uncapped artifact recovered 41,948 of the 49,105 (85.4%), which is what identifies the cap as the cause. -`NO_PEAK_GROUP` does not mean "never assembled `presence_min_fragments` distinct -fragments", although this document previously said so. `audit.rs` reads a +`DID_NOT_SURVIVE_EXTRACTION` (named `NO_PEAK_GROUP` until docs/29 #16) does not mean +"never assembled `presence_min_fragments` distinct fragments", although this document +previously said so. `audit.rs` reads a per-candidate audit table that `extract` does not write -- `emit_candidate_audit` -is unwired -- so the reason map is always empty and the `_ => NoPeakGroup` +is unwired -- so the reason map is always empty and the `_ => DidNotSurviveExtraction` catch-all absorbs presence failures, matched-fraction failures and every extraction-gate rejection alike. Read it as "did not survive extract". The attribution to the peak cap above stands on the replay experiment, not on the label. -When `NO_PEAK_GROUP` dominates an audit ladder, check the conversion cap before +When `DID_NOT_SURVIVE_EXTRACTION` dominates an audit ladder, check the conversion cap before tuning the presence thresholds or gates here; see docs/04_convert.md for the peak census and the cap dose-response. diff --git a/docs/12_quant_lfq_align_mbr_report_audit.md b/docs/12_quant_lfq_align_mbr_report_audit.md index b7cfbb8..9a2a05a 100644 --- a/docs/12_quant_lfq_align_mbr_report_audit.md +++ b/docs/12_quant_lfq_align_mbr_report_audit.md @@ -112,10 +112,11 @@ per accepted transfer (`mbr_worker.py:254`): `candidate_id`, `source`, `peptidof `transfer_q` (10 columns). When there are no transfer candidates at all the worker short-circuits and writes a placeholder table with a single empty `candidate_id` column (`pa_write_empty`, `mbr_worker.py:289`), so `.parquet` always exists. -Optionally writes an augmented scored table (`--out-psms-scored`) that lowers each accepted -transfer's `q_value` to `min(q_value, transfer_q)` on the matching `(candidate_id, -source)` row and adds an `is_transferred` bool (`mbr_worker.py:272`); this requires the -scored table to carry a `source` column. All MBR outputs are written by Python and +Optionally writes an augmented scored table (`--out-scored`) that lowers each accepted +transfer's PSM q columns to `min(q, transfer_q)` on the matching `(candidate_id, +source)` row and adds an `is_transferred` bool plus a `transfer_q` column (the accepted +q on transferred rows, NaN elsewhere; `mbr_worker.py`); this requires the scored table +to carry a `source` column. All MBR outputs are written by Python and have no `report.json`. The unreachable re-extraction tier (`--emit-transfer-targets `) instead writes @@ -130,8 +131,14 @@ Consumes `psms_scored.parquet` and optionally the two quant tables. The CLI take (`main.rs:796-797`); `run` passes the run out-dir and `q_threshold = quant.q_threshold` (`run.rs:484-491`). `peptides.tsv` header (`report.rs:95`): `precursor`, `stripped_sequence`, `charge`, `protein`, -`q_value`, `score`, `quantity`. `proteins.tsv` header (`report.rs:134`): -`protein_group`, `q_value`, `quantity`. No Parquet or `report.json` is written; +`q_value`, `score`, `quantity`, `is_transferred`, `transfer_q`. `proteins.tsv` header +(`report.rs:134`): `protein_group`, `q_value`, `quantity`, `is_transferred`, +`transfer_q`. The last two export the acceptance basis: a row is written when it is a +target and either its grouped q passes the threshold or it is an accepted +match-between-runs transfer, and a transferred row keeps its grouped q (usually 1.0) +next to the transfer q it was accepted at. A tighter threshold does not revoke a +transfer; a transfer is not protein-group confidence (docs/29 #19). No Parquet or +`report.json` is written; `report::run` returns `(n_precursors, n_protein_groups)`, and the `mumdia report` handler prints a one-line summary (`main.rs:806`). @@ -423,7 +430,8 @@ The worker (`scripts/mbr_worker.py`) implements two tiers: The M5 augmented scored output (`--out-scored`, `mbr_worker.py:272`) requires the scored table to have a `source` column and matches transfers on `(candidate_id, -source)`, taking `min(q_value, transfer_q)` and setting `is_transferred`. The worker +source)`, taking `min(q, transfer_q)` on the PSM q columns and setting `is_transferred` +and `transfer_q`. The worker prints a validation summary: accepted counts per run, how many permuted-RT null draws fall inside the accepted RT window, and the window `delta_star` at `q_transfer`. It no longer prints a "decoy fraction" among accepted transfers: every transfer candidate @@ -463,16 +471,20 @@ presentation value rounded to one decimal. The search space is every library precursor (`candidate_id`, `peptidoform`, `charge`, `label`, `protein`). Survivor sets are built as `HashSet` of `candidate_id` from -`psms` (extracted), `competed`, and `scored`; scored also yields `q_by_cid` and -optional `pepq_by_cid` (peptide-level q, only present in some scored schemas, -`audit.rs:90`). `load_extract_reasons` (`audit.rs:51`) tries to read a +`psms` (extracted), `competed`, and `scored`; scored also yields `q_by_cid`, read from +`precursor_q` (the unit the `passed_precursor_fdr` label names; the PSM `q_value` only +on an older table without that column, recorded as `q_unit` in the metrics JSON; +docs/29 #16) and optional `pepq_by_cid` (peptide-level q, only present in some scored +schemas). A scored table with more than one `source` is refused: the audit keys on +`candidate_id`, so a pooled table would attribute one run's fate to every run; audit +the per-run split tables. `load_extract_reasons` (`audit.rs:51`) tries to read a `.audit.parquet` sidecar to refine the extract-stage bucket, but nothing in the current chain writes that file (see the gotchas), so the map is empty and every -extract-stage loss buckets to `NO_PEAK_GROUP`. For each candidate the earliest +extract-stage loss buckets to `DID_NOT_SURVIVE_EXTRACTION`. For each candidate the earliest rejection reason is assigned along the ladder (`audit.rs:136`): - not in `extracted` -> refined from the sidecar (`NO_FRAGMENT_TRACES`/`NO_VALID_FRAGMENTS`/`PEAK_NOT_SELECTED`/`RT_PRUNED`/ - `WRONG_ISOLATION_WINDOW`, `audit.rs:139`) or the generic `NO_PEAK_GROUP` when no + `WRONG_ISOLATION_WINDOW`, `audit.rs:139`) or the generic `DID_NOT_SURVIVE_EXTRACTION` when no sidecar (the only outcome today); - extracted but not in `competed` -> `OUTCOMPETED_BY_DECOY` (decoy) or `OUTCOMPETED_BY_TARGET` (target); @@ -641,14 +653,14 @@ but do not affect the wired `mumdia mbr` path. requires the peptide gate. A candidate can therefore have `reported=true` yet `rejection_reason=FAILED_PEPTIDE_FDR`. Treat `rejection_reason` as authoritative. - **audit reason coverage.** In the current chain `audit.rs` can only ever emit - `NO_PEAK_GROUP`, `OUTCOMPETED_BY_TARGET`/`OUTCOMPETED_BY_DECOY`, `FAILED_PRECURSOR_FDR`, + `DID_NOT_SURVIVE_EXTRACTION`, `OUTCOMPETED_BY_TARGET`/`OUTCOMPETED_BY_DECOY`, `FAILED_PRECURSOR_FDR`, `FAILED_PEPTIDE_FDR`, and `REPORTED`. The five refined extract codes (`NO_FRAGMENT_TRACES`, `NO_VALID_FRAGMENTS`, `PEAK_NOT_SELECTED`, `RT_PRUNED`, `WRONG_ISOLATION_WINDOW`) are matched by `load_extract_reasons` (`audit.rs:139`) but need the absent sidecar to fire. The six remaining enum variants (`PEPTIDE_NOT_GENERATED`, `MODIFICATION_NOT_ALLOWED`, `CHARGE_OUT_OF_RANGE`, `PRECURSOR_MZ_OUT_OF_RANGE`, `CANDIDATE_CAP_REACHED`, `REMOVED_DURING_REPORTING`) have - no producer at all: even a sidecar string of that name falls to the `_ => NoPeakGroup` + no producer at all: even a sidecar string of that name falls to the `_ => DidNotSurviveExtraction` arm (`audit.rs:145`). - **report.json coverage is partial.** `quant` writes reports for `peptide_quant`, `protein_quant`, and emitted `fragment_quant`; the peak-bounds diff --git a/docs/13_sidecars.md b/docs/13_sidecars.md index e0bf32b..3b390ae 100644 --- a/docs/13_sidecars.md +++ b/docs/13_sidecars.md @@ -155,8 +155,9 @@ code. - OUT `.parquet`: one row per accepted transfer with `candidate_id, source, peptidoform, charge, protein_group, label, expected_rt, observed_rt, rt_delta, transfer_q` (`mbr_worker.py:254-265`). Optional - `--out-scored` writes the scored table with accepted transfers' `q_value` - lowered to `transfer_q` and an `is_transferred` flag added. Optional + `--out-scored` writes the scored table with accepted transfers' PSM q columns + lowered to `transfer_q`, an `is_transferred` flag and a `transfer_q` column (NaN + on non-transferred rows) added. Optional `--emit-transfer-targets` writes per-run `run_windows`-format tables (`candidate_id, rt_pred_cal, rt_lo, rt_hi, im_*`) plus a permuted-RT decoy file for the re-extraction tier (`mbr_worker.py:142-151`). diff --git a/docs/15_data_dictionary.md b/docs/15_data_dictionary.md index 14279fb..a260ef5 100644 --- a/docs/15_data_dictionary.md +++ b/docs/15_data_dictionary.md @@ -817,11 +817,19 @@ Each `ArtifactRecord`: Not Parquet, but these are the files most users actually read, so they are documented here rather than left to the source. Written by `report` -(`stages/report.rs`), which has two call sites: inside `run` (`run.rs:484`) and -the `mumdia report` handler (`main.rs:798`). `run-experiment` never calls it, so -an experiment output tree contains neither TSV at any level; take per-run counts -from the per-run split `scored.parquet` tables on `run_psm_q`, or invoke -`mumdia report` yourself. +(`stages/report.rs`): `report::run` inside `run` and the `mumdia report` handler, +and `report::run_experiment` inside `run-experiment` (also `mumdia report +--experiment-dir`), which writes one experiment-wide pair at the experiment root +with an `n_runs` column and one quantity column per run and no per-run TSVs; take +per-run counts from the per-run split `scored.parquet` tables on `run_psm_q`. + +**Acceptance rule, both writers.** A row is written when it is a target and either +its grouped q is at or under the threshold or it is flagged `is_transferred` by +`mumdia mbr`. The two columns `is_transferred` and `transfer_q` export that basis +(docs/29 #19): a transferred row keeps its grouped q, which is usually 1.0, and a +tighter `--q` does not revoke a transfer that already passed `mbr.q_transfer`. A +transfer is a peptide-level acceptance, not protein-group confidence; a protein group +admitted through a transferred row carries the flag so the reader can tell. Values are formatted for reading, not for analysis: `q_value` is printed to six decimals, `score` to four, and `quantity` to one @@ -843,6 +851,8 @@ the best q (`report.rs:90,105-108`). Targets only, filtered on | `q_value` | `peptide_q_value` | **base-peptide** q, six decimals | | `score` | `score` | rescorer score, four decimals | | `quantity` | `peptide_quant.parquet` joined on `(peptidoform, charge)` | empty when the precursor was not quantifiable or no `--peptide-quant` was passed (`report.rs:39-45,61-72`) | +| `is_transferred` | `is_transferred` (`mumdia mbr`), `false` without it | `true` when the row was admitted as a match-between-runs transfer rather than on its q | +| `transfer_q` | `transfer_q` (`mumdia mbr`) | the transfer q the row was accepted at, six decimals; empty on every other row | The row unit and the filter column deliberately disagree, and this is the single most common misreading of MuMDIA output: rows are precursors @@ -861,6 +871,8 @@ with a non-empty group, filtered on `pg_q_value <= --q` (`report.rs:137`). | `protein_group` | `protein_group` | accession set as grouped by rescore | | `q_value` | `pg_q_value` | protein-group q, six decimals | | `quantity` | `protein_group_quant.parquet` joined on `protein_group` | empty when not quantifiable or no `--protein-quant` was passed | +| `is_transferred` | `is_transferred` of the admitting row | `true` when the group entered through a transferred peptide row; a transfer is not protein-group confidence | +| `transfer_q` | `transfer_q` of the admitting row | empty unless `is_transferred` | Both grouped q columns (`peptide_q_value`, `pg_q_value`) are written only to each group's single winning row, so under an experiment-wide rescore the grouping is diff --git a/docs/17_troubleshooting.md b/docs/17_troubleshooting.md index f399d22..f12a41b 100644 --- a/docs/17_troubleshooting.md +++ b/docs/17_troubleshooting.md @@ -27,7 +27,7 @@ is the source of truth if it has moved. | Two runs of the same command give slightly different identification counts | A nondeterministic opt-in path is on: DeepLC fine-tune (no seed) or the PyTorch NN rescorer (approximately reproducible) | Accept the trade for the ID gain, or leave those paths off; set `MUMDIA_NN_SEEDS > 1` to average out NN variance | | Peptidoforms silently get iRT 0.0 | DeepLC returned no iRT for some peptidoforms; they are anchored at 0.0 with a warning | Check the DeepLC env and the `n_irt_missing` warning; unmatched peptidoforms are a known foot-gun | | Peptide count is a fraction of expectation while the empirical decoy fraction still sits at the target | `--top-peaks-ms2` truncated the MS2 peaks at convert time and baked the truncation into the spectra artifact (`convert.rs:76-79`) | Reconvert uncapped or with a much larger cap; the right cap is acquisition-specific, not a universal preset | -| `mumdia audit` attributes most missed peptides to `NO_PEAK_GROUP` at `candidate_generated` | The same peak truncation leaves too few surviving fragments to satisfy `extract.presence_min_fragments` (`extract.rs:1926-1931`) | Raise or remove `--top-peaks-ms2` and reconvert; do not lower `presence_min_fragments` to compensate | +| `mumdia audit` attributes most missed peptides to `DID_NOT_SURVIVE_EXTRACTION` at `candidate_generated` | The same peak truncation leaves too few surviving fragments to satisfy `extract.presence_min_fragments` (`extract.rs:1926-1931`) | Raise or remove `--top-peaks-ms2` and reconvert; do not lower `presence_min_fragments` to compensate | | Every modified form of a peptide is missing from the output; `precursor_q` reports exactly 1.000 precursors per peptide | `compete.group_by = base_peptide` keys the stripped sequence, not the precursor (`compete.rs:88`), and deletes all but the top-scoring sibling before rescore (`compete.rs:319-340`) | Set `compete.group_by = peptidoform_charge` (`compete.rs:93-98`); required for any PTM search | | `cal.json` reports a small `rt_residual_abs_median_s` but RT windows still miss peptides | The residual is measured on the same anchors the calibration was fitted to (`rt_im_train.rs:137`, `rt_im_train.rs:177-185`), so it is in-sample | Treat it as a fit diagnostic; size external RT tolerances from an out-of-sample comparison | | RT windows in a PTM search behave as if the modification were absent | The imported library gave every modform of a stripped peptide the same `predicted_irt` | Check per-stripped-peptide variance of `predicted_irt` in the library; re-predict iRT per peptidoform | @@ -222,7 +222,7 @@ The mechanism is peak-group formation, not scoring. With most peaks gone, `extract.presence_min_fragments` (default 3) cannot be met and the candidate returns no peak group (`extract.rs:1926-1931`). `mumdia audit` on the capped arm, restricted to peptides the reference search confirms are present, put -49,105 of 78,782 (62.3%) at `candidate_generated` with `NO_PEAK_GROUP`; only +49,105 of 78,782 (62.3%) at `candidate_generated` with `DID_NOT_SURVIVE_EXTRACTION`; only 5,380 were lost to FDR and 355 to competition. A counterfactual replay against the uncapped artifact recovered 41,948 of those 49,105 (85.4%). Do not lower `presence_min_fragments` to compensate: that trades a real fragment requirement diff --git a/docs/18_findings_and_decisions.md b/docs/18_findings_and_decisions.md index 741f1b1..c5e149c 100644 --- a/docs/18_findings_and_decisions.md +++ b/docs/18_findings_and_decisions.md @@ -150,7 +150,7 @@ The mechanism is peak-group formation, not scoring. With most peaks removed, (`rust/mumdia/crates/mumdia-core/src/config.rs:523`, default at `:690`) cannot be satisfied, so real peptides never form a peak group. `mumdia audit` on the capped arm, restricted to peptides DIA-NN confirms are present, shows 49,105 of 78,782 -(62.3 percent) stopped at `candidate_generated` with `NO_PEAK_GROUP`, against +(62.3 percent) stopped at `candidate_generated` with `DID_NOT_SURVIVE_EXTRACTION`, against only 5,380 lost to FDR and 355 lost to competition. A counterfactual replay on the uncapped artifact recovered 41,948 of those 49,105 (85.4 percent). diff --git a/docs/20_sensitivity_and_quantification_playbook.md b/docs/20_sensitivity_and_quantification_playbook.md index 8af4b22..ceba92b 100644 --- a/docs/20_sensitivity_and_quantification_playbook.md +++ b/docs/20_sensitivity_and_quantification_playbook.md @@ -131,7 +131,7 @@ real peptides never form a peak group. On a 50-window Orbitrap DIA run, a cap of 300 discarded 78.6% of all MS2 peaks and truncated 85.5% of spectra, and `mumdia audit` restricted to peptides DIA-NN 2.2.0 confirms are present showed 49,105 of 78,782 (62.3%) stopped at `candidate_generated` with -`NO_PEAK_GROUP`, versus 5,380 lost to FDR and 355 lost to competition. Peptides.tsv rows at `peptide_q_value` <= 0.01 +`DID_NOT_SURVIVE_EXTRACTION`, versus 5,380 lost to FDR and 355 lost to competition. Peptides.tsv rows at `peptide_q_value` <= 0.01 fell from 63,237 uncapped to 25,425 at cap 300 with the decoy fraction at 0.99% in both arms, so this is lost sensitivity and not a changed threshold. The canonical treatment of the flag, including the full peak census and the cap @@ -295,7 +295,7 @@ remaining lever, then rescore, then seed, all on faint signal. Before reading a presence/apex loss on any other run as an extraction problem, check the conversion cap. A cap that truncates most spectra produces exactly the -same signature (`NO_PEAK_GROUP` at `candidate_generated`) while the cause is the +same signature (`DID_NOT_SURVIVE_EXTRACTION` at `candidate_generated`) while the cause is the converted artifact, not the extraction thresholds. On the 50-window run above, 62.3% of the confirmed-present peptides were lost this way at a cap of 300. diff --git a/docs/29_code_review_2026-09-07.md b/docs/29_code_review_2026-09-07.md index 00462cd..a7263a0 100644 --- a/docs/29_code_review_2026-09-07.md +++ b/docs/29_code_review_2026-09-07.md @@ -257,3 +257,42 @@ The probes were written under `C:/Users/robbi/AppData/Local/Temp/mumdia_review_2 No full real-data DIA experiment, model retraining benchmark, desktop GUI end-to-end run, fresh release/container build, or current external vulnerability scan was completed for this review. Optional-sidecar-enabled test completion was not established, so it is not counted as validation. Static desktop findings remain source-confirmed rather than reproduced in a running GUI. The findings support the concrete corrections above; they do not quantify their prevalence or sensitivity effect on biological datasets. **Agreed delivery:** A–D plus the separate local hygiene action, with P1 data-integrity defects first. Add the regression fixtures where each finding is fixed. D's LOESS correction requires the HYE B01 before/after count before merge. Deferred research, architecture, resume, and export projects are not prerequisites for these four PRs. Future promotion of sensitivity or MBR defaults still requires the repository's scientific validation; this maintenance plan does not claim that validation has already happened. + +--- + +## Implementation status (2026-09-07) + +All four work packages are implemented as stacked pull requests, each on the previous one and +the first on #61 (`perf/ms2pip-worker-throughput`); merge in order. Each PR carries its tests +and CHANGELOG entries; CI was green on every one at the time of writing. + +| Package | PR | Findings | Notes | +|---|---|---|---| +| A. Engine data integrity | #62 `review/a-data-integrity` | 1, 2, 4, 9, 11, 17, 18, 21 | Predictor misses drop the candidate and its pair; numeric domains validated at load; scanners on `git ls-files`. | +| B. Workers | #63 `review/b-workers` | 3, 6, 7, 8, 12, 20 | Single-class fold is an error; MBR q uses the `+1` pseudocount and prints the permuted-null draws; `selected_peak_rank` joined. | +| C. Desktop and output ownership | #64 `review/c-desktop` | 5, 13, 14 | Frontend start guard plus backend results-folder reservation; converter probe with the request's configuration and the engine's msconvert fallback rule; single terminal-state writer; Dependabot and `cargo audit` for `desktop/Cargo.lock` (one documented ignore). | +| D. Provenance, reporting, calibration, docs | #65 `review/d-provenance` | 10, 15, 16, 19 | Boundary-continuous LOESS; experiment manifest with `config_json`, model identities, effective `q_filter`, inputs hashed at start; audit on `precursor_q`, single source, `NO_PEAK_GROUP` renamed `DID_NOT_SURVIVE_EXTRACTION`; `is_transferred` / `transfer_q` in both TSVs. | +| Local hygiene | (no PR) | 21 | 173 scratch copies moved to `C:/Users/robbi/mumdia_bench/untracked_backup_2026-09-06/covr2/`, nothing deleted. | + +**Finding 10 merge check.** Merge check, HYE B01 (`LFQ_Orbitrap_AIF_Condition_B_Sample_Alpha_01`, DIA-NN library with +its imported iRT as the RT source, `native_tda`, 32 threads, one binary changed): stripped +peptides at `peptide_q_value` 1% 45,946 before and 45,957 after; PSM-q 1% targets 50,332 +against 50,335 at an unchanged 1.0% decoy fraction; protein groups 6,410 both. 208,130 of +10.88 M candidates (1.9%) received a different RT window, 194,698 of them with iRT above the +anchor range (170), which the old global line had placed at 10,625 to 11,946 s in a 9,000 s +run; the continued local fit places them at 8,578 to 9,100 s. Extract accepted 454 more rows +and the run took 37:30 against 34:14. +Second pair, same file and settings with the DeepLC 4.1.1 base-model re-predicted precursor +table as the library (`rt_im_train.library_irt = auto` output of 2026-09-06, passed as +`--lib-precursors`; `w_rt` 414 s against 691 s, in-sample residual median 78 s against 130 s): +stripped peptides at 1% 48,533 in both arms; PSM-q 1% targets 53,124 against 53,127 at an +unchanged 1.0% decoy fraction; protein groups 6,519 in both; 22,779 of 10.88 M candidates +(0.2%) received a different window (median move 424 s); extract accepted 27 more rows; wall +26:06 against 25:46. The correction is neutral on both RT sources, which is the expected +result for a boundary continuity fix: the anchors are inside the range and the candidates it +moves are the ones the old line had placed outside the acquisition. + +Row and q-value unit of every count above: `peptides.tsv` rows selected on the stripped-peptide +`peptide_q_value` at 1% (one row per stripped sequence here, since `native_tda` was run with +the default `compete.group_by = peptidoform_charge` and the report deduplicates on the winning +row), PSM counts on the pooled `q_value`. diff --git a/rust/mumdia/crates/mumdia-core/src/rejection.rs b/rust/mumdia/crates/mumdia-core/src/rejection.rs index 06c1974..943e8c5 100644 --- a/rust/mumdia/crates/mumdia-core/src/rejection.rs +++ b/rust/mumdia/crates/mumdia-core/src/rejection.rs @@ -7,8 +7,15 @@ //! DIA-NN-only precursor first lost?" without conflating later stages. //! //! The serialized spelling is SCREAMING_SNAKE_CASE and matches the reason strings -//! in the specification exactly (e.g. `NO_PEAK_GROUP`). Use [`RejectionReason::code`] -//! for the stable string written to Parquet/JSON (no serde round-trip cost). +//! in the specification exactly (e.g. `NO_FRAGMENT_TRACES`). Use +//! [`RejectionReason::code`] for the stable string written to Parquet/JSON (no serde +//! round-trip cost). +//! +//! `DID_NOT_SURVIVE_EXTRACTION` was `NO_PEAK_GROUP` until docs/29 #16. The audit +//! assigns it to every candidate that has traces in no extracted row, and `extract` +//! does not write the per-candidate table that would separate presence failures, +//! matched-fraction failures and gate rejections, so the old name claimed a cause the +//! audit cannot see. The name now says what is known. use serde::{Deserialize, Serialize}; @@ -29,7 +36,9 @@ pub enum RejectionReason { CandidateCapReached, // --- extraction + peak formation (Stages C, D) --- NoFragmentTraces, - NoPeakGroup, + /// Extracted no accepted row, for any extraction-side reason; see the module + /// documentation for why this is not called a peak-group failure. + DidNotSurviveExtraction, // --- peak / peptide ranking (Stage E) --- PeakNotSelected, // --- competition (Stage G) --- @@ -59,7 +68,7 @@ impl RejectionReason { RtPruned => "RT_PRUNED", CandidateCapReached => "CANDIDATE_CAP_REACHED", NoFragmentTraces => "NO_FRAGMENT_TRACES", - NoPeakGroup => "NO_PEAK_GROUP", + DidNotSurviveExtraction => "DID_NOT_SURVIVE_EXTRACTION", PeakNotSelected => "PEAK_NOT_SELECTED", OutcompetedByTarget => "OUTCOMPETED_BY_TARGET", OutcompetedByDecoy => "OUTCOMPETED_BY_DECOY", @@ -85,7 +94,7 @@ impl RejectionReason { RtPruned => 6, CandidateCapReached => 7, NoFragmentTraces => 8, - NoPeakGroup => 9, + DidNotSurviveExtraction => 9, PeakNotSelected => 10, OutcompetedByTarget => 11, OutcompetedByDecoy => 12, @@ -119,7 +128,7 @@ mod tests { #[test] fn codes_match_spec_strings() { - assert_eq!(NoPeakGroup.code(), "NO_PEAK_GROUP"); + assert_eq!(DidNotSurviveExtraction.code(), "DID_NOT_SURVIVE_EXTRACTION"); assert_eq!(PeakNotSelected.code(), "PEAK_NOT_SELECTED"); assert_eq!(OutcompetedByDecoy.code(), "OUTCOMPETED_BY_DECOY"); assert_eq!(Reported.code(), "REPORTED"); @@ -136,8 +145,14 @@ mod tests { #[test] fn earliest_keeps_smaller_stage() { // an extraction loss precedes an FDR loss - assert_eq!(NoPeakGroup.earliest(FailedPrecursorFdr), NoPeakGroup); - assert_eq!(FailedPrecursorFdr.earliest(NoPeakGroup), NoPeakGroup); + assert_eq!( + DidNotSurviveExtraction.earliest(FailedPrecursorFdr), + DidNotSurviveExtraction + ); + assert_eq!( + FailedPrecursorFdr.earliest(DidNotSurviveExtraction), + DidNotSurviveExtraction + ); // Reported never wins against a real rejection assert_eq!(Reported.earliest(NoFragmentTraces), NoFragmentTraces); assert_eq!(NoFragmentTraces.earliest(Reported), NoFragmentTraces); @@ -145,7 +160,7 @@ mod tests { #[test] fn is_rejection_flags_only_losses() { - assert!(NoPeakGroup.is_rejection()); + assert!(DidNotSurviveExtraction.is_rejection()); assert!(!Reported.is_rejection()); } @@ -153,8 +168,8 @@ mod tests { fn stage_order_is_monotone_ladder() { // ladder ordering across the major stages assert!(PeptideNotGenerated.stage_order() < NoFragmentTraces.stage_order()); - assert!(NoFragmentTraces.stage_order() < NoPeakGroup.stage_order()); - assert!(NoPeakGroup.stage_order() < PeakNotSelected.stage_order()); + assert!(NoFragmentTraces.stage_order() < DidNotSurviveExtraction.stage_order()); + assert!(DidNotSurviveExtraction.stage_order() < PeakNotSelected.stage_order()); assert!(PeakNotSelected.stage_order() < OutcompetedByTarget.stage_order()); assert!(OutcompetedByTarget.stage_order() < FailedPrecursorFdr.stage_order()); assert!(FailedPrecursorFdr.stage_order() < RemovedDuringReporting.stage_order()); diff --git a/rust/mumdia/crates/mumdia/src/calibrate.rs b/rust/mumdia/crates/mumdia/src/calibrate.rs index 8594d96..4111157 100644 --- a/rust/mumdia/crates/mumdia/src/calibrate.rs +++ b/rust/mumdia/crates/mumdia/src/calibrate.rs @@ -28,12 +28,27 @@ pub fn linear_fit(xs: &[f64], ys: &[f64]) -> (f64, f64) { } /// A LOESS smoother evaluated on a precomputed grid for fast bulk application. +/// +/// Outside the training range the smoother continues the local fit at the nearest +/// boundary: the value at the boundary grid point and the slope of the local line +/// fitted there. It used to switch to the global least-squares line the moment `x` +/// left the grid, and the two models need not agree at the boundary: on +/// `y = 200 + 10 x^2` over `x in [0, 9.9]` (span 0.3) the prediction jumped from +/// 193.4 at `x = 1e-6` to 38.3 at `x = 0`, and from 1173.5 to 1018.4 at the top, +/// discontinuities of about 155 s in retention-time units that misplaced +/// gradient-edge peptides relative to a narrow extraction window (docs/29 #10). pub struct Loess { grid_x: Vec, grid_y: Vec, - // linear fallback for x outside the training range + /// The global least-squares line: the fit for fewer than four points, and the + /// value a degenerate local window (no spread) falls back to. Not the + /// extrapolation model. slope: f64, intercept: f64, + /// Slope of the local line at the first and at the last grid point; extrapolation + /// continues these from the boundary values, so `predict` is continuous there. + lo_slope: f64, + hi_slope: f64, } impl Loess { @@ -62,6 +77,8 @@ impl Loess { grid_y, slope, intercept, + lo_slope: slope, + hi_slope: slope, }; } let k = ((span * n as f64).ceil() as usize).clamp(3, n); @@ -69,28 +86,44 @@ impl Loess { let gn = grid_n.max(2); let mut grid_x = Vec::with_capacity(gn); let mut grid_y = Vec::with_capacity(gn); + let (mut lo_slope, mut hi_slope) = (slope, slope); for g in 0..gn { let x = lo + (hi - lo) * g as f64 / (gn - 1) as f64; + let (y, b) = local_linear(&sx, &sy, x, k, slope, intercept); grid_x.push(x); - grid_y.push(local_linear(&sx, &sy, x, k, slope, intercept)); + grid_y.push(y); + if g == 0 { + lo_slope = b; + } + if g == gn - 1 { + hi_slope = b; + } } Loess { grid_x, grid_y, slope, intercept, + lo_slope, + hi_slope, } } - /// Predict at x by linear interpolation over the grid; linear extrapolation - /// outside the training range. + /// Predict at x by linear interpolation over the grid. Outside the training range + /// the local fit at the nearest boundary is continued (its value there and its + /// slope), so the prediction is continuous across the boundary; see the type's + /// documentation for what the previous global-line switch did. pub fn predict(&self, x: f64) -> f64 { let g = &self.grid_x; - if x <= g[0] || g.len() < 2 { + if g.len() < 2 { return self.slope * x + self.intercept; } - if x >= g[g.len() - 1] { - return self.slope * x + self.intercept; + if x <= g[0] { + return self.grid_y[0] + self.lo_slope * (x - g[0]); + } + let last = g.len() - 1; + if x >= g[last] { + return self.grid_y[last] + self.hi_slope * (x - g[last]); } let j = g.partition_point(|&gx| gx < x); let (x0, x1) = (g[j - 1], g[j]); @@ -103,8 +136,16 @@ impl Loess { } /// Weighted local linear regression at `x0` over the `k` nearest points -/// (tricubic weights). Falls back to the global line if degenerate. -fn local_linear(sx: &[f64], sy: &[f64], x0: f64, k: usize, slope: f64, intercept: f64) -> f64 { +/// (tricubic weights): the fitted value at `x0` and the slope of the local line. +/// Falls back to the global line (value and slope) if degenerate. +fn local_linear( + sx: &[f64], + sy: &[f64], + x0: f64, + k: usize, + slope: f64, + intercept: f64, +) -> (f64, f64) { let n = sx.len(); // find window of k nearest by walking from the insertion point let mut lo = sx.partition_point(|&x| x < x0); @@ -145,11 +186,11 @@ fn local_linear(sx: &[f64], sy: &[f64], x0: f64, k: usize, slope: f64, intercept } let denom = sw * swxx - swx * swx; if sw < 1e-12 || denom.abs() < 1e-12 { - return slope * x0 + intercept; + return (slope * x0 + intercept, slope); } let b = (sw * swxy - swx * swy) / denom; let a = (swy - b * swx) / sw; - a + b * x0 + (a + b * x0, b) } /// The p-th percentile (0..1) of the values (a copy is sorted). @@ -184,6 +225,38 @@ mod tests { assert!((lo.predict(5.0) - 25.0).abs() < 2.0, "{}", lo.predict(5.0)); } + #[test] + fn loess_is_continuous_at_both_training_boundaries() { + // The docs/29 #10 reproduction: a curved map whose global line disagrees with + // the local fit at both ends. The old predictor jumped by about 155 at each + // boundary; the continued local fit must not. + let xs: Vec = (0..100).map(|i| i as f64 / 10.0).collect(); + let ys: Vec = xs.iter().map(|x| 200.0 + 10.0 * x * x).collect(); + let lo = Loess::fit(&xs, &ys, 0.3, 50); + let eps = 1e-6; + let (a0, a1) = (lo.predict(0.0), lo.predict(eps)); + assert!((a0 - a1).abs() < 1e-3, "lower boundary jump: {a0} vs {a1}"); + let (b0, b1) = (lo.predict(9.9), lo.predict(9.9 - eps)); + assert!((b0 - b1).abs() < 1e-3, "upper boundary jump: {b0} vs {b1}"); + // Extrapolation follows the LOCAL slope, not the global line's (99 here): the + // curve is nearly flat at the low end and steep (about 200) at the high end. + let low_step = lo.predict(0.0) - lo.predict(-1.0); + assert!( + (0.0..50.0).contains(&low_step), + "low-end extrapolation should be nearly flat, got a step of {low_step}" + ); + let high_step = lo.predict(10.9) - lo.predict(9.9); + assert!( + high_step > 120.0, + "high-end extrapolation should follow the steep local slope, got {high_step}" + ); + // Far from the data the extrapolation is a line through the boundary value. + assert!( + (lo.predict(-2.0) - (lo.predict(0.0) - 2.0 * low_step)).abs() < 1e-9, + "extrapolation must be linear in x" + ); + } + #[test] fn percentile_basic() { let v: Vec = (0..=100).map(|i| i as f64).collect(); diff --git a/rust/mumdia/crates/mumdia/src/stages/audit.rs b/rust/mumdia/crates/mumdia/src/stages/audit.rs index 10b9e36..2ab512f 100644 --- a/rust/mumdia/crates/mumdia/src/stages/audit.rs +++ b/rust/mumdia/crates/mumdia/src/stages/audit.rs @@ -35,7 +35,9 @@ pub struct AuditParams<'a> { pub scored: &'a str, /// Output `candidate_audit.parquet`. pub out: &'a str, - /// Precursor q-value threshold for `passed_precursor_fdr` / `reported`. + /// Precursor q-value threshold for `passed_precursor_fdr` / `reported`, applied to + /// the scored table's `precursor_q` (the PSM `q_value` only on an older table + /// without that column; `.metrics.json` records which, as `q_unit`). pub q_threshold: f64, /// Run identifier stamped on every row. pub run_id: &'a str, @@ -85,8 +87,36 @@ pub fn run(p: AuditParams) -> Result { .into_iter() .collect(); let scored_t = TableFile::open(p.scored)?; + // One run only. Every survivor set and q lookup here is keyed by candidate_id, so a + // pooled table (several `source` values) would overwrite one run's q with another's + // and attribute the last run's fate to every run (docs/29 #16). `run-experiment` + // writes per-run split tables; audit those. + if let Some(n_sources) = crate::stages::quant::pooled_source_count(&scored_t, p.scored)? { + if n_sources > 1 { + anyhow::bail!( + "audit: {} is a pooled scored table with {n_sources} sources; the audit is \ + per run and keys on candidate_id, so pass one run's split table \ + (//scored.parquet)", + p.scored + ); + } + } let scored_cid = scored_t.u32("candidate_id")?; - let scored_q = scored_t.f64("q_value")?; + // The gate is named `passed_precursor_fdr`, so it reads the precursor q. It read + // the PSM `q_value` (docs/29 #16), a different unit: a PSM can pass at 1% while + // its precursor group does not, and the other way round. Older scored tables have + // no `precursor_q`; there the PSM q is used and the metrics say so. + let (scored_q, q_unit) = match scored_t.f64("precursor_q") { + Ok(v) => (v, "precursor_q"), + Err(_) => { + tracing::warn!( + scored = p.scored, + "audit: no `precursor_q` column; the precursor gate falls back to the PSM \ + q_value, which is not the same unit" + ); + (scored_t.f64("q_value")?, "q_value") + } + }; // peptide-level q is optional (only present in some scored schemas). let scored_pep_q = scored_t.f64("peptide_q_value").ok(); let mut q_by_cid: HashMap = HashMap::with_capacity(scored_cid.len()); @@ -135,15 +165,16 @@ pub fn run(p: AuditParams) -> Result { // Earliest rejection reason along the ladder. let reason: RejectionReason = if !traces { - // Extraction produced no accepted peak for this candidate. Refine with - // the in-extract audit sidecar if present; otherwise the generic bucket. + // Extraction produced no accepted row for this candidate. Refine with the + // in-extract audit sidecar if present; otherwise the generic bucket, which + // says only that the candidate did not survive extraction. match extract_reasons.get(&c).map(String::as_str) { Some("NO_FRAGMENT_TRACES") => RejectionReason::NoFragmentTraces, Some("NO_VALID_FRAGMENTS") => RejectionReason::NoValidFragments, Some("PEAK_NOT_SELECTED") => RejectionReason::PeakNotSelected, Some("RT_PRUNED") => RejectionReason::RtPruned, Some("WRONG_ISOLATION_WINDOW") => RejectionReason::WrongIsolationWindow, - _ => RejectionReason::NoPeakGroup, + _ => RejectionReason::DidNotSurviveExtraction, } } else if !variant { if is_decoy { @@ -207,6 +238,7 @@ pub fn run(p: AuditParams) -> Result { let metrics = json!({ "run_id": p.run_id, "q_threshold": p.q_threshold, + "q_unit": q_unit, "search_space": n, "extracted": n_extracted, "competed": n_competed, @@ -281,9 +313,9 @@ mod tests { // 1 target -> reported (extract+compete+scored q<=0.01) // 2 target -> failed precursor (scored q=0.5) // 3 target -> outcompeted (extract yes, compete no) - // 4 target -> no peak group (not extracted) + // 4 target -> did not survive extraction (not extracted) // 5 decoy -> outcompeted decoy (extract yes, compete no) - // 6 decoy -> no peak group (not extracted) + // 6 decoy -> did not survive extraction (not extracted) let lib = tmp("lib.parquet"); let psms = tmp("psms.parquet"); let comp = tmp("comp.parquet"); @@ -332,9 +364,96 @@ mod tests { assert_eq!(by[&2].0, "FAILED_PRECURSOR_FDR"); assert!(!by[&2].1); assert_eq!(by[&3].0, "OUTCOMPETED_BY_TARGET"); - assert_eq!(by[&4].0, "NO_PEAK_GROUP"); + assert_eq!(by[&4].0, "DID_NOT_SURVIVE_EXTRACTION"); assert_eq!(by[&5].0, "OUTCOMPETED_BY_DECOY"); - assert_eq!(by[&6].0, "NO_PEAK_GROUP"); + assert_eq!(by[&6].0, "DID_NOT_SURVIVE_EXTRACTION"); + } + + #[test] + fn the_precursor_gate_reads_precursor_q_not_psm_q() { + // Two candidates whose PSM and precursor q lie on opposite sides of 1%. The + // label `passed_precursor_fdr` has to follow `precursor_q` (docs/29 #16). + let lib = tmp("lib_q.parquet"); + let psms = tmp("psms_q.parquet"); + let comp = tmp("comp_q.parquet"); + let scored = tmp("scored_q.parquet"); + let out = tmp("audit_q.parquet"); + write_lib(&lib, &[1, 2], &["target", "target"]); + write_cid_only(&psms, &[1, 2]); + write_cid_only(&comp, &[1, 2]); + write_table( + &scored, + vec![ + Col::U32("candidate_id".into(), vec![1, 2]), + Col::F64("q_value".into(), vec![0.001, 0.5]), + Col::F64("precursor_q".into(), vec![0.5, 0.001]), + ], + ) + .unwrap(); + run(AuditParams { + library_precursors: &lib, + psms: &psms, + competed: &comp, + scored: &scored, + out: &out, + q_threshold: 0.01, + run_id: "t", + entrapment_substr: "", + }) + .unwrap(); + let a = TableFile::open(&out).unwrap(); + let cid = a.u32("precursor_id").unwrap(); + let reason = a.str("rejection_reason").unwrap(); + let passed = a.bool("passed_precursor_fdr").unwrap(); + let by: std::collections::HashMap = cid + .iter() + .cloned() + .zip(reason.into_iter().zip(passed)) + .collect(); + // PSM q 0.001 but precursor q 0.5: fails the precursor gate. + assert_eq!(by[&1], ("FAILED_PRECURSOR_FDR".to_string(), false)); + // PSM q 0.5 but precursor q 0.001: passes it. + assert_eq!(by[&2], ("REPORTED".to_string(), true)); + let m: serde_json::Value = + mumdia_io::json::read_json(&format!("{out}.metrics.json")).unwrap(); + assert_eq!(m["q_unit"], "precursor_q"); + } + + #[test] + fn a_pooled_scored_table_is_refused() { + // Keyed by candidate_id alone, a two-source table would let the second run's q + // overwrite the first's. The experiment writes split tables; audit those. + let lib = tmp("lib_pool.parquet"); + let psms = tmp("psms_pool.parquet"); + let comp = tmp("comp_pool.parquet"); + let scored = tmp("scored_pool.parquet"); + let out = tmp("audit_pool.parquet"); + write_lib(&lib, &[1], &["target"]); + write_cid_only(&psms, &[1]); + write_cid_only(&comp, &[1]); + write_table( + &scored, + vec![ + Col::U32("candidate_id".into(), vec![1, 1]), + Col::U32("source".into(), vec![0, 1]), + Col::F64("q_value".into(), vec![0.001, 0.5]), + Col::F64("precursor_q".into(), vec![0.001, 0.5]), + ], + ) + .unwrap(); + let e = run(AuditParams { + library_precursors: &lib, + psms: &psms, + competed: &comp, + scored: &scored, + out: &out, + q_threshold: 0.01, + run_id: "t", + entrapment_substr: "", + }) + .unwrap_err() + .to_string(); + assert!(e.contains("pooled") && e.contains("2 sources"), "{e}"); } #[test] diff --git a/rust/mumdia/crates/mumdia/src/stages/quant.rs b/rust/mumdia/crates/mumdia/src/stages/quant.rs index 67dc8e8..30a45c2 100644 --- a/rust/mumdia/crates/mumdia/src/stages/quant.rs +++ b/rust/mumdia/crates/mumdia/src/stages/quant.rs @@ -2198,7 +2198,7 @@ mod tests { /// ran on the engine's own output: a pooled table quantified against one run's /// chromatograms produced one identical row per run without a word (docs/29 #1). A /// present column that is neither u32 nor i32 is an error now, not a shrug. -fn pooled_source_count(ps: &TableFile, path: &str) -> Result> { +pub(crate) fn pooled_source_count(ps: &TableFile, path: &str) -> Result> { if !ps.has_column("source") { return Ok(None); } @@ -2213,8 +2213,8 @@ fn pooled_source_count(ps: &TableFile, path: &str) -> Result> { .collect::>() .len(), Err(_) => anyhow::bail!( - "quant: column `source` in {path} is present but neither u32 (what rescore \ - writes) nor i32: {u32_err:#}" + "column `source` in {path} is present but neither u32 (what rescore writes) \ + nor i32: {u32_err:#}" ), }, }; diff --git a/rust/mumdia/crates/mumdia/src/stages/report.rs b/rust/mumdia/crates/mumdia/src/stages/report.rs index b705abb..091aeb3 100644 --- a/rust/mumdia/crates/mumdia/src/stages/report.rs +++ b/rust/mumdia/crates/mumdia/src/stages/report.rs @@ -71,6 +71,28 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { }; let accepted = |i: usize, q: &[f64]| label[i] == "target" && (is_transferred[i] || q[i] <= p.q_threshold); + // The acceptance basis is exported with every row (`is_transferred`, `transfer_q`), + // because the rule above is otherwise invisible in the TSV: a transferred row keeps + // its grouped q, usually 1.0, and a tighter threshold does not revoke a transfer that + // already passed `mbr.q_transfer` (docs/29 #19). A transfer is a peptide-level + // acceptance and not protein-group confidence; a group admitted through one carries + // the flag so a reader can tell. + let transfer_q: Vec = match t.f64("transfer_q") { + Ok(v) => v, + Err(_) => vec![f64::NAN; n], + }; + let transfer_cells = |i: usize| -> String { + if is_transferred[i] { + let q = if transfer_q[i].is_nan() { + String::new() + } else { + format!("{:.6}", transfer_q[i]) + }; + format!("\ttrue\t{q}") + } else { + "\tfalse\t".to_string() + } + }; let pep_quant: HashMap<(String, i32), f64> = match p.peptide_quant { Some(path) => { @@ -111,7 +133,8 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { // sequence; the header and the returned count reflect that. writeln!( w, - "precursor\tstripped_sequence\tcharge\tprotein\tq_value\tscore\tquantity" + "precursor\tstripped_sequence\tcharge\tprotein\tq_value\tscore\tquantity\t\ + is_transferred\ttransfer_q" )?; let mut npep = 0u64; for &i in &order { @@ -126,14 +149,15 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { let qv = pep_quant.get(&key).copied().unwrap_or(f64::NAN); writeln!( w, - "{}\t{}\t{}\t{}\t{:.6}\t{:.4}\t{}", + "{}\t{}\t{}\t{}\t{:.6}\t{:.4}\t{}{}", pform[i], strip(&pform[i]), charge[i], protein[i], pep_q[i], score[i], - qcell(qv) + qcell(qv), + transfer_cells(i) )?; npep += 1; } @@ -147,7 +171,10 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { let mut pseen: HashSet = HashSet::new(); let prot_target = mumdia_io::table::AtomicPath::new(p.out_proteins)?; let mut w2 = std::io::BufWriter::new(std::fs::File::create(prot_target.tmp())?); - writeln!(w2, "protein_group\tq_value\tquantity")?; + writeln!( + w2, + "protein_group\tq_value\tquantity\tis_transferred\ttransfer_q" + )?; let mut nprot = 0u64; for &i in &porder { if !accepted(i, &pg_q) || pg[i].is_empty() { @@ -157,7 +184,14 @@ pub fn run(p: ReportParams) -> Result<(u64, u64)> { continue; } let qv = prot_quant.get(&pg[i]).copied().unwrap_or(f64::NAN); - writeln!(w2, "{}\t{:.6}\t{}", pg[i], pg_q[i], qcell(qv))?; + writeln!( + w2, + "{}\t{:.6}\t{}{}", + pg[i], + pg_q[i], + qcell(qv), + transfer_cells(i) + )?; nprot += 1; } w2.flush()?; @@ -229,6 +263,28 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { }; let accepted = |i: usize, q: &[f64]| label[i] == "target" && (is_transferred[i] || q[i] <= p.q_threshold); + // The acceptance basis is exported with every row (`is_transferred`, `transfer_q`), + // because the rule above is otherwise invisible in the TSV: a transferred row keeps + // its grouped q, usually 1.0, and a tighter threshold does not revoke a transfer that + // already passed `mbr.q_transfer` (docs/29 #19). A transfer is a peptide-level + // acceptance and not protein-group confidence; a group admitted through one carries + // the flag so a reader can tell. + let transfer_q: Vec = match t.f64("transfer_q") { + Ok(v) => v, + Err(_) => vec![f64::NAN; n], + }; + let transfer_cells = |i: usize| -> String { + if is_transferred[i] { + let q = if transfer_q[i].is_nan() { + String::new() + } else { + format!("{:.6}", transfer_q[i]) + }; + format!("\ttrue\t{q}") + } else { + "\tfalse\t".to_string() + } + }; // Runs in which each precursor / protein group was identified on its own per-run FDR. let mut pep_runs: HashMap<(String, i32), Vec> = HashMap::new(); @@ -305,6 +361,7 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { write!(w, "\tquantity_{name}")?; } } + write!(w, "\tis_transferred\ttransfer_q")?; writeln!(w)?; let mut npep = 0u64; for &i in &order { @@ -330,6 +387,7 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { for m in &pep_quant { write!(w, "\t{}", qcell(m.get(&key).copied().unwrap_or(f64::NAN)))?; } + write!(w, "{}", transfer_cells(i))?; writeln!(w)?; npep += 1; } @@ -349,6 +407,7 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { write!(w2, "\tlfq_{name}")?; } } + write!(w2, "\tis_transferred\ttransfer_q")?; writeln!(w2)?; let mut nprot = 0u64; for &i in &porder { @@ -371,6 +430,7 @@ pub fn run_experiment(p: ExperimentReportParams) -> Result<(u64, u64)> { write!(w2, "\t{}", qcell(*x))?; } } + write!(w2, "{}", transfer_cells(i))?; writeln!(w2)?; nprot += 1; } @@ -510,7 +570,7 @@ mod tests { let lines: Vec<&str> = text.lines().collect(); assert_eq!( lines[0], - "precursor\tstripped_sequence\tcharge\tprotein\tq_value\tscore\tn_runs\tquantity_a\tquantity_b" + "precursor\tstripped_sequence\tcharge\tprotein\tq_value\tscore\tn_runs\tquantity_a\tquantity_b\tis_transferred\ttransfer_q" ); assert_eq!(lines.len(), 3, "header plus two precursors:\n{text}"); assert!( @@ -531,7 +591,10 @@ mod tests { let prot = std::fs::read_to_string(&proteins).unwrap(); let plines: Vec<&str> = prot.lines().collect(); - assert_eq!(plines[0], "protein_group\tq_value\tn_runs\tlfq_a\tlfq_b"); + assert_eq!( + plines[0], + "protein_group\tq_value\tn_runs\tlfq_a\tlfq_b\tis_transferred\ttransfer_q" + ); assert_eq!( plines.len(), 2, @@ -642,11 +705,13 @@ mod tests { "stripped sequence missing:\n{text}" ); // Quantity is empty when no quant table was supplied, not zero: absence of - // a measurement is not a measurement of zero. + // a measurement is not a measurement of zero. Nothing here is a transfer. let data_line = text.lines().nth(1).unwrap(); - assert!( - data_line.ends_with('\t'), - "expected an empty quantity cell: {data_line:?}" + let cells: Vec<&str> = data_line.split('\t').collect(); + assert_eq!( + (cells[6], cells[7], cells[8]), + ("", "false", ""), + "expected an empty quantity cell and no transfer: {data_line:?}" ); let prot = std::fs::read_to_string(&proteins).unwrap(); @@ -656,6 +721,121 @@ mod tests { assert_eq!(prot.lines().count(), 2, "header plus one group:\n{prot}"); } + /// The same table after `mumdia mbr`: the above-threshold target is an accepted + /// transfer (`is_transferred`, with the transfer q the worker accepted it at). + fn scored_table_with_transfer() -> String { + let path = tmp("scored_mbr.parquet"); + write_table( + &path, + vec![ + Col::Str( + "peptidoform".into(), + vec![ + "PEPTIDEK".into(), + "M[Oxidation]EGVDGHK".into(), + "DECOY_KEDITPEP".into(), + "LATEPEPTIDEK".into(), + ], + ), + Col::I32("charge".into(), vec![2, 3, 2, 2]), + Col::Str( + "protein".into(), + vec!["P1".into(), "P2".into(), "DECOY_P1".into(), "P1".into()], + ), + Col::Str( + "label".into(), + vec![ + "target".into(), + "target".into(), + "decoy".into(), + "target".into(), + ], + ), + Col::F64("peptide_q_value".into(), vec![0.001, 0.005, 0.0001, 0.5]), + Col::Str( + "protein_group".into(), + vec!["PG1".into(), "PG2".into(), "DECOY_PG1".into(), "".into()], + ), + Col::F64("pg_q_value".into(), vec![0.002, 0.9, 0.0001, 1.0]), + Col::F64("score".into(), vec![3.5, 2.5, 9.9, 0.1]), + Col::Bool("is_transferred".into(), vec![false, false, false, true]), + Col::F64( + "transfer_q".into(), + vec![f64::NAN, f64::NAN, f64::NAN, 0.004], + ), + ], + ) + .unwrap(); + path + } + + #[test] + fn report_exports_the_transfer_basis_and_a_tighter_threshold_does_not_revoke_it() { + // docs/29 #19: the acceptance rule admits any target flagged `is_transferred` + // whatever its grouped q, and the TSV has to show that basis rather than a + // `q_value` of 0.5 next to a 1% threshold with no explanation. + let scored = scored_table_with_transfer(); + let out = tmp("report_mbr"); + std::fs::create_dir_all(&out).unwrap(); + let peptides = format!("{out}/peptides.tsv"); + let proteins = format!("{out}/proteins.tsv"); + let report = |q: f64| { + run(ReportParams { + scored: &scored, + peptide_quant: None, + protein_quant: None, + out_peptides: &peptides, + out_proteins: &proteins, + q_threshold: q, + }) + .unwrap() + }; + + let (n_pep, _) = report(0.01); + assert_eq!(n_pep, 3, "two confident targets plus the transfer"); + let text = std::fs::read_to_string(&peptides).unwrap(); + let header = text.lines().next().unwrap(); + assert!( + header.ends_with("\tquantity\tis_transferred\ttransfer_q"), + "{header}" + ); + let late = text + .lines() + .find(|l| l.starts_with("LATEPEPTIDEK")) + .expect("the transfer is reported"); + let cells: Vec<&str> = late.split('\t').collect(); + // The grouped q is preserved as it is (0.5), and the basis is next to it. + assert_eq!( + (cells[4], cells[7], cells[8]), + ("0.500000", "true", "0.004000") + ); + let first = text.lines().nth(1).unwrap(); + let c: Vec<&str> = first.split('\t').collect(); + assert_eq!( + (c[7], c[8]), + ("false", ""), + "a scored row is not a transfer: {first}" + ); + + // Protein rows carry the same two columns; PG1 was admitted on its own q. + let prot = std::fs::read_to_string(&proteins).unwrap(); + let plines: Vec<&str> = prot.lines().collect(); + assert_eq!( + plines[0], + "protein_group\tq_value\tquantity\tis_transferred\ttransfer_q" + ); + let pg: Vec<&str> = plines[1].split('\t').collect(); + assert_eq!((pg[0], pg[3], pg[4]), ("PG1", "false", "")); + + // Tighter threshold: the 0.005 target drops, the transfer stays. This is the + // rule the columns exist to make visible, not a new one. + let (n_pep, _) = report(0.001); + assert_eq!(n_pep, 2); + let text = std::fs::read_to_string(&peptides).unwrap(); + assert!(text.contains("LATEPEPTIDEK\t"), "{text}"); + assert!(!text.contains("M[Oxidation]EGVDGHK"), "{text}"); + } + #[test] fn report_writes_header_only_when_nothing_passes() { // A run that identifies nothing at the threshold still has to produce the @@ -715,7 +895,8 @@ mod tests { let mut quantified = 0; let mut empty = 0; for line in text.lines().skip(1) { - let cell = line.rsplit('\t').next().unwrap(); + // The quantity column, by position: the transfer columns follow it. + let cell = line.split('\t').nth(6).unwrap(); if cell.is_empty() { empty += 1; } else { diff --git a/rust/mumdia/crates/mumdia/src/stages/run_experiment.rs b/rust/mumdia/crates/mumdia/src/stages/run_experiment.rs index 69282f5..ca78064 100644 --- a/rust/mumdia/crates/mumdia/src/stages/run_experiment.rs +++ b/rust/mumdia/crates/mumdia/src/stages/run_experiment.rs @@ -368,6 +368,34 @@ pub fn run(p: RunExperimentParams) -> Result<()> { std::fs::create_dir_all(p.out_dir).ok(); let d = |name: &str| format!("{}/{}", p.out_dir, name); let n_runs = p.mzmls.len(); + + // Provenance: the identity of the code, of the configuration and of the INPUTS, + // hashed now, before anything reads them for compute (docs/29 #15). Hashing at the + // end recorded whatever bytes were on disk after a multi-hour experiment, which is + // not necessarily what the search read; the single-run orchestrator has always + // hashed first, and the two now agree. + let mut prov = Manifest::new(cfg.canonical_json(), ch.clone()); + for (i, m) in p.mzmls.iter().enumerate() { + if let (Ok(bytes), Ok(hash)) = ( + std::fs::metadata(m).map(|x| x.len()), + mumdia_io::hash::blake3_file(m), + ) { + prov.record_input(&format!("mzml[{i}]"), m, bytes, hash); + } + } + for (role, path) in [ + ("fasta", p.fasta), + ("lib_precursors", p.lib_precursors), + ("lib_fragments", p.lib_fragments), + ] { + let Some(path) = path else { continue }; + if let (Ok(bytes), Ok(hash)) = ( + std::fs::metadata(path).map(|x| x.len()), + mumdia_io::hash::blake3_file(path), + ) { + prov.record_input(role, path, bytes, hash); + } + } // Reject a bad --run-names rather than silently substituting r0..rN-1. // // The old `_ =>` arm swallowed any count mismatch with no warning, and accepted @@ -603,6 +631,20 @@ pub fn run(p: RunExperimentParams) -> Result<()> { cfg: &cfg.rescore, config_hash: &ch, })?; + // The classifier that actually ran, from the rescore artifact report: the source of + // truth, since the configured enum can differ from it under a compatibility path. + let rescore_report: mumdia_io::report::ArtifactReport = + mumdia_io::json::read_json(&format!("{scored_combined}.report.json"))?; + let actual_rescorer = rescore_report + .params + .get("classifier") + .and_then(|v| v.as_str()) + .unwrap_or("unknown") + .to_string(); + let actual_rescorer_model = rescore_report + .model_identity + .clone() + .unwrap_or_else(|| actual_rescorer.clone()); // --- optional rescuable-tier MBR transfer --- let scored_for_quant = if cfg.mbr.strategy != mumdia_core::config::MbrStrategy::None { @@ -741,7 +783,15 @@ pub fn run(p: RunExperimentParams) -> Result<()> { "runs": names, "scored_combined": scored_combined, "scored_for_quant": scored_for_quant, + "rescorer": actual_rescorer, "mbr": format!("{:?}", cfg.mbr.strategy), + // Per-run quant gates on the pooled q_value whatever the configuration says + // (see the warning above). The configuration hash is of the configuration as + // given, so the substitution has to be recorded here or it is recorded nowhere. + "quant_q_filter": { + "configured": format!("{:?}", cfg.quant.q_filter), + "effective": format!("{:?}", qcfg.q_filter), + }, "lfq": lfq, "peptide_quants": peptide_quants, "report": { @@ -754,7 +804,8 @@ pub fn run(p: RunExperimentParams) -> Result<()> { }, }); // Provenance parity with the single-run manifest: the identity of the code, of - // the inputs, and of every artifact this stage produced. + // the inputs (hashed at the start of the run, see above), and of every artifact + // this stage produced. // // The per-artifact records were the gap. The experiment manifest listed output // PATHS in its `experiment` block and nothing else, so an experiment result had @@ -763,28 +814,44 @@ pub fn run(p: RunExperimentParams) -> Result<()> { // needed to tell whether two experiment outputs are the same data. Files written // by the per-run chains are not covered here (those chains do not thread a shared // manifest); what is covered is everything `run-experiment` itself writes. - let mut prov = Manifest::new(cfg.canonical_json(), ch.clone()); - for (i, m) in p.mzmls.iter().enumerate() { - if let (Ok(bytes), Ok(hash)) = ( - std::fs::metadata(m).map(|x| x.len()), - mumdia_io::hash::blake3_file(m), - ) { - prov.record_input(&format!("mzml[{i}]"), m, bytes, hash); - } - } - for (role, path) in [ - ("fasta", p.fasta), - ("lib_precursors", p.lib_precursors), - ("lib_fragments", p.lib_fragments), - ] { - let Some(path) = path else { continue }; - if let (Ok(bytes), Ok(hash)) = ( - std::fs::metadata(path).map(|x| x.len()), - mumdia_io::hash::blake3_file(path), - ) { - prov.record_input(role, path, bytes, hash); + // + // Model identities reflect the path that produced the downstream artifacts, as + // in the single-run manifest (docs/29 #15): which RT source the library carried, + // which fragment predictor, and the classifier that actually ran. + let library_input = p.lib_precursors.is_some(); + let rt_identity = if cfg.rt_im_train.finetune_deeplc { + if matches!(cfg.experiment.finetune_scope, FinetuneScope::FirstRunOnly) { + "deeplc-finetuned-first-run".to_string() + } else { + "deeplc-finetuned-per-run".to_string() } - } + } else if cfg + .rt_im_train + .repredicts_library_irt(library_input, cfg.predict_frag.deeplc_python.is_some()) + { + "deeplc-base-model".to_string() + } else if library_input { + "imported-library".to_string() + } else { + format!("{:?}", cfg.predict_frag.rt_predictor) + }; + let fragment_identity = if library_input { + "imported-library".to_string() + } else { + format!("{:?}", cfg.predict_frag.predictor) + }; + prov.model_identities + .insert("rt_predictor".into(), rt_identity); + prov.model_identities + .insert("fragment_predictor".into(), fragment_identity); + prov.model_identities + .insert("rescorer".into(), actual_rescorer_model); + prov.model_identities.insert( + "feature_schema_id".into(), + features::feature_schema_id(&features::active_features(cfg.features.set)), + ); + prov.model_identities + .insert("mbr".into(), format!("{:?}", cfg.mbr.strategy)); // Recorded in a fixed order, and every record hashes its file. `Manifest` // stores them in a BTreeMap, so the serialised order is by logical name and // does not depend on this sequence. @@ -841,12 +908,18 @@ pub fn run(p: RunExperimentParams) -> Result<()> { &ch, )?); + // The resolved configuration itself, not only its hash: a hash identifies a + // configuration but cannot replay one (docs/29 #15). let manifest = serde_json::json!({ "mumdia_version": prov.mumdia_version, "git_sha": prov.git_sha, "commit_date": prov.commit_date, "cli_args": prov.cli_args, + "config_hash": prov.config_hash, + "config_json": prov.config_json, + "model_identities": prov.model_identities, "inputs": prov.inputs, + "inputs_hashed_at": "start", "artifacts": prov.artifacts, "experiment": manifest, }); diff --git a/scripts/mbr_worker.py b/scripts/mbr_worker.py index 948fb77..53796ab 100644 --- a/scripts/mbr_worker.py +++ b/scripts/mbr_worker.py @@ -337,6 +337,10 @@ def cos(pa_, pb): if col in full.columns: full[col] = np.minimum(full[col].to_numpy(dtype=float), tq) full["is_transferred"] = is_tr + # The q the transfer was accepted at, on the transferred rows only, so the + # report can print the acceptance basis next to the untouched grouped q + # (docs/29 #19). NaN on every other row: no transfer, no transfer q. + full["transfer_q"] = np.where(is_tr, tq, np.nan) write_engine_parquet(full, a.out_scored) print(f"wrote {a.out_scored} (augmented scored; {int(is_tr.sum())} rows flagged transferred)") diff --git a/tests/python/test_mbr_worker.py b/tests/python/test_mbr_worker.py index 734cab3..a8d6c5d 100644 --- a/tests/python/test_mbr_worker.py +++ b/tests/python/test_mbr_worker.py @@ -334,6 +334,18 @@ def test_m5_sets_is_transferred_on_exactly_the_accepted_rows(mbr_dataset, mbr_re got = np.asarray(after["is_transferred"], dtype=bool) assert np.array_equal(got, expect) + # The acceptance basis travels with the flag: `transfer_q` is the accepted q on + # exactly the flagged rows and NaN everywhere else (docs/29 #19). + tq = np.asarray(after["transfer_q"], dtype=float) + assert np.array_equal(np.isfinite(tq), expect) + tr = read_columns(mbr_result["transferred"]) + accepted_q = {k: q for k, q in zip( + zip((int(c) for c in tr["candidate_id"]), (int(s) for s in tr["source"])), + (float(q) for q in tr["transfer_q"]))} + for k, flagged, q in zip(key, expect, tq): + if flagged: + assert q == accepted_q[k] + def test_m5_touches_only_the_matching_candidate_id_and_source(mbr_dataset, mbr_result): """Only the `(candidate_id, source)` row of an accepted transfer may change.