From 476f8e90fce4b11a99ab3c1cd681676e7fa1e158 Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Sat, 12 Sep 2026 20:39:56 -0700 Subject: [PATCH 1/2] Score sequence-free libraries with library-wide chemistry fallbacks --- docs/development.md | 20 +- rust/timsquery/src/models/capabilities.rs | 14 +- rust/timsquery/src/models/query_handle.rs | 6 +- rust/timsquery/src/models/target_columns.rs | 21 +- rust/timsquery/src/serde/library_file.rs | 10 +- rust/timsquery/src/traits/fragment_label.rs | 13 ++ rust/timsquery_cli/src/commands.rs | 4 +- .../src/data_sources/reference_library.rs | 207 +++++++++++------- rust/timsseek/src/ml/qvalues.rs | 156 ++++++++++--- rust/timsseek/src/scoring/blocks/lazy.rs | 13 +- rust/timsseek/src/scoring/parquet_writer.rs | 83 ++++++- rust/timsseek/src/scoring/pipeline.rs | 106 +++++++-- rust/timsseek/src/scoring/plan.rs | 72 +++++- rust/timsseek/src/scoring/timings.rs | 2 + rust/timsseek_cli/src/processing.rs | 51 +++++ rust/timsseek_cli/tests/build_library.rs | 2 +- 16 files changed, 614 insertions(+), 166 deletions(-) diff --git a/docs/development.md b/docs/development.md index d8f3a096..551aa762 100644 --- a/docs/development.md +++ b/docs/development.md @@ -22,11 +22,21 @@ and target JSON. Format detection also inspects contents; the `timsseek` and `timsquery_viewer` use the same registry through `ReferenceLibrary`. Every loaded scoring library contains geometry usable for query extraction. -The current scoring bridge requires ion-annotated fragments and reference -fragment intensities; extraction also accepts opaque fragment labels. The -query reader's target-list JSON schemas (`Target` and `ElutionGroupInput` arrays) -currently load geometry without a reference-intensity sidecar. This describes -those reader paths, not a restriction of JSON as an encoding. +Scoring requires retained fragments and aligned reference intensities. mzSpecLib +can supply these without sequence or fragment annotations. Opaque peaks disable +fragment-isotope scores and automatic mass-shift decoy generation library-wide; +supplied decoys remain usable. Without decoys, search scores the full acquisition +RT range and writes raw scores, omitting competition, discriminant-score and +q-value columns (`result_mode=raw` Parquet metadata). It bypasses calibration, +rescoring and q-value filtering. Supplied decoys do not by themselves validate a +decoy strategy for a new analyte class. + +Programmatic `TargetTable::Str` accepts an optional intensity sidecar. Scoring +assigns unique packed unknown keys (currently at most 255 peaks per entry), +preserving the original opaque labels separately; even a string such as `y3` +is not interpreted as chemistry. The query reader's target-list JSON schemas +(`Target` and `ElutionGroupInput` arrays) still supply geometry without reference +intensities, so those reader routes remain extraction-only. Standalone `calib_dash` reads saved `calibration.json`, not a spectral library. diff --git a/rust/timsquery/src/models/capabilities.rs b/rust/timsquery/src/models/capabilities.rs index f88115d3..e8e09d81 100644 --- a/rust/timsquery/src/models/capabilities.rs +++ b/rust/timsquery/src/models/capabilities.rs @@ -10,10 +10,10 @@ pub struct TargetCapabilities { pub decoys: DecoyStrategy, } -/// Runtime reflection of whether this arena's label carries ion chemistry -/// (`FragmentLabel`). `IonAnnot` arenas => `Available`; string-labelled -/// arenas => `Unavailable`. Sequence-operation eligibility is resolved from -/// stored analyte facts by the scoring library. +/// Whether the label representation can carry ion chemistry. `IonAnnot` may +/// still contain unknown placeholders; scoring inspects actual labels across +/// the whole library before enabling operations. Sequence-operation eligibility +/// is independently resolved from stored analyte facts. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum FragmentFeatureState { Available, @@ -126,9 +126,9 @@ impl std::str::FromStr for DecoyPolicy { /// What to do with a peak this reader cannot annotate. /// /// A kept peak lands at the m/z the file measured, where an annotated one lands -/// at the m/z its annotation implies. Nothing downstream can tell the two apart -/// once they are in the arena, so which of them a library is made of is the -/// caller's decision rather than a fallback the reader picks. +/// at the m/z its annotation implies. Unknown keys remain recognizable downstream +/// and cannot establish fragment chemistry. The policy chooses which peaks to +/// retain; operation eligibility is resolved over all retained labels. /// /// Asked only of a library that annotates nothing at all, with /// [`KeepAll`](Self::KeepAll) the exception that asks nothing. mzPAF spells diff --git a/rust/timsquery/src/models/query_handle.rs b/rust/timsquery/src/models/query_handle.rs index 7afceb1b..a33b7852 100644 --- a/rust/timsquery/src/models/query_handle.rs +++ b/rust/timsquery/src/models/query_handle.rs @@ -300,12 +300,12 @@ mod tests { analyte: analyte::Analyte::from_sequence("PEP").as_input(), ..Default::default() }); - // Sealed `IfMissing` with no shipped decoy, so the arena derives ± - // variants -- and variant 1 is where an ion-annotated label would shift. + // Opaque labels disable generation for the whole library. let c = c .seal(DecoyPolicy::IfMissing) .expect("fixture ids are usable"); - let q = Query::new(&c, c.flat_for(first_row(&c), 1)); + assert_eq!(c.variants_per_row(), 1); + let q = Query::new(&c, c.flat_for(first_row(&c), 0)); let frags: Vec<_> = q.iter_fragments_refs().collect(); assert!((frags[0].1 - 300.0).abs() < 1e-9); } diff --git a/rust/timsquery/src/models/target_columns.rs b/rust/timsquery/src/models/target_columns.rs index c40ff117..23ce6dda 100644 --- a/rust/timsquery/src/models/target_columns.rs +++ b/rust/timsquery/src/models/target_columns.rs @@ -287,7 +287,10 @@ impl TargetColumnsBuilder { self.inner.n_fragments() } - pub fn seal(self, decoys: DecoyPolicy) -> Result, TargetBuildError> { + pub fn seal(self, decoys: DecoyPolicy) -> Result, TargetBuildError> + where + L: DecoyShift, + { self.inner.seal(decoys) } } @@ -547,7 +550,10 @@ impl TargetColumns { self.frag_off[tgt] as usize..self.frag_off[tgt + 1] as usize } - fn seal(mut self, decoys: DecoyPolicy) -> Result { + fn seal(mut self, decoys: DecoyPolicy) -> Result + where + L: DecoyShift, + { assert_eq!(self.analytes.len(), self.n_rows()); assert_eq!(self.entry_names.len(), self.n_rows()); self.analytes @@ -566,6 +572,17 @@ impl TargetColumns { self.build_source_ids()?; let ships_decoys = self.is_decoy.iter().any(|&d| d); self.caps.decoys = decoys.strategy(ships_decoys); + if matches!(self.caps.decoys, DecoyStrategy::MassShift { .. }) + && self + .frag_labels + .iter() + .any(|label| !label.supports_decoy_shift()) + { + tracing::info!( + "Mass-shift decoy generation disabled library-wide: missing fragment chemistry; using retained rows only" + ); + self.caps.decoys = DecoyStrategy::Stored; + } // Nothing to build when no groups were declared: a row is then its own // group, which `decoy_group_code` derives. Only the case that actually // loses information is worth a word. diff --git a/rust/timsquery/src/serde/library_file.rs b/rust/timsquery/src/serde/library_file.rs index 1e52d0c2..1c456b35 100644 --- a/rust/timsquery/src/serde/library_file.rs +++ b/rust/timsquery/src/serde/library_file.rs @@ -178,7 +178,7 @@ impl ElutionGroupCollection { /// `geom.frag_labels`/`geom.frag_mzs` (same length). The columnar store itself /// stores query geometry; readers that supply reference intensities populate /// the sidecar for scoring. Geometry-only `Mzpaf` inputs leave it `None`; -/// the `Str` variant has no intensity sidecar. Extraction ignores intensities. +/// both variants can carry the sidecar. Target-list JSON supplies none. Extraction ignores intensities. pub enum TargetTable { Mzpaf { geom: TargetColumns, @@ -186,6 +186,7 @@ pub enum TargetTable { }, Str { geom: TargetColumns>, + frag_intens: Option>, }, } @@ -339,7 +340,10 @@ impl TargetTable { }); } let geom = geom.seal(DecoyPolicy::Never)?; - Ok(TargetTable::Str { geom }) + Ok(TargetTable::Str { + geom, + frag_intens: None, + }) } } } @@ -578,7 +582,7 @@ mod tests { } match super::read_targets(path).expect("fixture loads") { super::TargetTable::Mzpaf { geom, .. } => ids(&geom), - super::TargetTable::Str { geom } => ids(&geom), + super::TargetTable::Str { geom, .. } => ids(&geom), } } diff --git a/rust/timsquery/src/traits/fragment_label.rs b/rust/timsquery/src/traits/fragment_label.rs index fe79ba1a..07cbf99e 100644 --- a/rust/timsquery/src/traits/fragment_label.rs +++ b/rust/timsquery/src/traits/fragment_label.rs @@ -9,6 +9,12 @@ use std::sync::Arc; /// label-generic flyweight compiles. Labels without ion chemistry apply the /// identity shift. pub trait DecoyShift { + /// Whether this label supplies the facts required for a mass-shift decoy. + /// The arena requires this for every retained fragment before generating any. + fn supports_decoy_shift(&self) -> bool { + false + } + /// Decoy m/z shift. Identity for labels without ion chemistry. fn decoy_shift_mz(&self, mz: f64, shift: f64) -> f64; } @@ -31,6 +37,13 @@ pub trait FragmentLabel: KeyLike + DecoyShift { } impl DecoyShift for IonAnnot { + fn supports_decoy_shift(&self) -> bool { + !matches!( + self.series_ordinal(), + micromzpaf::IonSeriesOrdinal::unknown { .. } + ) && self.get_charge() > 0 + } + fn decoy_shift_mz(&self, mz: f64, shift: f64) -> f64 { // Verbatim `create_mass_shifted_decoy` rule: shift ordinal>2 (and // ordinal-less) fragments by `shift / charge`; leave ordinal<=2 alone. diff --git a/rust/timsquery_cli/src/commands.rs b/rust/timsquery_cli/src/commands.rs index e61ec471..86c6bfe3 100644 --- a/rust/timsquery_cli/src/commands.rs +++ b/rust/timsquery_cli/src/commands.rs @@ -94,7 +94,7 @@ pub fn main_query_index(args: QueryIndexArgs) -> Result<(), CliError> { &put_path, batch_size, ), - TargetTable::Str { geom } => stream_process_batches( + TargetTable::Str { geom, .. } => stream_process_batches( &geom, aggregator_use, &index, @@ -493,7 +493,7 @@ mod tests { let arena = read_query_elution_groups(tmp_file.path()).unwrap(); let geom = match arena { - TargetTable::Str { geom } => geom, + TargetTable::Str { geom, .. } => geom, TargetTable::Mzpaf { .. } => { panic!("string labels must land in the Str arena, not Mzpaf") } diff --git a/rust/timsseek/src/data_sources/reference_library.rs b/rust/timsseek/src/data_sources/reference_library.rs index d03110ec..9cb93b08 100644 --- a/rust/timsseek/src/data_sources/reference_library.rs +++ b/rust/timsseek/src/data_sources/reference_library.rs @@ -33,6 +33,8 @@ pub struct ReferenceLibrary { /// Parallel to `geom.frag_labels` / `geom.frag_mzs`; same `frag_off` ranges. frag_intens: Vec, plan: plan::ScoringPlan, + /// Original opaque labels; aligned with packed extraction keys when present. + opaque_labels: Option>>, } pub trait ExpectedIntensity { @@ -54,6 +56,12 @@ impl ReferenceLibrary { &self.geom } + /// Source opaque labels in fragment order, when this library used string keys. + /// Packed geometry keys are internal extraction handles, not chemical annotations. + pub fn opaque_fragment_labels(&self) -> Option<&[std::sync::Arc]> { + self.opaque_labels.as_deref() + } + /// Reference intensities in the geometry's fragment order. pub fn fragment_intensities(&self) -> &[f32] { &self.frag_intens @@ -107,19 +115,11 @@ impl ReferenceLibrary { }) } - /// Narrow a label-generic [`TargetTable`] (timsquery's one library funnel) - /// into the ion-annotated `ReferenceLibrary` timsseek scores against. - /// - /// Scoring requires an `Mzpaf` arena with reference fragment intensities. - /// Rejects string labels, missing/misaligned intensities, repeated fragment - /// keys within a row, unsupported isotope counts, and nonempty libraries - /// with no annotated fragments. These are scoring restrictions; - /// query extraction also accepts geometry with opaque fragment labels. - /// Readers that supply reference intensities preserve them in the sidecar; - /// this includes the DIA-NN TSV/Parquet adapters. - /// - /// Resolves the operation plan from the validated arena. The public load - /// boundary additionally reports decoys, operation coverage, and isotope fallbacks. + /// Validate retained fragments and their aligned reference intensities, then + /// resolve library-wide operations. String keys are opaque: preserve their + /// source labels while assigning packed unknown keys for extraction. + /// Missing/misaligned intensities, duplicate keys, unsupported isotope counts + /// and nonempty libraries without retained fragments are rejected. pub(crate) fn from_arena(arena: TargetTable) -> Result { match arena { TargetTable::Mzpaf { geom, frag_intens } => { @@ -169,9 +169,7 @@ impl ReferenceLibrary { if geom.n_rows() > 0 && geom.n_fragments() == 0 { return Err(TargetReadingError::UnsupportedFormat { message: format!( - "library has {} entries and not one annotated fragment; timsseek \ - scores annotated peptide fragments, so there is nothing here to \ - score against", + "library has {} entries and no retained fragments; there is nothing to score against", geom.n_rows() ), }); @@ -182,12 +180,57 @@ impl ReferenceLibrary { geom, frag_intens, plan, + opaque_labels: None, }) } - TargetTable::Str { .. } => Err(TargetReadingError::UnsupportedFormat { - message: "timsseek requires ion-annotated fragments (mzpaf); got string labels" - .to_string(), - }), + TargetTable::Str { geom, frag_intens } => { + use timsquery::models::{ + Row, + TargetColumnsBuilder, + }; + let mut builder = + TargetColumnsBuilder::with_capabilities(geom.capabilities().clone()); + let mut opaque_labels = Vec::with_capacity(geom.n_fragments()); + for row in geom.rows() { + let mut seen = HashSet::new(); + let mut counter = timsquery::ion::UnknownIonCounter::default(); + let mut fragments = Vec::with_capacity(geom.frag_labels(row).len()); + for (label, &mz) in geom.frag_labels(row).iter().zip(geom.frag_mzs(row)) { + if !seen.insert(label) { + return Err(TargetReadingError::InvalidLibrary { + message: format!( + "source {:?}: duplicate fragment key {label:?}", + geom.output_id(row) + ), + }); + } + // Packed unknown keys only identify peaks. Their placeholder + // charge is never evidence for a chemistry-dependent operation. + let key = counter.next_unknown(1).map_err(|e| TargetReadingError::UnsupportedFormat { message: format!("source {:?}: opaque scoring currently supports at most 255 peaks per entry: {e}", geom.output_id(row)) })?; + fragments.push((key, mz)); + opaque_labels.push(label.clone()); + } + let analyte = geom.analyte(row).to_owned(); + builder.push_row(Row { + precursor_mz: geom.precursor_mz(row), + charge: geom.charge(row), + rt_seconds: geom.rt_seconds(row), + mobility: geom.mobility(row), + frags: &fragments, + analyte: analyte.as_input(), + entry_name: geom.entry_name(row), + is_decoy: geom.is_decoy(row), + id: Some(geom.output_id(row).to_owned_id()), + decoy_group: Some(geom.decoy_group(row).to_owned_id()), + }); + } + let geom = builder + .seal(crate::models::DecoyPolicy::Never) + .map_err(timsquery::serde::TargetReadingError::from)?; + let mut library = Self::from_arena(TargetTable::Mzpaf { geom, frag_intens })?; + library.opaque_labels = Some(opaque_labels); + Ok(library) + } } } } @@ -209,6 +252,10 @@ impl<'a> RefQuery<'a> { } } + pub fn scoring_plan(&self) -> &plan::ScoringPlan { + &self.lib.plan + } + pub fn geom(&self) -> &Query<&'a TargetColumns, IonAnnot> { &self.geom } @@ -444,8 +491,7 @@ impl ReferenceLibrary { /// What will actually be scored on the decoy side of the FDR estimate. /// /// Read off the counts rather than off `caps.decoys`: what matters at load - /// time is whether anything is there, and the one case worth a warning -- - /// nothing derived and nothing shipped -- is not a strategy. + /// time is whether any decoys will actually be scored. fn report_decoys(&self) { let n_rows = self.geom.n_rows(); let n_stored_decoys = self.geom.n_stored_decoys(); @@ -459,10 +505,9 @@ impl ReferenceLibrary { scored entries", ); } else if n_stored_decoys == 0 { - tracing::warn!( - "Library ships no decoys and none will be derived; scoring {n_rows} stored \ - rows as-is. FDR would be estimated with nothing to estimate it from. Use \ - --decoy-strategy if-missing to derive mass-shift decoys.", + tracing::info!( + "Library ships no decoys and none will be derived; search will write raw \ + scores for {n_rows} stored rows without FDR estimates", ); } else { tracing::info!( @@ -652,10 +697,51 @@ mod tests { let sgeom = sgeom .seal(crate::models::DecoyPolicy::Never) .expect("an empty arena seals"); - let s = TargetTable::Str { geom: sgeom }; + let s = TargetTable::Str { + geom: sgeom, + frag_intens: None, + }; assert!(ReferenceLibrary::try_from(s).is_err()); } + #[test] + fn opaque_labels_keep_identity_geometry_intensity_and_disable_generation() { + use std::sync::Arc; + use timsquery::models::OwnedSourceId; + let mut builder = TargetColumnsBuilder::>::with_capabilities( + TargetCapabilities::default_unlabeled(), + ); + let labels: [Arc; 2] = [Arc::from("arbitrary product"), Arc::from("y3")]; + builder.push_row(Row { + id: Some(OwnedSourceId::Text("original-id".into())), + entry_name: Some("independent name"), + precursor_mz: 500.0, + charge: 2, + frags: &[(labels[0].clone(), 200.0), (labels[1].clone(), 300.0)], + ..Default::default() + }); + let geom = builder.seal(crate::models::DecoyPolicy::IfMissing).unwrap(); + assert_eq!(geom.variants_per_row(), 1); + let lib = ReferenceLibrary::try_from(TargetTable::Str { + geom, + frag_intens: Some(vec![1.0, 0.5]), + }) + .unwrap(); + assert_eq!(lib.opaque_fragment_labels().unwrap(), &labels); + assert!(!lib.scoring_plan().fragment_isotopes().enabled); + assert_eq!(lib.scoring_plan().width(), 0); + let q = lib.iter().next().unwrap(); + assert_eq!(q.output_id().to_string(), "original-id"); + let row = q.geom().row(); + assert_eq!(lib.geometry().entry_name(row), Some("independent name")); + assert_eq!(lib.geometry().frag_mzs(row), &[200.0, 300.0]); + assert_eq!(lib.fragment_intensities(), &[1.0, 0.5]); + assert!(lib.geometry().frag_labels(row).iter().all(|label| matches!( + label.series_ordinal(), + timsquery::ion::IonSeriesOrdinal::unknown { .. } + ))); + } + #[test] fn expected_fragments_pair_labels_with_intensities() { let lib = tiny_ref_lib(); @@ -1334,41 +1420,14 @@ mod load_tests { assert_ne!(groups[0], groups[3], "separate rows do not"); } - /// `Force` over a library that ships its own decoys, the one state a - /// hand-built arena cannot reach. - /// - /// The drop is `DecoyPolicy::accepts`, applied by the reader as it pushes - /// rows, so `seal` never sees the decoys at all; a fixture that dropped them - /// itself would be asserting on the test's own arithmetic. Read through - /// `from_file` for that reason. timsquery's - /// `skipping_shipped_decoys_leaves_only_targets` covers the arena side of - /// the same load; this is what scoring then reads off it. #[test] - fn forcing_mass_shift_decoys_replaces_the_ones_a_library_shipped() { + fn placeholder_fragment_charges_disable_mass_shift_decoys() { let path = timsquery_fixture("mzspeclib_files/target_decoy_attribute_set.mzspeclib.txt"); - let lib = + let library = ReferenceLibrary::from_file(&path, deciding_decoys(crate::models::DecoyPolicy::Force)) - .expect("the fixture loads"); - - assert_eq!( - lib.geom.n_rows(), - 5, - "the five targets, without their decoys" - ); - assert_eq!(lib.geom.n_stored_decoys(), 0, "the shipped decoys are gone"); - assert_eq!(lib.geom.variants_per_row(), 3); - - let marks = marks_of(&lib); - assert_eq!(marks.len(), 15, "five rows, three variants each"); - assert!( - marks.chunks(3).all(|row| row - == [ - DecoyMarking::Target, - DecoyMarking::MassShiftedDecoy, - DecoyMarking::MassShiftedDecoy - ]), - "every row is a target with two derived decoys, got {marks:?}" - ); + .unwrap(); + assert_eq!(library.geom.variants_per_row(), 1); + assert_eq!(library.geom.n_stored_decoys(), 0); } /// The exclusion `decoy_marking` reads as "variant 0 and not a target means @@ -1512,32 +1571,24 @@ mod load_tests { panic!("expected an unsupported-format refusal, got {err:?}"); }; assert!( - message.contains("not one annotated fragment"), + message.contains("no retained fragments"), "the refusal has to say what is missing: {message}" ); } - /// The same library under the default policy, which keeps an unannotated - /// peak at the m/z the file measured: the arena has fragments, so the guard - /// above has nothing to refuse and the library is searchable. - /// - /// Five peaks, one row, and no sequence -- a small molecule has none, so the - /// plan disables sequence operations while preserving spectrum features. + /// Opaque peaks can serve common scoring/viewing without invented chemistry. + /// A library without suitable decoys cannot use the supervised CLI route. #[test] - fn the_default_policy_makes_that_same_library_searchable() { - let lib = ReferenceLibrary::from_file( - &timsquery_fixture("mzspeclib_files/small_molecule.mzspeclib.txt"), - crate::models::LoadPolicy::default(), - ) - .expect("a library whose peaks were kept has something to score against"); - + fn sequence_free_library_preserves_peaks_and_disables_chemistry_operations() { + let path = timsquery_fixture("mzspeclib_files/small_molecule.mzspeclib.txt"); + let lib = ReferenceLibrary::from_file(&path, Default::default()).unwrap(); + assert_eq!(lib.geom.variants_per_row(), 1); assert_eq!(lib.geom.n_rows(), 1); assert_eq!(lib.geom.n_fragments(), 5); - assert_eq!(lib.frag_intens.len(), 5, "the sidecar stays parallel"); - assert!( - !lib.all_sequence_counts_enabled(), - "a small molecule has no sequence" - ); + assert_eq!(lib.frag_intens.len(), 5); + assert!(!lib.all_sequence_counts_enabled()); + assert!(!lib.scoring_plan().fragment_isotopes().enabled); + assert_eq!(lib.scoring_plan().fragment_isotopes().usable_fragments, 0); } /// Geometry-only arenas can serve extraction, but scoring requires the diff --git a/rust/timsseek/src/ml/qvalues.rs b/rust/timsseek/src/ml/qvalues.rs index 14c15797..a4facb67 100644 --- a/rust/timsseek/src/ml/qvalues.rs +++ b/rust/timsseek/src/ml/qvalues.rs @@ -419,7 +419,7 @@ pub fn rescore(mut data: Vec, library: &ReferenceLibrary) -> /// Selected via the `rescore_model` config field / `--rescore-model` CLI flag /// ([`crate::ml::RescoreModel::Lda`]). /// See `ml::lda` for the fit details. -pub fn rescore_lda(mut data: Vec, _library: &ReferenceLibrary) -> RescoreResult { +pub fn rescore_lda(mut data: Vec, library: &ReferenceLibrary) -> RescoreResult { // Canonical sort + seeded shuffle -- the same helper, key and seed as every // other rescorer. canonicalize_and_shuffle(&mut data); @@ -429,17 +429,15 @@ pub fn rescore_lda(mut data: Vec, _library: &ReferenceLibrary // emit time by the grammar, so there is no data-dependent normalization step // here -- the only remaining data-dependent op is LDA's own standardization. // Built after the shuffle, per `canonicalize_and_shuffle`. - let names: Vec> = linear_feature_name_set(); + let names: Vec> = linear_feature_name_set(library); let nrows = data.len(); - let ncols = LINEAR_NCOLS; + let ncols = linear_ncols(library); debug_assert_eq!(names.len(), ncols); - let dataset = StreamingDataset::new( - &data, - names.clone(), - N_RESCORE_FOLDS as usize, - &write_competed_linear_row, - ); + let write_row = |candidate: &CompetedCandidate, out: &mut [f64]| { + write_competed_linear_row(library, candidate, out) + }; + let dataset = StreamingDataset::new(&data, names.clone(), N_RESCORE_FOLDS as usize, &write_row); let cf = crossfit_lda(&dataset)?; debug!("LDA cross-fit scored {nrows} candidates across {N_RESCORE_FOLDS} folds"); @@ -459,12 +457,16 @@ pub fn rescore_lda(mut data: Vec, _library: &ReferenceLibrary /// /// The fold count and assignment match the GBM's /// [`CrossValidatedScorer`] partition. -fn hybrid_linear_dataset(data: &[CompetedCandidate]) -> StreamingDataset<'_, CompetedCandidate> { +fn hybrid_linear_dataset<'a>( + library: &ReferenceLibrary, + data: &'a [CompetedCandidate], + write_row: &'a (dyn Fn(&CompetedCandidate, &mut [f64]) + Sync), +) -> StreamingDataset<'a, CompetedCandidate> { StreamingDataset::new( data, - linear_feature_name_set(), + linear_feature_name_set(library), N_RESCORE_FOLDS as usize, - &write_competed_linear_row, + write_row, ) } @@ -536,7 +538,10 @@ pub fn rescore_hybrid( // `canonicalize_and_shuffle`), so row `i` is the same candidate in the // streamed linear source, nonlinear frame, lda_score, responses, and moved data. let responses: Vec = data.iter().map(|c| c.get_y()).collect(); - let lin_dataset = hybrid_linear_dataset(&data); + let write_row = |candidate: &CompetedCandidate, out: &mut [f64]| { + write_competed_linear_row(library, candidate, out) + }; + let lin_dataset = hybrid_linear_dataset(library, &data, &write_row); let lda_score = crossfit::<_, LdaModel>(&lin_dataset, &LdaConfig::default(), "LDA")?.scores; @@ -635,8 +640,11 @@ const BASE_NONLINEAR_NCOLS: usize = fn nonlinear_ncols(library: &ReferenceLibrary) -> usize { BASE_NONLINEAR_NCOLS + library.scoring_plan().width() } +fn linear_ncols(library: &ReferenceLibrary) -> usize { + library.scoring_plan().linear_indices().len() + ResultMeta::LINEAR_LEN + Derived::LINEAR_LEN +} fn all_ncols(library: &ReferenceLibrary) -> usize { - LINEAR_NCOLS + nonlinear_ncols(library) + linear_ncols(library) + nonlinear_ncols(library) } const _: () = assert!(LINEAR_NCOLS > 0 && BASE_NONLINEAR_NCOLS > 0); @@ -686,12 +694,21 @@ impl ValueSink for SliceSink<'_> { /// Project one row's linear values into a streaming or retained sink. fn project_linear_row( + library: &ReferenceLibrary, scoring: &ScoringFields, meta: &ResultMeta, derived: &Derived, out: &mut impl ValueSink, ) { - out.push(&scoring.linear_feature_array()); + let values = scoring.linear_feature_array(); + let indices = library.scoring_plan().linear_indices(); + if indices.len() == values.len() { + out.push(&values); + } else { + for &i in indices { + out.push(&values[i..i + 1]); + } + } out.push(&meta.linear_feature_array()); out.push(&derived.linear_feature_array()); } @@ -726,29 +743,39 @@ fn write_competed_all_row( let meta = candidate.result_meta(); let derived = Derived::compute(scoring); let mut sink = SliceSink::new(out); - project_linear_row(scoring, &meta, &derived, &mut sink); + project_linear_row(library, scoring, &meta, &derived, &mut sink); project_nonlinear_row(library, scoring, &meta, &derived, &mut sink); sink.finish(); } /// Matrix-free LINEAR-lane projection for the standalone and hybrid LDA paths. -fn write_competed_linear_row(candidate: &CompetedCandidate, out: &mut [f64]) { - assert_eq!(out.len(), LINEAR_NCOLS); +fn write_competed_linear_row( + library: &ReferenceLibrary, + candidate: &CompetedCandidate, + out: &mut [f64], +) { + assert_eq!(out.len(), linear_ncols(library)); let scoring = &candidate.scoring; let meta = candidate.result_meta(); let derived = Derived::compute(scoring); let mut sink = SliceSink::new(out); - project_linear_row(scoring, &meta, &derived, &mut sink); + project_linear_row(library, scoring, &meta, &derived, &mut sink); sink.finish(); } /// The LINEAR-lane matrix for `data` in its current order. #[cfg(test)] -fn build_linear_matrix(data: &[CompetedCandidate]) -> Vec { - let mut out = Vec::with_capacity(data.len() * LINEAR_NCOLS); +fn build_linear_matrix(library: &ReferenceLibrary, data: &[CompetedCandidate]) -> Vec { + let mut out = Vec::with_capacity(data.len() * linear_ncols(library)); for c in data { let meta = c.result_meta(); - project_linear_row(&c.scoring, &meta, &Derived::compute(&c.scoring), &mut out); + project_linear_row( + library, + &c.scoring, + &meta, + &Derived::compute(&c.scoring), + &mut out, + ); } out } @@ -781,7 +808,7 @@ fn build_all_matrix<'a>( let mut out = Vec::with_capacity(rows.len() * all_ncols(library)); for (s, meta) in rows { let derived = Derived::compute(s); - project_linear_row(s, &meta, &derived, &mut out); + project_linear_row(library, s, &meta, &derived, &mut out); project_nonlinear_row(library, s, &meta, &derived, &mut out); } out @@ -810,12 +837,20 @@ pub fn feature_frame( } /// LINEAR-lane feature names (LDA), in `project_linear_row`'s order. -pub fn linear_feature_name_set() -> Vec> { +pub fn linear_feature_name_set(library: &ReferenceLibrary) -> Vec> { let mut n = NameSink::new(); ::linear_feature_names(&mut n); + let all = n.into_names(); + let mut n = NameSink::new(); ::linear_feature_names(&mut n); ::linear_feature_names(&mut n); - n.into_names() + library + .scoring_plan() + .linear_indices() + .iter() + .map(|&i| all[i].clone()) + .chain(n.into_names()) + .collect() } /// NONLINEAR-lane feature names, in `project_nonlinear_row`'s order. The @@ -833,7 +868,7 @@ pub fn nonlinear_feature_name_set(library: &ReferenceLibrary) -> Vec> { /// The ALL-lane feature names (GBM) = linear ++ nonlinear, matching /// `build_all_matrix`'s column order. pub fn all_feature_name_set(library: &ReferenceLibrary) -> Vec> { - let mut v = linear_feature_name_set(); + let mut v = linear_feature_name_set(library); v.extend(nonlinear_feature_name_set(library)); v } @@ -988,6 +1023,12 @@ mod feature_tests { )) } fn library_with_analyte(analyte: timsquery::chemistry::analyte::Analyte) -> ReferenceLibrary { + library_with_labels(analyte, &["y1"]) + } + fn library_with_labels( + analyte: timsquery::chemistry::analyte::Analyte, + labels: &[&str], + ) -> ReferenceLibrary { use timsquery::models::{ Row, TargetColumnsBuilder, @@ -996,7 +1037,10 @@ mod feature_tests { let mut builder = TargetColumnsBuilder::with_capabilities( timsquery::models::TargetCapabilities::default_diann(), ); - let frags = [(timsquery::ion::IonAnnot::try_from("y1").unwrap(), 300.0)]; + let frags: Vec<_> = labels + .iter() + .map(|label| (timsquery::ion::IonAnnot::try_from(*label).unwrap(), 300.0)) + .collect(); for _ in 0..1024 { builder.push_row(Row { analyte: analyte.as_input(), @@ -1007,8 +1051,10 @@ mod feature_tests { }); } ReferenceLibrary::from_sealed_arena(TargetTable::Mzpaf { - geom: builder.seal(Default::default()).unwrap(), - frag_intens: Some(vec![1.0; 1024]), + geom: builder + .seal(timsquery::models::capabilities::DecoyPolicy::Never) + .unwrap(), + frag_intens: Some(vec![1.0; 1024 * labels.len()]), }) .unwrap() } @@ -1061,6 +1107,40 @@ mod feature_tests { } } + #[test] + fn one_unknown_fragment_disables_isotope_features_in_every_projection() { + let full = library(); + let mixed = library_with_labels( + timsquery::chemistry::analyte::Analyte::from_sequence("PEPTIDEK"), + &["y1", "?1"], + ); + assert!(!mixed.scoring_plan().fragment_isotopes().enabled); + assert_eq!( + mixed.scoring_plan().fragment_isotopes().usable_fragments, + 1024 + ); + let data = vec![sample_competed_candidate()]; + let full_names = all_feature_name_set(full); + let full_values = build_all_matrix(full, competed_rows(&data)); + let names = all_feature_name_set(&mixed); + let values = build_all_matrix(&mixed, competed_rows(&data)); + assert_eq!(full_names.len() - names.len(), 2); + assert_eq!(values.len(), names.len()); + for (name, value) in names.iter().zip(&values) { + assert!(!name.starts_with("ms2_isotope_lazyscore")); + let i = full_names.iter().position(|n| n == name).unwrap(); + assert_eq!(value.to_bits(), full_values[i].to_bits()); + } + let mut streamed = vec![0.0; values.len()]; + write_competed_all_row(&mixed, &data[0], &mut streamed); + assert_eq!(streamed, values); + let linear = build_linear_matrix(&mixed, &data); + let mut streamed = vec![0.0; linear.len()]; + write_competed_linear_row(&mixed, &data[0], &mut streamed); + assert_eq!(streamed, linear); + assert_eq!(linear.len(), linear_feature_name_set(&mixed).len()); + } + // --- Lane walks (the live ML path) --- /// LANE WIDTH PARITY: a lane matrix's row width MUST equal that lane's @@ -1075,14 +1155,14 @@ mod feature_tests { fn lane_matrix_widths_match_name_sets() { for context in [library_with_sequence("PEPTIDEK"), library_with_sequence("")] { let data = vec![sample_competed_candidate()]; - assert_eq!(linear_feature_name_set().len(), LINEAR_NCOLS); + assert_eq!(linear_feature_name_set(library()).len(), LINEAR_NCOLS); assert_eq!( nonlinear_feature_name_set(&context).len(), nonlinear_ncols(&context) ); assert_eq!(all_feature_name_set(&context).len(), all_ncols(&context)); - assert_eq!(build_linear_matrix(&data).len(), LINEAR_NCOLS); + assert_eq!(build_linear_matrix(library(), &data).len(), LINEAR_NCOLS); assert_eq!( build_nonlinear_matrix(&context, &data).len(), nonlinear_ncols(&context) @@ -1100,7 +1180,7 @@ mod feature_tests { #[test] fn all_matrix_is_linear_then_nonlinear_per_row() { let data = vec![sample_competed_candidate(), sample_competed_candidate()]; - let lin = build_linear_matrix(&data); + let lin = build_linear_matrix(library(), &data); let nl = build_nonlinear_matrix(library(), &data); let all = build_all_matrix(library(), competed_rows(&data)); @@ -1108,7 +1188,7 @@ mod feature_tests { for i in 0..data.len() { let row = &all[i * ALL_NCOLS..(i + 1) * ALL_NCOLS]; let mut streamed_linear = vec![0.0; LINEAR_NCOLS]; - write_competed_linear_row(&data[i], &mut streamed_linear); + write_competed_linear_row(library(), &data[i], &mut streamed_linear); assert_eq!( bits(&lin[i * LINEAR_NCOLS..(i + 1) * LINEAR_NCOLS]), bits(&streamed_linear), @@ -1358,7 +1438,7 @@ mod feature_tests { } let lane: std::collections::HashSet> = - linear_feature_name_set().into_iter().collect(); + linear_feature_name_set(library()).into_iter().collect(); let reported: std::collections::HashSet> = stats[0] .feature_importance .iter() @@ -1674,7 +1754,9 @@ mod feature_tests { fn crossfit_rejects_incomplete_prediction_vectors() { let mut data = synthetic_competed(12); canonicalize_and_shuffle(&mut data); - let dataset = hybrid_linear_dataset(&data); + let write_row = + |c: &CompetedCandidate, out: &mut [f64]| write_competed_linear_row(library(), c, out); + let dataset = hybrid_linear_dataset(library(), &data, &write_row); let error = match crossfit::<_, ShortPredict>(&dataset, &(), "short predictor") { Err(error) => error, @@ -1698,7 +1780,9 @@ mod feature_tests { canonicalize_and_shuffle(&mut data); let responses: Vec = data.iter().map(|c| c.get_y()).collect(); - let lin_dataset = hybrid_linear_dataset(&data); + let write_row = + |c: &CompetedCandidate, out: &mut [f64]| write_competed_linear_row(library(), c, out); + let lin_dataset = hybrid_linear_dataset(library(), &data, &write_row); let cf = crossfit::<_, FoldSpy>(&lin_dataset, &(), "spy").expect("the spy cannot fail"); let (precomputed, names) = diff --git a/rust/timsseek/src/scoring/blocks/lazy.rs b/rust/timsseek/src/scoring/blocks/lazy.rs index 0e6dcfdb..e8dcc8d8 100644 --- a/rust/timsseek/src/scoring/blocks/lazy.rs +++ b/rust/timsseek/src/scoring/blocks/lazy.rs @@ -30,6 +30,13 @@ pub struct ApexLazyScores { pub struct SecondaryLazyScores { #[feat(ln1p)] pub ms2_lazyscore: f32, + #[block] + pub isotopes: FragmentIsotopeScores, +} + +/// Requires a usable fragment charge and representable +1 isotope for every peak. +#[derive(Debug, Clone, Copy, Serialize, ScoreBlock)] +pub struct FragmentIsotopeScores { #[feat(ln1p)] pub ms2_isotope_lazyscore: f32, #[feat(raw)] @@ -40,8 +47,10 @@ impl From for SecondaryLazyScores { fn from(s: SecondaryLazyScoresRaw) -> Self { Self { ms2_lazyscore: s.lazyscore, - ms2_isotope_lazyscore: s.iso_lazyscore, - ms2_isotope_lazyscore_log_diff: s.isotope_lazyscore_log_diff, + isotopes: FragmentIsotopeScores { + ms2_isotope_lazyscore: s.iso_lazyscore, + ms2_isotope_lazyscore_log_diff: s.isotope_lazyscore_log_diff, + }, } } } diff --git a/rust/timsseek/src/scoring/parquet_writer.rs b/rust/timsseek/src/scoring/parquet_writer.rs index b4656287..7c8914db 100644 --- a/rust/timsseek/src/scoring/parquet_writer.rs +++ b/rust/timsseek/src/scoring/parquet_writer.rs @@ -184,6 +184,30 @@ pub fn build_record_batch( .map_err(|e| std::io::Error::new(std::io::ErrorKind::InvalidData, e)) } +fn output_batch( + results: &[FinalResult], + geom: &TargetColumns, + raw: bool, +) -> std::io::Result { + let batch = build_record_batch(results, geom)?; + if !raw { + return Ok(batch); + } + let mut schema = SchemaSink::new(); + super::blocks::result_meta::ResultMeta::column_schema(&mut schema); + let omitted = schema.into_fields(); + let keep: Vec<_> = batch + .schema() + .fields() + .iter() + .enumerate() + .filter_map(|(i, field)| { + (!omitted.iter().any(|omit| omit.name() == field.name())).then_some(i) + }) + .collect(); + batch.project(&keep).map_err(std::io::Error::other) +} + // --------------------------------------------------------------------------- // Buffered Parquet writer // --------------------------------------------------------------------------- @@ -195,6 +219,7 @@ pub struct ResultParquetWriter<'a> { buffer: Vec, row_group_size: usize, geom: &'a TargetColumns, + raw: bool, } impl<'a> ResultParquetWriter<'a> { @@ -202,6 +227,24 @@ impl<'a> ResultParquetWriter<'a> { path: impl AsRef, row_group_size: usize, library: &'a crate::data_sources::reference_library::ReferenceLibrary, + ) -> std::io::Result { + Self::with_mode(path, row_group_size, library, false) + } + + /// Common scores only: no competition, discriminant score or q-value columns. + pub fn raw( + path: impl AsRef, + row_group_size: usize, + library: &'a crate::ReferenceLibrary, + ) -> std::io::Result { + Self::with_mode(path, row_group_size, library, true) + } + + fn with_mode( + path: impl AsRef, + row_group_size: usize, + library: &'a crate::ReferenceLibrary, + raw: bool, ) -> std::io::Result { let geom = library.geometry(); let file = match File::create_new(path.as_ref()) { @@ -213,10 +256,14 @@ impl<'a> ResultParquetWriter<'a> { }; // Build schema from a zero-row batch - let empty_batch = build_record_batch(&[], geom)?; + let empty_batch = output_batch(&[], geom, raw)?; let schema = empty_batch.schema(); let kv = vec![ + KeyValue { + key: "result_mode".into(), + value: Some(if raw { "raw" } else { "rescored" }.into()), + }, KeyValue { key: "scoring_plan".into(), value: Some( @@ -245,6 +292,18 @@ impl<'a> ResultParquetWriter<'a> { buffer: Vec::with_capacity(row_group_size), row_group_size, geom, + raw, + }) + } + + pub fn add_raw(&mut self, result: super::results::ScoredCandidate) -> std::io::Result<()> { + assert!(self.raw, "raw scores require the raw output schema"); + self.add(FinalResult { + scoring: result.scoring, + delta_group_ln1p_diff: f32::NAN, + delta_group_ln1p_ratio: f32::NAN, + discriminant_score: f32::NAN, + qvalue: f32::NAN, }) } @@ -261,7 +320,7 @@ impl<'a> ResultParquetWriter<'a> { return Ok(()); } debug!("Flushing {} results to parquet", self.buffer.len()); - let batch = build_record_batch(&self.buffer, self.geom)?; + let batch = output_batch(&self.buffer, self.geom, self.raw)?; self.writer.write(&batch).map_err(std::io::Error::other)?; self.buffer.clear(); Ok(()) @@ -290,6 +349,26 @@ mod tests { /// A sealed arena with one row per `(sequence, id)`, for the writer to /// resolve ids against. `None` for an id leaves the row to be minted. + #[test] + fn raw_output_omits_fdr_and_competition_columns_including_empty_files() { + let geom = one_row_arena(); + let row = sample_in(&geom); + for rows in [&[][..], std::slice::from_ref(&row)] { + let batch = output_batch(rows, &geom, true).unwrap(); + assert_eq!(batch.num_rows(), rows.len()); + for name in [ + "qvalue", + "discriminant_score", + "delta_group_ln1p_diff", + "delta_group_ln1p_ratio", + ] { + assert!(batch.schema().field_with_name(name).is_err()); + } + assert!(batch.schema().field_with_name("main_score").is_ok()); + assert!(batch.schema().field_with_name("library_id").is_ok()); + } + } + fn arena_of(rows: &[(&str, Option<&str>)]) -> TargetColumns { let mut geom = TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); for (seq, id) in rows { diff --git a/rust/timsseek/src/scoring/pipeline.rs b/rust/timsseek/src/scoring/pipeline.rs index 67454465..46f6b85f 100644 --- a/rust/timsseek/src/scoring/pipeline.rs +++ b/rust/timsseek/src/scoring/pipeline.rs @@ -482,14 +482,16 @@ pub struct SecondaryLazyScoresRaw { /// Compute lazyscores from inner and isotope collectors. fn compute_secondary_lazyscores( inner: &SpectralCollector, - isotope: &SpectralCollector, + isotope: Option<&SpectralCollector>, ) -> SecondaryLazyScoresRaw { let lazyscore = single_lazyscore( inner .iter_fragments() .map(|((_k, _mz), v)| v.weight() as f32), ); - let iso_lazyscore = single_lazyscore(isotope.iter_fragments().map(|((_k, _mz), v)| *v)); + let iso_lazyscore = isotope.map_or(f32::NAN, |isotope| { + single_lazyscore(isotope.iter_fragments().map(|((_k, _mz), v)| *v)) + }); let isotope_lazyscore_log_diff = lazyscore.ln_1p() - iso_lazyscore.ln_1p(); SecondaryLazyScoresRaw { lazyscore, @@ -590,6 +592,14 @@ impl Scorer { let inner = worker.inner_collector.as_mut().expect("init above"); inner.reset_with_overrides(query, Some(new_rt_seconds), Some(mobility as f32)); + self.index.add_query(inner, isotope_tol); + if !query.scoring_plan().fragment_isotopes().enabled { + // Workers may be reused with another library. Never consume a stale + // isotope collector when this library disables the operation. + worker.isotope_collector = None; + return; + } + // Isotope query holds `query` with +1 neutron offset applied // (buffer-override -- reuses Vec capacity after warm-up). let isotope_query = worker.isotope_query.get_or_insert_with(Target::empty_like); @@ -605,8 +615,6 @@ impl Scorer { isotope.reset_with_overrides(isotope_query, Some(new_rt_seconds), Some(mobility as f32)); // Both queries share the same isotope_tol per existing logic. - let inner = worker.inner_collector.as_mut().expect("init above"); - self.index.add_query(inner, isotope_tol); let isotope = worker.isotope_collector.as_mut().expect("init above"); self.index.add_query(isotope, isotope_tol); } @@ -621,7 +629,7 @@ impl Scorer { nqueries: u8, apex: ApexBlocks, inner_collector: &SpectralCollector, - isotope_collector: &SpectralCollector, + isotope_collector: Option<&SpectralCollector>, ) -> Result { let offsets = MzMobilityOffsets::new(inner_collector, metadata.ref_mobility_ook0 as f64); let rel_inten = RelativeIntensityCollector::new(inner_collector); @@ -690,17 +698,17 @@ impl Scorer { }) } - /// Phase 3: Score a peptide using calibrated extraction window. - /// Expects narrow calibrated extraction (from CalibrationResult). + /// Score with calibrated extraction, or full-RT extraction when calibration + /// is absent. Both modes share apex scoring and the enabled secondary queries. #[cfg_attr( feature = "instrumentation", tracing::instrument(skip_all, level = "trace") )] - fn score_calibrated_extraction( + fn score_extraction( &self, query: &RefQuery<'_>, identity: CandidateIdentity, - calibration: &CalibrationResult, + calibration: Option<&CalibrationResult>, worker: &mut ScoringWorker, timings: &mut ScoreTimings, ) -> Result { @@ -708,9 +716,36 @@ impl Scorer { let metadata = timed!( timings.extraction, - tracing::span!(tracing::Level::TRACE, "score_calibrated::extraction").in_scope( - || self.build_calibrated_extraction_into(query, identity, calibration, worker) - ) + tracing::span!(tracing::Level::TRACE, "score_calibrated::extraction").in_scope(|| { + match calibration { + Some(calibration) => { + self.build_calibrated_extraction_into(query, identity, calibration, worker) + } + None => { + let tolerance = self.broad_tolerance.clone().with_rt_tolerance( + timsquery::models::tolerance::RtTolerance::Unrestricted, + ); + super::extraction::build_extraction_into( + &mut worker.extraction, + query, + None, + &self.index, + &tolerance, + Some(TOP_N_FRAGMENTS), + )?; + Ok(super::apex_finding::CandidateMetadata { + is_target: identity.is_target, + charge: query.precursor_charge(), + handles: identity.handles, + source_id: identity.source_id, + library_rt: query.rt_seconds(), + calibrated_rt_seconds: f32::NAN, + ref_mobility_ook0: query.mobility_ook0(), + ref_precursor_mz: query.mono_precursor_mz(), + }) + } + } + }) )?; let scoring_ctx = worker @@ -731,8 +766,21 @@ impl Scorer { )?; timed!(timings.spectral_query, { - let spectral_tol = calibration.get_spectral_tolerance(); - let isotope_tol = calibration.get_isotope_tolerance(); + let spectral_tol = calibration.map_or_else( + || { + self.broad_tolerance.clone().with_rt_tolerance( + timsquery::models::tolerance::RtTolerance::Minutes(( + 0.5 / 60.0, + 0.5 / 60.0, + )), + ) + }, + CalibrationResult::get_spectral_tolerance, + ); + let isotope_tol = calibration.map_or_else( + || spectral_tol.clone(), + CalibrationResult::get_isotope_tolerance, + ); tracing::span!(tracing::Level::TRACE, "score_calibrated::secondary_query").in_scope( || { self.execute_secondary_query( @@ -746,7 +794,7 @@ impl Scorer { ) }); let inner_collector = worker.inner_collector.as_ref().expect("set by secondary"); - let isotope_collector = worker.isotope_collector.as_ref().expect("set by secondary"); + let isotope_collector = worker.isotope_collector.as_ref(); let scoring_ctx = worker .extraction @@ -781,6 +829,25 @@ impl Scorer { lib: &ReferenceLibrary, flats: &[FlatIdx], calibration: &CalibrationResult, + ) -> (Vec, ScoreTimings, SkipCounts) { + self.score_batch(lib, flats, Some(calibration)) + } + + /// Score spectra over the acquisition's full RT range, without calibration + /// or assumptions about the library's RT units. No FDR is assigned here. + pub fn score_raw_batch( + &self, + lib: &ReferenceLibrary, + flats: &[FlatIdx], + ) -> (Vec, ScoreTimings, SkipCounts) { + self.score_batch(lib, flats, None) + } + + fn score_batch( + &self, + lib: &ReferenceLibrary, + flats: &[FlatIdx], + calibration: Option<&CalibrationResult>, ) -> (Vec, ScoreTimings, SkipCounts) { let get_item = |f| lib.item_at(f); let num_cycles = self.num_cycles(); @@ -833,13 +900,8 @@ impl Scorer { source_id: q.output_id().to_owned_id(), is_target: q.is_target(), }; - let result = self.score_calibrated_extraction( - &q, - identity, - calibration, - &mut worker, - &mut t, - ); + let result = + self.score_extraction(&q, identity, calibration, &mut worker, &mut t); (worker, acc.fold((result, t))) }, |(wa, a), (_wb, b)| (wa, a.reduce(b)), diff --git a/rust/timsseek/src/scoring/plan.rs b/rust/timsseek/src/scoring/plan.rs index d1b7ddd9..e7e07236 100644 --- a/rust/timsseek/src/scoring/plan.rs +++ b/rust/timsseek/src/scoring/plan.rs @@ -11,6 +11,7 @@ //! indicators. Existing observation-derived blocks retain their fixed widths; //! their all-NaN checks remain data-dependent guards. Sequence operations have //! no Parquet score columns; their decisions and coverage are file metadata. +use super::blocks::lazy::FragmentIsotopeScores; use super::blocks::sequence_counts::{ ModificationCounts, ResidueCounts, @@ -19,6 +20,7 @@ use super::blocks::{ NameSink, ScoreBlock, }; +use super::results::ScoringFields; use crate::fragment_mass::isotope_plan::IsotopePlan; use serde::Serialize; use std::sync::Arc; @@ -27,7 +29,10 @@ use timsquery::chemistry::analyte::{ PeptideRef, PropertyRef, }; -use timsquery::ion::IonAnnot; +use timsquery::ion::{ + IonAnnot, + IonSeriesOrdinal, +}; use timsquery::models::TargetColumns; #[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] @@ -121,10 +126,22 @@ impl std::fmt::Display for OperationDecision { } } +/// Resolved over every retained fragment of every stored target and decoy. +#[derive(Debug, Clone, Serialize)] +pub struct FragmentIsotopeDecision { + pub enabled: bool, + pub usable_fragments: usize, + pub total_fragments: usize, + pub columns: Vec>, +} + /// Owned by its reference library; callers cannot install a plan from another library. #[derive(Debug, Clone, Serialize)] pub struct ScoringPlan { isotopes: IsotopePlan, + fragment_isotopes: FragmentIsotopeDecision, + #[serde(skip)] + linear_indices: Vec, rows: usize, unmodified_rows: usize, operations: Vec, @@ -132,6 +149,34 @@ pub struct ScoringPlan { impl ScoringPlan { pub(crate) fn resolve(geom: &TargetColumns) -> Result { let isotopes = IsotopePlan::resolve(geom)?; + let total_fragments = geom.n_fragments(); + let usable_fragments = geom + .rows() + .flat_map(|row| geom.frag_labels(row)) + .filter(|label| { + !matches!(label.series_ordinal(), IonSeriesOrdinal::unknown { .. }) + && label.get_charge() > 0 + && label.try_with_offset_neutrons(1).is_ok() + }) + .count(); + let enabled = total_fragments > 0 && usable_fragments == total_fragments; + let mut isotope_names = NameSink::new(); + FragmentIsotopeScores::linear_feature_names(&mut isotope_names); + let isotope_names = isotope_names.into_names(); + let mut linear_names = NameSink::new(); + ScoringFields::linear_feature_names(&mut linear_names); + let linear_indices = linear_names + .into_names() + .iter() + .enumerate() + .filter_map(|(i, name)| (enabled || !isotope_names.contains(name)).then_some(i)) + .collect(); + let fragment_isotopes = FragmentIsotopeDecision { + enabled, + usable_fragments, + total_fragments, + columns: if enabled { isotope_names } else { Vec::new() }, + }; let mut unmodified_rows = 0; let mut residues = Coverage::default(); let mut modifications = Coverage::default(); @@ -177,6 +222,8 @@ impl ScoringPlan { .collect(); Ok(Self { isotopes, + fragment_isotopes, + linear_indices, rows, unmodified_rows, operations, @@ -186,11 +233,30 @@ impl ScoringPlan { /// Plan-level counts, shared by CLI and viewer reporting. pub fn summary(&self) -> String { format!( - "{} library entries; {} unmodified (known empty modification list); {}", - self.rows, self.unmodified_rows, self.isotopes + "{} library entries; {} unmodified (known empty modification list); {}; fragment isotope scoring: {} ({}/{} usable fragments)", + self.rows, + self.unmodified_rows, + self.isotopes, + if self.fragment_isotopes.enabled { + "enabled" + } else { + "disabled" + }, + self.fragment_isotopes.usable_fragments, + self.fragment_isotopes.total_fragments ) } + pub fn fragment_isotopes(&self) -> &FragmentIsotopeDecision { + &self.fragment_isotopes + } + + /// Positions in the derive-generated ScoringFields linear lane. Names and + /// values use the same selection; disabled scores never enter a model. + pub(crate) fn linear_indices(&self) -> &[usize] { + &self.linear_indices + } + pub fn isotopes(&self) -> &IsotopePlan { &self.isotopes } diff --git a/rust/timsseek/src/scoring/timings.rs b/rust/timsseek/src/scoring/timings.rs index eb7e44b4..8fe720e0 100644 --- a/rust/timsseek/src/scoring/timings.rs +++ b/rust/timsseek/src/scoring/timings.rs @@ -137,6 +137,8 @@ impl std::ops::AddAssign for PrescoreTimings { /// All timing fields are in milliseconds. #[derive(Debug, Default, Serialize)] pub struct PipelineReport { + /// No calibration, competition, rescoring or q-value filtering was applied. + pub raw_scores: bool, // Per-file: index loading (ms) pub load_index_ms: u64, diff --git a/rust/timsseek_cli/src/processing.rs b/rust/timsseek_cli/src/processing.rs index e258bd19..762b2675 100644 --- a/rust/timsseek_cli/src/processing.rs +++ b/rust/timsseek_cli/src/processing.rs @@ -105,6 +105,49 @@ fn write_feature_stats_sidecar( Ok(()) } +/// Decoy-free libraries use full-RT spectrum scoring and a schema without FDR fields. +fn execute_raw_pipeline( + library: &ReferenceLibrary, + pipeline: &Scorer, + options: &PipelineOptions<'_>, +) -> Result { + info!( + "No decoys: raw spectrum scoring over full acquisition RT; calibration, competition, rescoring and q-value filtering disabled" + ); + let path = std::path::Path::new(&options.output.uri).join(RESULTS_PARQUET); + let io_error = |source| TimsSeekError::Io { + source, + path: Some(path.clone()), + }; + let mut writer = + timsseek::scoring::parquet_writer::ResultParquetWriter::raw(&path, 20_000, library) + .map_err(io_error)?; + let mut report = PipelineReport { + raw_scores: true, + ..Default::default() + }; + let mut timings = ScoreTimings::default(); + for batch in library.chunks(options.chunk_size) { + let (rows, timing, skips) = pipeline.score_raw_batch(library, &batch); + timings += timing; + report.phase3_skips += skips; + report.total_scored += rows.len(); + for row in rows { + writer.add_raw(row).map_err(io_error)?; + } + } + writer.close().map_err(io_error)?; + report.phase3_extraction_thread_ms = timings.extraction.as_millis() as u64; + report.phase3_scoring_thread_ms = timings.scoring.as_millis() as u64; + report.phase3_spectral_query_thread_ms = timings.spectral_query.as_millis() as u64; + report.phase3_assembly_thread_ms = timings.assembly.as_millis() as u64; + println!( + "{} raw scored candidates; no FDR estimates", + report.total_scored + ); + Ok(report) +} + #[cfg_attr( feature = "instrumentation", tracing::instrument(skip_all, level = "trace") @@ -124,6 +167,13 @@ pub fn execute_pipeline( pipeline: &Scorer, options: &PipelineOptions<'_>, ) -> std::result::Result { + if !matches!( + speclib.geometry().capabilities().decoys, + timsquery::models::capabilities::DecoyStrategy::MassShift { .. } + ) && speclib.geometry().n_stored_decoys() == 0 + { + return execute_raw_pipeline(speclib, pipeline, options); + } let PipelineOptions { chunk_size, output: out_path, @@ -369,6 +419,7 @@ pub fn execute_pipeline( println!("{} targets at 1% FDR", targets_at_1pct_qval); Ok(PipelineReport { + raw_scores: false, load_index_ms: 0, // set by caller after return phase1_prescore_ms: phase1_ms, phase1_detail: phase1_timings, diff --git a/rust/timsseek_cli/tests/build_library.rs b/rust/timsseek_cli/tests/build_library.rs index ad0e60d4..35efe339 100644 --- a/rust/timsseek_cli/tests/build_library.rs +++ b/rust/timsseek_cli/tests/build_library.rs @@ -35,7 +35,7 @@ fn build_library_writes_a_readable_library_and_sidecar() { .expect("generated library reads back"); let rows = match table { timsquery::serde::TargetTable::Mzpaf { geom, .. } => geom.n_rows(), - timsquery::serde::TargetTable::Str { geom } => geom.n_rows(), + timsquery::serde::TargetTable::Str { geom, .. } => geom.n_rows(), }; assert!(rows > 0); From d517f862a6a25c115e6c56ebd1f5ab79de14af0f Mon Sep 17 00:00:00 2001 From: "J. Sebastian Paez" Date: Sat, 12 Sep 2026 21:27:39 -0700 Subject: [PATCH 2/2] Include resolved library fallbacks in run reports --- docs/development.md | 8 ++++ rust/timsquery/src/models/capabilities.rs | 48 ++++++++++++++++++- rust/timsquery/src/models/target_columns.rs | 27 ++++++++--- .../src/data_sources/reference_library.rs | 28 ++++++++++- rust/timsseek/src/scoring/plan.rs | 2 + rust/timsseek/src/scoring/timings.rs | 8 +++- rust/timsseek_cli/src/search.rs | 2 + 7 files changed, 113 insertions(+), 10 deletions(-) diff --git a/docs/development.md b/docs/development.md index 551aa762..f8441048 100644 --- a/docs/development.md +++ b/docs/development.md @@ -31,6 +31,14 @@ q-value columns (`result_mode=raw` Parquet metadata). It bypasses calibration, rescoring and q-value filtering. Supplied decoys do not by themselves validate a decoy strategy for a new analyte class. +`run_report.json` records the resolved `scoring_plan` once for the search library, +plus `calibration_scoring_plan` when a separate calibration library is supplied. +These are the same plans serialized in Parquet metadata: sequence-operation +coverage, precursor-isotope method/reasons, fragment-isotope availability, and +`decoys` (requested policy, resolved strategy, reason, stored-decoy count). +Per-file `pipeline.raw_scores` records whether raw scoring was used. + + Programmatic `TargetTable::Str` accepts an optional intensity sidecar. Scoring assigns unique packed unknown keys (currently at most 255 peaks per entry), preserving the original opaque labels separately; even a string such as `y3` diff --git a/rust/timsquery/src/models/capabilities.rs b/rust/timsquery/src/models/capabilities.rs index e8e09d81..539413cd 100644 --- a/rust/timsquery/src/models/capabilities.rs +++ b/rust/timsquery/src/models/capabilities.rs @@ -34,7 +34,8 @@ pub enum IsotopeStrategy { /// and generates none is [`Stored`](Self::Stored) with /// `n_stored_decoys() == 0`; that is a property of the rows, not of the /// strategy, and reporting it belongs to whoever counts rows. -#[derive(Debug, Clone, Copy, PartialEq)] +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +#[serde(tag = "method", rename_all = "snake_case")] pub enum DecoyStrategy { /// Score the stored rows as they are. Any decoys the file shipped are the /// decoys; nothing is derived. @@ -45,6 +46,51 @@ pub enum DecoyStrategy { MassShift { offset: f64 }, } +/// Why sealing selected its decoy strategy. Preserved for downstream reports. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] +#[serde(rename_all = "snake_case")] +pub enum DecoyResolutionReason { + PolicyNever, + SuppliedDecoys, + MissingFragmentChemistry, + FragmentChemistryAvailable, +} + +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct DecoyResolution { + pub requested: DecoyPolicy, + pub strategy: DecoyStrategy, + pub reason: DecoyResolutionReason, + pub stored_decoys: usize, +} + +impl DecoyResolution { + pub(crate) fn resolve( + requested: DecoyPolicy, + stored_decoys: usize, + generation_supported: bool, + ) -> Self { + let mut strategy = requested.strategy(stored_decoys > 0); + let reason = match (requested, strategy) { + (DecoyPolicy::Never, _) => DecoyResolutionReason::PolicyNever, + (_, DecoyStrategy::Stored) => DecoyResolutionReason::SuppliedDecoys, + (_, DecoyStrategy::MassShift { .. }) if generation_supported => { + DecoyResolutionReason::FragmentChemistryAvailable + } + _ => { + strategy = DecoyStrategy::Stored; + DecoyResolutionReason::MissingFragmentChemistry + } + }; + Self { + requested, + strategy, + reason, + stored_decoys, + } + } +} + /// What the caller wants done about decoys: the whole decoy decision, stated /// once, before the file is read. /// diff --git a/rust/timsquery/src/models/target_columns.rs b/rust/timsquery/src/models/target_columns.rs index 23ce6dda..d4dec02b 100644 --- a/rust/timsquery/src/models/target_columns.rs +++ b/rust/timsquery/src/models/target_columns.rs @@ -6,6 +6,8 @@ use crate::chemistry::analyte::{ }; use crate::models::capabilities::{ DecoyPolicy, + DecoyResolution, + DecoyResolutionReason, DecoyStrategy, MASS_SHIFT_VARIANTS, TargetCapabilities, @@ -222,6 +224,7 @@ pub enum TargetBuildError { #[derive(Debug, Clone)] pub struct TargetColumns { pub(crate) caps: TargetCapabilities, + decoy_resolution: Option, // per-target scalars, len = n_rows; addressed by `RowIdx`, never by a // caller-supplied id (those live in `source_ids`) pub(crate) precursor_mz: Vec, @@ -299,6 +302,7 @@ impl TargetColumns { fn empty(caps: TargetCapabilities) -> Self { Self { caps, + decoy_resolution: None, precursor_mz: Vec::new(), charge: Vec::new(), rt_seconds: Vec::new(), @@ -571,18 +575,20 @@ impl TargetColumns { self.build_decoy_groups()?; self.build_source_ids()?; let ships_decoys = self.is_decoy.iter().any(|&d| d); - self.caps.decoys = decoys.strategy(ships_decoys); - if matches!(self.caps.decoys, DecoyStrategy::MassShift { .. }) - && self - .frag_labels + let resolution = DecoyResolution::resolve( + decoys, + stored_decoys, + self.frag_labels .iter() - .any(|label| !label.supports_decoy_shift()) - { + .all(|label| label.supports_decoy_shift()), + ); + if resolution.reason == DecoyResolutionReason::MissingFragmentChemistry { tracing::info!( "Mass-shift decoy generation disabled library-wide: missing fragment chemistry; using retained rows only" ); - self.caps.decoys = DecoyStrategy::Stored; } + self.caps.decoys = resolution.strategy; + self.decoy_resolution = Some(resolution); // Nothing to build when no groups were declared: a row is then its own // group, which `decoy_group_code` derives. Only the case that actually // loses information is worth a word. @@ -602,6 +608,13 @@ impl TargetColumns { Ok(self) } + /// The decision made during sealing, including why generation was disabled. + pub fn decoy_resolution(&self) -> &DecoyResolution { + self.decoy_resolution + .as_ref() + .expect("sealed arenas have a decoy resolution") + } + pub fn capabilities(&self) -> &TargetCapabilities { &self.caps } diff --git a/rust/timsseek/src/data_sources/reference_library.rs b/rust/timsseek/src/data_sources/reference_library.rs index 9cb93b08..4335201e 100644 --- a/rust/timsseek/src/data_sources/reference_library.rs +++ b/rust/timsseek/src/data_sources/reference_library.rs @@ -225,7 +225,7 @@ impl ReferenceLibrary { }); } let geom = builder - .seal(crate::models::DecoyPolicy::Never) + .seal(geom.decoy_resolution().requested) .map_err(timsquery::serde::TargetReadingError::from)?; let mut library = Self::from_arena(TargetTable::Mzpaf { geom, frag_intens })?; library.opaque_labels = Some(opaque_labels); @@ -728,6 +728,32 @@ mod tests { }) .unwrap(); assert_eq!(lib.opaque_fragment_labels().unwrap(), &labels); + let calibration_library = tiny_ref_lib(); + let report = crate::scoring::RunReport { + scoring_plan: Some(lib.scoring_plan()), + calibration_scoring_plan: Some(calibration_library.scoring_plan()), + ..Default::default() + }; + let report = serde_json::to_value(report).unwrap(); + assert_eq!(report["scoring_plan"]["decoys"]["requested"], "if_missing"); + assert_eq!( + report["scoring_plan"]["decoys"]["strategy"]["method"], + "stored" + ); + assert_eq!( + report["scoring_plan"]["decoys"]["reason"], + "missing_fragment_chemistry" + ); + assert_eq!( + report["calibration_scoring_plan"]["decoys"]["strategy"]["method"], + "mass_shift" + ); + assert!( + report["scoring_plan"]["isotopes"] + .get("envelopes") + .is_none() + ); + assert!(report["scoring_plan"].get("linear_indices").is_none()); assert!(!lib.scoring_plan().fragment_isotopes().enabled); assert_eq!(lib.scoring_plan().width(), 0); let q = lib.iter().next().unwrap(); diff --git a/rust/timsseek/src/scoring/plan.rs b/rust/timsseek/src/scoring/plan.rs index e7e07236..230c8e18 100644 --- a/rust/timsseek/src/scoring/plan.rs +++ b/rust/timsseek/src/scoring/plan.rs @@ -138,6 +138,7 @@ pub struct FragmentIsotopeDecision { /// Owned by its reference library; callers cannot install a plan from another library. #[derive(Debug, Clone, Serialize)] pub struct ScoringPlan { + decoys: timsquery::models::capabilities::DecoyResolution, isotopes: IsotopePlan, fragment_isotopes: FragmentIsotopeDecision, #[serde(skip)] @@ -221,6 +222,7 @@ impl ScoringPlan { }) .collect(); Ok(Self { + decoys: *geom.decoy_resolution(), isotopes, fragment_isotopes, linear_indices, diff --git a/rust/timsseek/src/scoring/timings.rs b/rust/timsseek/src/scoring/timings.rs index 8fe720e0..7d13487c 100644 --- a/rust/timsseek/src/scoring/timings.rs +++ b/rust/timsseek/src/scoring/timings.rs @@ -181,7 +181,13 @@ pub enum RunStatus { /// Top-level report for an entire CLI invocation. /// Contains shared loading costs and per-file pipeline reports. #[derive(Debug, Default, Serialize)] -pub struct RunReport { +pub struct RunReport<'a> { + /// The owning library's resolved decisions; serialized once for all raw files. + #[serde(skip_serializing_if = "Option::is_none")] + pub scoring_plan: Option<&'a super::plan::ScoringPlan>, + /// Present only when a separate calibration library was supplied. + #[serde(skip_serializing_if = "Option::is_none")] + pub calibration_scoring_plan: Option<&'a super::plan::ScoringPlan>, /// Terminal status of the run. See [`RunStatus`]. pub status: RunStatus, /// When `status == Aborted`, a short human-readable reason. `None` on diff --git a/rust/timsseek_cli/src/search.rs b/rust/timsseek_cli/src/search.rs index 3fbf717c..b15663ae 100644 --- a/rust/timsseek_cli/src/search.rs +++ b/rust/timsseek_cli/src/search.rs @@ -548,6 +548,8 @@ pub(crate) fn search(args: &SearchArgs) -> std::result::Result<(), errors::CliEr None => (None, None, 0), }; + run_report.scoring_plan = Some(speclib.scoring_plan()); + run_report.calibration_scoring_plan = calib_lib.as_ref().map(|library| library.scoring_plan()); run_report.load_speclib_ms = load_speclib_ms; run_report.load_calib_lib_ms = load_calib_lib_ms; run_report.speclib_entries = speclib.len();