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
17 changes: 13 additions & 4 deletions docs/07_search_seed.md
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ a candidate that always ranks below `report_psms` gets no row even if it cleared
| `matched_peaks` | i32 | matched-fragment count at the best scan |
| `scan_index` | u32 | `scan_index` of the best-scoring scan |

**`<out>.masscal.json`** (written at `search_seed.rs:217-227`):
**`<out>.masscal.json`** (fitted and serialised by `masscal::MassCal`, `masscal.rs`):

| key | type | meaning |
|---|---|---|
Expand Down Expand Up @@ -232,7 +232,16 @@ two-element average), then sets `tol = max(5.0, 1.5 * P95(|dev - offset|))` usin
`>= 20` survive, re-fit to `(o2, t2, 2)`; otherwise keep the single-pass result.
The second pass rejects random-match outliers so they cannot bias the median.

The result is written to `<out>.masscal.json` (`:217-227`).
The estimator is `masscal::MassCal::fit_from` and the JSON body its `to_json`, both in
`masscal.rs`, because a grouped search fits the very same estimator a second time: each
band's seed would otherwise learn a tolerance from its own ~2,000 deviations and
`seed-pool` would average the bands' scalars, which measured 35% wider and cost 3.4% of the
peptides (`docs/33_window_groups.md` section 4a). Under `groups.window_groups > 1` the stage
is additionally asked (`SearchSeedParams::emit_calibrants`) to write
`<out>.masscal.parquet`: `candidate_id` (library-wide), `scan_index`, `frag_mz`, `ppm`, one
row per deviation, plus the band's best-scoring targets down to
`masscal::CALIBRANT_OFFER_PSMS` so the pooled q rather than the band's own q chooses the
calibrants. An ungrouped run writes no sidecar and is byte-identical.

**7. Write** `seed_psms.parquet` (`write_table`, `:233`) and the `ArtifactReport`
(`:256-274`), then log `psms`, `confident`, `elapsed_ms`.
Expand Down Expand Up @@ -380,8 +389,8 @@ takes `--ms2`, `--lib-precursors`, `--lib-fragments`, `--out`, and
branch in the `if let Some(idx) = fidx` dispatch (`search_seed.rs:63-108`); mirror
the deterministic merge if the new path is parallel.
- **Charge-2 / m/z-binned tolerance.** The current fit produces one global offset
and tolerance. To make them charge- or m/z-dependent, partition `devs` before the
`fit` closure and emit per-bin entries in `masscal.json`, then teach `extract` to
and tolerance. To make them charge- or m/z-dependent, partition `devs` before
`masscal::fit` and emit per-bin entries in `masscal.json`, then teach `extract` to
pick the matching bin.
- **Change the score.** `hyperscore` (`search_seed.rs:413`) is a free function; keep
it monotone in matched-fragment count and observed intensity so the target-decoy q
Expand Down
77 changes: 69 additions & 8 deletions docs/33_window_groups.md
Original file line number Diff line number Diff line change
Expand Up @@ -108,14 +108,75 @@ ids are local.
`seed-pool` (`stages/seed_pool.rs`) reads every band's seed, maps the ids to library-wide
ones, keeps the higher-scoring row where two bands searched the same candidate, re-estimates
`spectrum_q` over the union with the seed's own target-decoy kernel (`fdr::target_decoy_q`),
and writes the run-level `seed_psms.parquet`. Beside it goes `seed_psms.parquet.masscal.json`:
the bands' scalar ppm offsets and learned tolerances combined by calibrant count (`n_dev`).
The optional m/z-dependent grids are not combined, because they need the per-fragment
deviations the seed does not keep; extract then applies the scalar offset. A band with no
calibrants wrote the configured tolerance in place of a learned one (search-seed's failure
branch), and when no band calibrated the pool keeps that tolerance rather than averaging
nothing into a zero, with a warning. It also writes each band's view of the pooled seed,
`groups/gNN/seed_psms_pooled.parquet`: the band's own rows, local ids, pooled q.
and writes the run-level `seed_psms.parquet`. Beside it goes
`seed_psms.parquet.masscal.json`, the run's fragment mass calibration, fitted once over the
bands' calibrant deviations (next subsection). It also writes each band's view of the
pooled seed, `groups/gNN/seed_psms_pooled.parquet`: the band's own rows, local ids,
pooled q.

### 4a. The mass calibration is fitted once, not averaged

The tolerance `search-seed` learns is `1.5 * p95(|dev - median|)` over the ppm deviations of
the matched fragments of its confident target PSMs (`docs/07_search_seed.md`). Combining the
bands' fitted SCALARS is not that estimator, and it is systematically wider. Measured on the
six-file HYE Astral benchmark, the same library and the same retention-time model, 100 bands
against none:

| | unbanded | 100 bands, scalars combined |
|---|---|---|
| `frag_ppm_offset` | -1.8486 | -1.8834 |
| `frag_ppm_sigma` | 8.452 | 11.400 |
| `n_dev` | 181,196 | 200,257 |
| `ppm_residual_mad` | 0.907 | 1.035 |
| candidates accepted by extract | 4,986,153 | 6,609,984 |
| peptides at 1% | 113,860 | 110,006 |
| precursors at 1% | 126,436 | 121,966 |

The offsets agree and the tolerance is 35% wider, extract then accepts 33% more candidates
and the run returns 3.4% fewer peptides. Causal rather than correlated: one band extracted
twice, identical in every input but the calibration file it was handed, accepted 12,414
candidates at 11.40 ppm and 8,791 at 8.45; and an unbanded control on the same adapted
library returned the unbanded arm's 113,860 peptides exactly.

Two things compound. A p95 estimated on one band's ~2,000 deviations has a heavier tail than
the p95 of the union, and a mean of those p95s does not recover it. And each band selects
its calibrants on its OWN `spectrum_q`, which is not the pooled one: on that benchmark it is
looser (106,088 band-confident PSMs against 97,584 pooled).

So the bands write their deviations and the pool fits them, exactly as the retention-time
calibration already uses the pooled anchors under `groups.calibration = global`:

- `search-seed`, asked for it (`SearchSeedParams::emit_calibrants`, set only by the grouped
path), writes `<seed>.masscal.parquet` beside the masscal: `candidate_id` (LIBRARY-WIDE,
the band's local id plus its fragment offset), `scan_index`, `frag_mz` and `ppm`, one row
per matched fragment. 16 B per deviation, a few MB for a whole run.
- `seed-pool` reads them, keeps the deviations whose PSM the POOLED q accepts at
`search_seed.fdr_seed` (and whose scan is the one the pool kept, so an overlap candidate
contributes one PSM's fragments as an ungrouped seed would), and fits
`masscal::MassCal::fit_from` -- the same function `search-seed` calls -- once.
`two_pass_mass_cal` and the `mass_cal_loess` grid work on this path, because both read the
deviations rather than a scalar.
- The band's own `<seed>.masscal.json` is unchanged: still fitted on that band's confident
targets alone, and still what `groups.calibration = per_group` extracts with.

The band's q is not the pooled q in either direction, so the sidecar carries more than the
band's own selection: the band's best-scoring targets down to
`masscal::CALIBRANT_OFFER_PSMS` (2,000) are offered whatever its own q says, and the pooled
q decides. That matters where the band q is STRICTER, which is the `1/T` case above: on the
CI fixture at three bands, every band's own q rejects every one of its targets, so a
strictly q-selected sidecar would be empty in all three. With the offer, the three bands
contribute 437, 401 and 230 deviations, the pool selects all 1,068 of them and fits
`frag_tol_ppm` 5.0 at offset 0.0 -- the same calibration, to the digit, that the ungrouped
search of the same spectra produces. Before this the same run calibrated nothing and
extracted at the configured 20 ppm.

`masscal.json` gains `masscal_source`, which reads `pooled_deviations` or `band_scalars`. A
band directory seeded before the sidecar existed has none, and the pool then combines the
scalars as it always did, naming the missing groups in a warning; if `mass_cal_loess` is set
it warns separately that the grid cannot be recovered from scalars. A band with no
calibrants wrote the configured tolerance in place of a learned one, and when the pooled fit
has fewer than `masscal::MIN_CALIBRANTS` deviations the pool keeps that tolerance rather
than fitting a percentile of a handful of points, with a warning.

`groups.calibration` decides which anchors each band's calibration sees:

Expand Down
1 change: 1 addition & 0 deletions rust/mumdia/crates/mumdia/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ pub mod calibrate;
pub mod fdr;
pub mod groups;
pub mod index;
pub mod masscal;
pub mod matchers;
pub mod memlog;
pub mod peaks;
Expand Down
1 change: 1 addition & 0 deletions rust/mumdia/crates/mumdia/src/main.rs
Original file line number Diff line number Diff line change
Expand Up @@ -1203,6 +1203,7 @@ fn real_main() -> Result<()> {
let ch = mumdia_io::hash::blake3_str(&cfg.canonical_json());
stages::search_seed::run(stages::search_seed::SearchSeedParams {
fragment_offset: None,
emit_calibrants: false,
ms2: &ms2,
library_precursors: &lib_precursors,
library_fragments: &lib_fragments,
Expand Down
Loading
Loading