Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 41 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand Down Expand Up @@ -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
Expand Down
35 changes: 21 additions & 14 deletions CLAUDE.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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<Vec<f64>>`, 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

Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
15 changes: 10 additions & 5 deletions ci/check_smoke.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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")
Expand Down
13 changes: 7 additions & 6 deletions desktop/src-tauri/src/run.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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);
Expand All @@ -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.
Expand Down
2 changes: 1 addition & 1 deletion docs/02_config_and_data_model.md
Original file line number Diff line number Diff line change
Expand Up @@ -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) |
Expand Down
4 changes: 2 additions & 2 deletions docs/04_convert.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
21 changes: 12 additions & 9 deletions docs/08_rt_im_train.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand All @@ -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
Expand Down Expand Up @@ -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. |
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading