diff --git a/docs/development.md b/docs/development.md index 5a123c46..d8f3a096 100644 --- a/docs/development.md +++ b/docs/development.md @@ -33,6 +33,15 @@ Standalone `calib_dash` reads saved `calibration.json`, not a spectral library. [Analyte module documentation](../rust/timsquery/src/chemistry/analyte.rs) documents chemistry storage, reader mappings, library-wide sequence eligibility, and results format version 4. +Precursor isotope envelopes retain the three-bin C/S approximation. Scoring +finalization includes known modification C/S deltas or an explicitly based +molecular formula. If any stored target or decoy lacks usable counts, every +entry uses mass-estimated C/S. Generated mass-shift decoys reuse the parent's +envelope. The selected method and unavailable-count reasons appear in the +shared CLI/viewer plan report and Parquet scoring-plan metadata. Other elements +are ignored by this approximation; isotope-labelled C/S and unspecified formula +bases cannot supply its composition counts. + ## Cargo features | Feature | Crate | Effect | Use case | Enable | diff --git a/rust/timsquery/src/chemistry.rs b/rust/timsquery/src/chemistry.rs index 95da7744..d051ee13 100644 --- a/rust/timsquery/src/chemistry.rs +++ b/rust/timsquery/src/chemistry.rs @@ -2,7 +2,7 @@ //! //! Two crates need it and neither owns it: the mzSpecLib reader parses each //! analyte's ProForma through mzannotate, which takes an `&Ontologies`, and -//! timsseek's fallback sequence parser needs the same indexes. It lives here +//! timsseek's modification C/S resolver needs the same indexes. It lives here //! because timsseek depends on timsquery and not the reverse, so this is the //! lowest crate both can reach. //! @@ -15,10 +15,8 @@ use std::sync::OnceLock; /// The modification ontologies, built once on first use. /// /// GNOme is omitted. It adds 191,529 entries and 26.4 MB, taking -/// initialization from about 48 ms to 2.6 s, and a GNO modification already -/// produces no usable parse downstream: `count_carbon_sulphur_in_sequence` -/// rejects it during formula counting. PSI-MOD, XL-MOD, Unimod and RESID are -/// loaded because formula counts need them. +/// initialization from about 48 ms to 2.6 s. GNO accessions remain unresolved +/// with this ontology set. PSI-MOD, XL-MOD, Unimod and RESID are loaded. /// /// Note this is PSI-**MOD**, the protein-modification vocabulary, and not /// PSI-**MS**. Nothing here carries `MS:` terms, which is why the reader's @@ -27,8 +25,9 @@ use std::sync::OnceLock; /// mislead. /// /// Lazily built, and which libraries pay for it differs by format. A DIA-NN -/// library whose sequences all match the shared explicit-sequence parser never -/// reaches here. An mzSpecLib library always does: mzannotate takes +/// library can avoid it while parsing explicit sequences; scoring resolves +/// modification C/S contributions through these indexes. An mzSpecLib library +/// always reaches here: mzannotate takes /// `&Ontologies` to parse an analyte at all. pub fn ontologies() -> &'static mzcore::ontology::Ontologies { static ONTOLOGIES: OnceLock = OnceLock::new(); diff --git a/rust/timsquery/src/chemistry/analyte.rs b/rust/timsquery/src/chemistry/analyte.rs index a9bcf6d4..9a28697b 100644 --- a/rust/timsquery/src/chemistry/analyte.rs +++ b/rust/timsquery/src/chemistry/analyte.rs @@ -31,9 +31,10 @@ //! silently selecting one. Unsupported peptide structures preserve their annotation. //! Molecular formulas retain their declared basis (or `Unspecified`) and signed //! electron counts. Peptide and formula facts can coexist: sealing rejects conflicting -//! neutral formulas for unmodified canonical peptides. Comparison of modified or -//! ambiguous structures and ion-basis formulas is deferred to composition support; -//! coexistence alone does not certify chemical consistency. +//! neutral formulas for unmodified canonical peptides. Scoring finalization also +//! checks comparable neutral-formula C/S counts against modified peptides. Full +//! formula comparison for modified/ambiguous structures and ion-basis formulas +//! remains unsupported; coexistence alone does not certify chemical consistency. //! //! Search candidates carry a row and competition metadata, not copied sequences. //! Rescorers and the dashboard receive the owning `ReferenceLibrary`. Its library-wide @@ -42,8 +43,12 @@ //! one row lacking residues disables residue features for all targets and decoys. //! Disabled operations have no ML projections. Global/labile/ambiguous modifications //! remain unresolved where their complete set cannot be represented. The isotope -//! model still uses residues and its existing averagine fallback, not modification -//! or declared-formula composition. +//! model in timsseek includes modification C/S deltas or explicitly based formulas. +//! It selects composition-derived counts only when every stored target and decoy +//! supports them; otherwise all entries use mass-estimated C/S. Both paths retain +//! the same approximate C/S calculator. Synthetic variants reuse parent envelopes. +//! Labelled C/S isotopes and unspecified formula bases are unresolved for this model; +//! other elements do not contribute. This is not a full elemental isotope model. //! //! Peptide deduplication compares complete structural keys plus charge and m/z, //! then chooses the highest score, preferring a target on a tie. Equivalent diff --git a/rust/timsquery/src/models/capabilities.rs b/rust/timsquery/src/models/capabilities.rs index 83ed8d57..f88115d3 100644 --- a/rust/timsquery/src/models/capabilities.rs +++ b/rust/timsquery/src/models/capabilities.rs @@ -22,7 +22,8 @@ pub enum FragmentFeatureState { #[derive(Debug, Clone, Copy, PartialEq)] pub enum IsotopeStrategy { - /// Per-peptide: C/S countable -> composition envelope; else -> averagine. + /// Requested isotope bins. The scoring library resolves one C/S method + /// from whole-library coverage; this legacy name does not select a method. FromComposition { n_isotopes: u8 }, } diff --git a/rust/timsquery/src/models/target_columns.rs b/rust/timsquery/src/models/target_columns.rs index 3c77de19..c40ff117 100644 --- a/rust/timsquery/src/models/target_columns.rs +++ b/rust/timsquery/src/models/target_columns.rs @@ -358,6 +358,11 @@ impl TargetColumns { self.charge.len() } + /// Build a sidecar addressed by this arena's stored-row handles. + pub fn map_rows(&self, f: impl FnMut(RowIdx) -> T) -> RowValues { + RowValues(self.rows().map(f).collect()) + } + /// The rows of this arena, in storage order. pub fn rows(&self) -> impl Iterator + use { (0..self.n_rows() as u32).map(RowIdx::new) @@ -727,6 +732,26 @@ impl TargetColumns { } } +/// Dense sidecar for stored rows. Use only with handles from its owning arena. +#[derive(Clone)] +pub struct RowValues(Vec); + +impl std::ops::Index for RowValues { + type Output = T; + + fn index(&self, row: RowIdx) -> &T { + &self.0[row.get()] + } +} + +impl std::fmt::Debug for RowValues { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("RowValues") + .field("len", &self.0.len()) + .finish_non_exhaustive() + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/rust/timsquery_viewer/src/file_loader.rs b/rust/timsquery_viewer/src/file_loader.rs index 9e5e7f50..00e10d7d 100644 --- a/rust/timsquery_viewer/src/file_loader.rs +++ b/rust/timsquery_viewer/src/file_loader.rs @@ -167,7 +167,7 @@ impl FileLoader { /// `RefQuery` flyweights: the geometry feeds `build_extraction` + `TraceScorer` /// and the reference intensities + isotope envelope come from the flyweight's /// `ExpectedIntensity` impl, which routes the envelope through -/// `isotope_dist_or_averagine` (averagine fallback). There is no private +/// library-wide C/S plan (composition or mass-estimated counts). There is no private /// isotope model here -- the viewer shows exactly what the CLI scores. #[derive(Debug)] pub struct ElutionGroupData { @@ -445,10 +445,10 @@ mod tests { TargetColumnsBuilder, }; use timsquery::utils::constants::PROTON_MASS; - use timsseek::fragment_mass::isotope_dist_or_averagine; + use timsseek::fragment_mass::isotope_dist_from_mass; /// Build a one-entry `ReferenceLibrary` whose STRIPPED sequence has an - /// uncountable composition (`B` is not a real residue), forcing the + /// uncountable composition (`X` has unknown composition), forcing the /// averagine isotope path. fn uncountable_lib() -> ReferenceLibrary { let mut geom = TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); @@ -461,7 +461,7 @@ mod tests { (IonAnnot::try_from("y3").unwrap(), 300.0), (IonAnnot::try_from("y5").unwrap(), 500.0), ], - analyte: timsquery::chemistry::analyte::Analyte::from_sequence("PEPBK").as_input(), + analyte: timsquery::chemistry::analyte::Analyte::from_sequence("PEPXK").as_input(), ..Default::default() }); let geom = geom @@ -476,16 +476,16 @@ mod tests { /// The envelope the viewer DISPLAYS (obtained via `get_elem`, i.e. the /// flyweight's `expected_precursor_envelope`) must equal what the scoring - /// path computes via `isotope_dist_or_averagine` -- NOT the deleted + /// path computes via `isotope_dist_from_mass` -- NOT the deleted /// `[1.0, 0.0, 0.0]` fallback. #[test] - fn displayed_envelope_matches_isotope_dist_or_averagine() { + fn displayed_envelope_matches_isotope_dist_from_mass() { let data = ElutionGroupData::new(uncountable_lib()); let (_eg, expected) = data.get_elem(0).unwrap(); let charge = 1.0_f64; let neutral = 600.0 * charge - charge * PROTON_MASS; - let (_src, env) = isotope_dist_or_averagine("PEPBK", neutral); + let env = isotope_dist_from_mass(neutral); // Every displayed precursor isotope equals the shared model output. for (iso_idx, ref_intensity) in env.iter().enumerate() { diff --git a/rust/timsseek/src/data_sources/reference_library.rs b/rust/timsseek/src/data_sources/reference_library.rs index ae408c56..d03110ec 100644 --- a/rust/timsseek/src/data_sources/reference_library.rs +++ b/rust/timsseek/src/data_sources/reference_library.rs @@ -19,12 +19,11 @@ use timsquery::models::{ TargetColumns, }; use timsquery::traits::QueryGeom; +#[cfg(test)] use timsquery::utils::constants::PROTON_MASS; -use crate::fragment_mass::{ - IsotopeSource, - isotope_dist_or_averagine, -}; +#[cfg(test)] +use crate::fragment_mass::isotope_dist_from_mass; use crate::models::DecoyMarking; use crate::scoring::plan; @@ -177,7 +176,8 @@ impl ReferenceLibrary { ), }); } - let plan = plan::ScoringPlan::resolve(&geom); + let plan = plan::ScoringPlan::resolve(&geom) + .map_err(|message| TargetReadingError::InvalidLibrary { message })?; Ok(ReferenceLibrary { geom, frag_intens, @@ -251,16 +251,7 @@ impl<'a> ExpectedIntensity for RefQuery<'a> { fn expected_precursor_envelope(&self) -> SmallVec<[(i8, f32); 3]> { let tgt = self.geom.row(); let IsotopeStrategy::FromComposition { n_isotopes } = self.lib.geom.capabilities().isotopes; - let seq = self - .lib - .geom - .analyte(tgt) - .peptide - .known() - .map_or("", |p| p.residues); - let charge = self.lib.geom.charge(tgt) as f64; - let neutral = self.lib.geom.precursor_mz(tgt) * charge - charge * PROTON_MASS; - let (_src, env) = isotope_dist_or_averagine(seq, neutral); + let env = self.lib.plan.isotopes().envelope(tgt); (0..n_isotopes as usize) .map(|i| (i as i8, env[i])) .collect() @@ -379,7 +370,7 @@ fn retired_format(path: &Path) -> Option<&'static str> { /// Loading: the one path from a path on disk to a scored-against arena. impl ReferenceLibrary { /// Narrow a sealed [`TargetTable`] and finish it: decoy reporting, the - /// operation report, and the averagine tally. + /// operation and isotope-method report. /// /// The one definition of a finished library. `TargetTable`'s variants and /// fields are public, so a caller outside timsseek can assemble an arena @@ -481,36 +472,12 @@ impl ReferenceLibrary { } } - /// Report independently resolved operations and the existing isotope fallback tally. + /// Report the same library-wide decisions serialized with result metadata. fn report_scoring_plan(&self) { - let n_rows = self.geom.n_rows(); - let mut n_averagine_fallback = 0usize; - for tgt in self.geom.rows() { - let stripped = self - .geom - .analyte(tgt) - .peptide - .known() - .map_or("", |p| p.residues); - let charge = self.geom.charge(tgt) as f64; - let neutral_mass = self.geom.precursor_mz(tgt) * charge - charge * PROTON_MASS; - let (isotope_src, _envelope) = isotope_dist_or_averagine(stripped, neutral_mass); - if isotope_src == IsotopeSource::Averagine { - n_averagine_fallback += 1; - } - } - tracing::info!("{}", self.plan.summary()); for operation in self.plan.operations() { tracing::info!("{operation}"); } - if n_averagine_fallback > 0 { - tracing::warn!( - "{}/{} library entries used averagine isotope fallback", - n_averagine_fallback, - n_rows - ); - } } /// Log a one-line summary of the arena's shape at load time. @@ -723,7 +690,7 @@ mod tests { assert!(matches!(peptide.modifications, PropertyRef::Missing)); assert!(peptide.sequence().is_none()); assert!(!lib.all_sequence_counts_enabled()); - let expected = isotope_dist_or_averagine("PEPTIDE", 0.0).1; + let expected = isotope_dist_from_mass(500.0 * 2.0 - 2.0 * PROTON_MASS); for query in lib.iter() { for (i, intensity) in query.expected_precursor_envelope() { assert_eq!(intensity, expected[i as usize]); diff --git a/rust/timsseek/src/fragment_mass/averagine.rs b/rust/timsseek/src/fragment_mass/averagine.rs index 021c53a6..8f0ea65f 100644 --- a/rust/timsseek/src/fragment_mass/averagine.rs +++ b/rust/timsseek/src/fragment_mass/averagine.rs @@ -1,13 +1,5 @@ -use super::elution_group_converter::count_carbon_sulphur_in_sequence; use crate::isotopes::peptide_isotopes; -/// Which model produced an isotope envelope, for load-time reporting. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum IsotopeSource { - Composition, - Averagine, -} - // Senko averagine residue (avg amino acid): C4.9384 H7.7583 N1.3577 O1.4773 S0.0417, // average residue mass ~111.1054 Da. Per-Dalton element counts: const C_PER_DA: f64 = 4.9384 / 111.1054; @@ -22,26 +14,12 @@ pub fn averagine_cs_from_mass(neutral_mass: f64) -> (u16, u16) { /// Averagine isotope envelope: relative intensity, tallest peak == 1.0. /// -/// Matches `peptide_isotopes`'s own normalization (max peak, not sum), which -/// is what the `Composition` branch of `isotope_dist_or_averagine` returns -/// verbatim. Both isotope sources must share this scale so they're -/// interchangeable at scoring time. +/// Uses the same C/S calculator and normalization as composition-derived counts. pub fn isotope_dist_from_mass(neutral_mass: f64) -> [f32; 3] { let (c, s) = averagine_cs_from_mass(neutral_mass); peptide_isotopes(c, s) } -/// Composition envelope when the sequence is countable, else averagine from mass. -pub fn isotope_dist_or_averagine(seq: &str, neutral_mass: f64) -> (IsotopeSource, [f32; 3]) { - match count_carbon_sulphur_in_sequence(seq) { - Ok((c, s)) => (IsotopeSource::Composition, peptide_isotopes(c, s)), - Err(_) => ( - IsotopeSource::Averagine, - isotope_dist_from_mass(neutral_mass), - ), - } -} - #[cfg(test)] mod tests { use super::*; @@ -65,22 +43,4 @@ mod tests { "env {env:?} has a value outside [0, 1]" ); } - - #[test] - fn or_averagine_uses_composition_for_standard_peptide() { - let (src, _env) = isotope_dist_or_averagine("PEPTIDEK", 900.4); - assert_eq!(src, IsotopeSource::Composition); - } - - #[test] - fn or_averagine_falls_back_on_nonstandard() { - // mzcore treats `B`, also called Asx, as ambiguous between Asp and Asn. - // That gives the formula path multiple results and makes it return an - // error. `X` does not exercise this path because mzcore assigns it a - // zero-C/S formula. - let (src, env) = isotope_dist_or_averagine("PEPBK", 600.0); - assert_eq!(src, IsotopeSource::Averagine); - let max = env.iter().copied().fold(f32::MIN, f32::max); - assert!((max - 1.0).abs() < 1e-4, "env {env:?} max is {max}"); - } } diff --git a/rust/timsseek/src/fragment_mass/elution_group_converter.rs b/rust/timsseek/src/fragment_mass/elution_group_converter.rs index 8f63d06b..ab4e7510 100644 --- a/rust/timsseek/src/fragment_mass/elution_group_converter.rs +++ b/rust/timsseek/src/fragment_mass/elution_group_converter.rs @@ -1,4 +1,4 @@ -use crate::isotopes::peptide_isotopes; +#[cfg(test)] use mzcore::prelude::{ AmbiguousMolecule, Element, @@ -33,6 +33,7 @@ pub fn supersimpleprediction(mz: f64, charge: i32) -> f64 { + (1.176651e-01 * charge as f64) } +#[cfg(test)] fn count_carbon_sulphur(form: &MolecularFormula) -> (u16, u16) { let mut ncarbon = 0; let mut nsulphur = 0; @@ -91,7 +92,7 @@ const RESIDUE_CS: [Option<(u16, u16)>; 26] = { /// Fast (C, S) tally over a bare amino-acid sequence via [`RESIDUE_CS`]. /// Empty input and non-standard residues return `None`, so mzcore handles them. -fn count_cs_fast(sequence: &str) -> Option<(u16, u16)> { +pub(crate) fn count_cs_fast(sequence: &str) -> Option<(u16, u16)> { if sequence.is_empty() { return None; } @@ -100,22 +101,13 @@ fn count_cs_fast(sequence: &str) -> Option<(u16, u16)> { for &b in sequence.as_bytes() { let idx = b.wrapping_sub(b'A') as usize; let (c, s) = RESIDUE_CS.get(idx).copied().flatten()?; - ncarbon += c; - nsulphur += s; + ncarbon = ncarbon.checked_add(c)?; + nsulphur = nsulphur.checked_add(s)?; } Some((ncarbon, nsulphur)) } -/// (C, S) counts for `sequence` (a bare, mod-stripped peptide on the hot path). -/// Try the allocation-free table first. Use mzcore for empty or non-standard -/// input, where it remains the authority. -pub fn count_carbon_sulphur_in_sequence(sequence: &str) -> Result<(u16, u16), String> { - if let Some(cs) = count_cs_fast(sequence) { - return Ok(cs); - } - count_carbon_sulphur_in_sequence_mzcore(sequence) -} - +#[cfg(test)] fn count_carbon_sulphur_in_sequence_mzcore(sequence: &str) -> Result<(u16, u16), String> { let peptide = crate::models::sequence::parse_proforma(sequence) .map_err(|e| format!("Error parsing peptide sequence {sequence}: {e}"))?; @@ -133,11 +125,6 @@ fn count_carbon_sulphur_in_sequence_mzcore(sequence: &str) -> Result<(u16, u16), Ok(count_carbon_sulphur(&form)) } -pub fn isotope_dist_from_seq(sequence: &str) -> Result<[f32; 3], String> { - let (ncarbon, nsulphur) = count_carbon_sulphur_in_sequence(sequence)?; - Ok(peptide_isotopes(ncarbon, nsulphur)) -} - #[cfg(test)] mod tests { use super::*; diff --git a/rust/timsseek/src/fragment_mass/isotope_plan.rs b/rust/timsseek/src/fragment_mass/isotope_plan.rs new file mode 100644 index 00000000..8ba7368e --- /dev/null +++ b/rust/timsseek/src/fragment_mass/isotope_plan.rs @@ -0,0 +1,538 @@ +//! One C/S envelope method for the entire library, resolved before extraction. +//! +//! The existing three-bin C/S approximation is retained; this is not a full +//! elemental isotope calculation. Modifications contribute signed C/S deltas. +//! A mass-only annotation cannot establish those deltas. Generated mass shifts +//! reuse their stored parent's envelope, including when counts come from mass. +use std::collections::{ + BTreeMap, + HashMap, +}; + +use mzcore::prelude::{ + AmbiguousMolecule, + AminoAcid, + Element, + MolecularFormula, + Molecule, +}; +use mzcore::sequence::SimpleModificationInner; +use serde::Serialize; +use timsquery::IonAnnot; +use timsquery::chemistry::analyte::{ + AnalyteRef, + FormulaBasis, + PeptideRef, +}; +use timsquery::chemistry::ontologies; +use timsquery::models::target_columns::RowValues; +use timsquery::models::{ + RowIdx, + TargetColumns, +}; +use timsquery::utils::constants::PROTON_MASS; + +use super::averagine::isotope_dist_from_mass; +use crate::isotopes::peptide_isotopes; + +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] +#[serde(rename_all = "snake_case")] +pub enum IsotopeMethod { + CompositionCs, + MassEstimatedCs, +} + +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Serialize)] +#[serde(rename_all = "snake_case")] +pub enum UnavailableReason { + MissingChemistry, + IncompleteModifications, + UnresolvedModification, + UnresolvedResidues, + UnsupportedFormula, + InvalidCounts, + ModelRange, +} + +type Counts = (i64, i64); +type Resolution = Result; + +/// Reported with the scoring plan; cached envelopes never enter feature metadata. +#[derive(Debug, Clone, Serialize)] +pub struct IsotopePlan { + pub method: IsotopeMethod, + pub composition_rows: usize, + pub total_rows: usize, + pub unavailable: BTreeMap, + #[serde(skip)] + envelopes: RowValues<[f32; 3]>, +} + +impl IsotopePlan { + pub(crate) fn resolve(geom: &TargetColumns) -> Result { + let mut resolver = Resolver::default(); + let mut unavailable = BTreeMap::new(); + let mut composition_rows = 0; + let mut conflict = None; + let counts = geom.map_rows(|row| { + let result = resolver.analyte(geom.analyte(row)); + match result { + Ok(Ok(cs)) => { + composition_rows += 1; + Some(cs) + } + Ok(Err(reason)) => { + *unavailable.entry(reason).or_insert(0) += 1; + None + } + Err(message) => { + conflict = Some(format!("entry {:?}: {message}", geom.output_id(row))); + None + } + } + }); + if let Some(message) = conflict { + return Err(message); + } + let total_rows = geom.n_rows(); + let method = if composition_rows == total_rows { + IsotopeMethod::CompositionCs + } else { + IsotopeMethod::MassEstimatedCs + }; + let mut invalid_mass = None; + let envelopes = geom.map_rows(|row| match method { + IsotopeMethod::CompositionCs => { + let (c, s) = counts[row].expect("library plan promised C/S counts"); + peptide_isotopes(c as u16, s as u16) + } + IsotopeMethod::MassEstimatedCs => { + // Stored geometry, never the synthetic variant's shifted mass. + let z = f64::from(geom.charge(row)); + let envelope = isotope_dist_from_mass(geom.precursor_mz(row) * z - z * PROTON_MASS); + if !envelope.iter().all(|v| v.is_finite()) { + invalid_mass = Some(format!( + "entry {}: precursor mass exceeds the C/S isotope model's numerical range", + geom.output_id(row) + )); + } + envelope + } + }); + if let Some(message) = invalid_mass { + return Err(message); + } + Ok(Self { + method, + composition_rows, + total_rows, + unavailable, + envelopes, + }) + } + + pub(crate) fn envelope(&self, row: RowIdx) -> &[f32; 3] { + &self.envelopes[row] + } +} + +impl std::fmt::Display for IsotopePlan { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let method = match self.method { + IsotopeMethod::CompositionCs => "composition-derived C/S", + IsotopeMethod::MassEstimatedCs => "mass-estimated C/S (averagine)", + }; + write!( + f, + "Precursor isotopes: {method} for all {} entries ({}/{} with usable composition counts)", + self.total_rows, self.composition_rows, self.total_rows + )?; + for (reason, count) in &self.unavailable { + let reason = match reason { + UnavailableReason::MissingChemistry => "missing or unresolved chemistry", + UnavailableReason::IncompleteModifications => "incomplete modifications", + UnavailableReason::UnresolvedModification => "unresolved modification composition", + UnavailableReason::UnresolvedResidues => "unresolved residue C/S counts", + UnavailableReason::UnsupportedFormula => { + "unsupported formula basis or labelled C/S" + } + UnavailableReason::InvalidCounts => "invalid C/S counts", + UnavailableReason::ModelRange => "outside the isotope model's numerical range", + }; + write!(f, "; {count} {reason}")?; + } + Ok(()) + } +} + +#[derive(Default)] +struct Resolver { + modifications: HashMap, +} + +impl Resolver { + // Conflicting comparable declarations are invalid input, not absent chemistry. + fn analyte(&mut self, analyte: AnalyteRef<'_>) -> Result { + let peptide = analyte.peptide.known().map(|p| self.peptide(p)); + let formula = analyte.formula.known().map(|f| { + if f.basis == FormulaBasis::Unspecified { + return Err(UnavailableReason::UnsupportedFormula); + } + elements_cs(&f.elements).and_then(valid_counts) + }); + if let (Some(Ok(a)), Some(Ok(b))) = (peptide, formula) + && analyte + .formula + .known() + .is_some_and(|f| f.basis == FormulaBasis::NeutralMolecule) + && a != b + { + return Err( + "declared neutral formula disagrees with peptide/modification C/S counts".into(), + ); + } + // An explicitly supplied formula can establish composition without sequence. + // An observed-ion formula already includes adduct atoms; don't add them again. + if let Some(result) = formula { + return Ok(result); + } + Ok(peptide.unwrap_or(Err(UnavailableReason::MissingChemistry))) + } + + fn peptide(&mut self, peptide: PeptideRef<'_>) -> Resolution { + let modifications = peptide + .modifications + .known() + .ok_or(UnavailableReason::IncompleteModifications)?; + let mut counts = (0i64, 0i64); + if let Some((c, s)) = super::elution_group_converter::count_cs_fast(peptide.residues) { + counts = (i64::from(c), i64::from(s)); + } else { + for residue in peptide.residues.bytes() { + // X's placeholder zero formula is not evidence of zero C/S. + if residue == b'X' { + return Err(UnavailableReason::UnresolvedResidues); + } + let aa = AminoAcid::try_from(residue as char) + .map_err(|_| UnavailableReason::UnresolvedResidues)?; + let formulas = aa.formulas(); + let mut alternatives = formulas.iter().map(formula_cs); + let cs = alternatives + .next() + .ok_or(UnavailableReason::UnresolvedResidues)??; + if !alternatives.all(|v| v == Ok(cs)) { + return Err(UnavailableReason::UnresolvedResidues); + } + counts.0 += cs.0; + counts.1 += cs.1; + } + } + for (_, modification) in modifications.iter() { + let annotation = modification.annotation(); + let delta = if let Some(&delta) = self.modifications.get(annotation) { + delta + } else { + let delta = modification_cs(annotation); + self.modifications.insert(annotation.to_owned(), delta); + delta + }; + let (c, s) = delta?; + counts.0 += c; + counts.1 += s; + } + valid_counts(counts) + } +} + +fn modification_cs(annotation: &str) -> Resolution { + let ((parsed, _), _) = SimpleModificationInner::pro_forma( + annotation, + &mut Default::default(), + &mut Default::default(), + ontologies(), + ) + .map_err(|_| UnavailableReason::UnresolvedModification)?; + let modification = parsed + .defined() + .ok_or(UnavailableReason::UnresolvedModification)?; + if matches!( + modification.as_ref(), + SimpleModificationInner::Mass(..) | SimpleModificationInner::Info(..) + ) { + return Err(UnavailableReason::UnresolvedModification); + } + formula_cs(&modification.formula()) +} + +fn formula_cs(formula: &MolecularFormula) -> Resolution { + if *formula.additional_mass() != 0.0 { + return Err(UnavailableReason::UnresolvedModification); + } + elements_cs(formula.elements()) +} + +fn elements_cs(elements: &[(Element, Option, i32)]) -> Resolution { + let mut counts = (0, 0); + for &(element, isotope, count) in elements { + if matches!(element, Element::C | Element::S) && isotope.is_some() { + return Err(UnavailableReason::UnsupportedFormula); + } + match element { + Element::C => counts.0 += i64::from(count), + Element::S => counts.1 += i64::from(count), + _ => {} + } + } + Ok(counts) +} + +fn valid_counts(cs: Counts) -> Resolution { + if !(0..=i64::from(u16::MAX)).contains(&cs.0) || !(0..=i64::from(u16::MAX)).contains(&cs.1) { + Err(UnavailableReason::InvalidCounts) + } else if !peptide_isotopes(cs.0 as u16, cs.1 as u16) + .iter() + .all(|v| v.is_finite()) + { + Err(UnavailableReason::ModelRange) + } else { + Ok(cs) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::data_sources::reference_library::{ + ExpectedIntensity, + ReferenceLibrary, + }; + use timsquery::chemistry::analyte::{ + Analyte, + Formula, + Property, + }; + use timsquery::models::capabilities::{ + DecoyPolicy, + TargetCapabilities, + }; + use timsquery::models::{ + Row, + TargetColumnsBuilder, + }; + use timsquery::serde::TargetTable; + + fn library( + analytes: &[Analyte], + shipped_decoy: bool, + ) -> Result { + let mut builder = + TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); + for (i, analyte) in analytes.iter().enumerate() { + builder.push_row(Row { + precursor_mz: 500.0 + i as f64 * 100.0, + charge: 2, + frags: &[(IonAnnot::try_from("y3").unwrap(), 300.0)], + analyte: analyte.as_input(), + is_decoy: shipped_decoy && i == 1, + ..Default::default() + }); + } + ReferenceLibrary::try_from(TargetTable::Mzpaf { + geom: builder + .seal(if shipped_decoy { + DecoyPolicy::Never + } else { + DecoyPolicy::Force + }) + .unwrap(), + frag_intens: Some(vec![1.0; analytes.len()]), + }) + } + + fn envelope(lib: &ReferenceLibrary, row: usize, variant: u8) -> Vec { + let geom = lib.geometry(); + let row = geom.rows().nth(row).unwrap(); + lib.item_at(geom.flat_for(row, variant)) + .expected_precursor_envelope() + .iter() + .map(|(_, i)| *i) + .collect() + } + + #[test] + fn modifications_contribute_signed_cs_and_other_atoms_leave_model_unchanged() { + // PEPTCIDEK has C43 S1; carbamidomethyl adds C2, sulfur adds S1. + for (sequence, expected) in [ + ("PEPTCIDEK", (43, 1)), + ("PEPTC[UNIMOD:4]IDEK", (45, 1)), + ("PEPTC[Formula:C2H3NO]IDEK", (45, 1)), + ("PEPTS[UNIMOD:21]IDE", (37, 0)), + ("PEPTS[MOD:00046]IDE", (37, 0)), + ("PEPBK", (25, 0)), // D/N alternatives agree on C/S. + ("PEPDK", (25, 0)), + ("[UNIMOD:1]-PEPTCIDEK", (45, 1)), + ("PEPTC[Formula:S-1]IDEK", (43, 0)), + ("PEPTC[Formula:S]IDEK", (43, 2)), + ("PEPTC[Formula:C-1]IDEK", (42, 1)), + ("PEPTC[Formula:O]IDEK", (43, 1)), + ] { + let lib = library(&[Analyte::from_sequence(sequence)], false).unwrap(); + assert_eq!( + lib.scoring_plan().isotopes().method, + IsotopeMethod::CompositionCs, + "{sequence}" + ); + assert_eq!( + envelope(&lib, 0, 0), + peptide_isotopes(expected.0, expected.1) + ); + } + } + + #[test] + fn one_mass_only_entry_selects_mass_for_every_row_and_variant() { + for shipped in [false, true] { + let lib = library( + &[ + Analyte::from_sequence("PEPTC[UNIMOD:4]IDEK"), + Analyte::from_sequence("PEPTC[+57.021]IDEK"), + ], + shipped, + ) + .unwrap(); + let plan = lib.scoring_plan().isotopes(); + assert_eq!(plan.method, IsotopeMethod::MassEstimatedCs); + assert_eq!(plan.composition_rows, 1); + assert_eq!( + plan.unavailable[&UnavailableReason::UnresolvedModification], + 1 + ); + for row in 0..2 { + let expected = + isotope_dist_from_mass((500.0 + row as f64 * 100.0) * 2.0 - 2.0 * PROTON_MASS); + for variant in 0..if shipped { 1 } else { 3 } { + assert_eq!(envelope(&lib, row, variant), expected); + } + } + } + } + + #[test] + fn formula_only_neutral_and_ion_inputs_need_no_sequence() { + for basis in [FormulaBasis::NeutralMolecule, FormulaBasis::ObservedIon] { + let analyte = Analyte::from_formula(&mzcore::molecular_formula!(C 12 H 20 S 2), basis); + let lib = library(&[analyte], false).unwrap(); + assert_eq!( + lib.scoring_plan().isotopes().method, + IsotopeMethod::CompositionCs + ); + assert_eq!(envelope(&lib, 0, 0), peptide_isotopes(12, 2)); + assert_eq!(lib.scoring_plan().width(), 0); + } + } + + #[test] + fn invalid_unknown_and_labelled_composition_select_mass() { + for analyte in [ + Analyte::default(), + Analyte::from_sequence("PEPXIDEK"), + Analyte::from_sequence("<13C>PEPTIDEK"), + Analyte::from_sequence("PEPTC[Formula:C-100]IDEK"), + Analyte::from_sequence("PEPTC[Formula:[13C2]]IDEK"), + Analyte::from_formula( + &mzcore::molecular_formula!(C 10 H 20), + FormulaBasis::Unspecified, + ), + Analyte { + formula: Property::Known(Formula { + elements: vec![(Element::C, std::num::NonZeroU16::new(13), 2)], + basis: FormulaBasis::NeutralMolecule, + }), + ..Default::default() + }, + ] { + let lib = library(std::slice::from_ref(&analyte), false).unwrap(); + assert_eq!( + lib.scoring_plan().isotopes().method, + IsotopeMethod::MassEstimatedCs, + "{analyte:?}" + ); + } + } + + #[test] + fn comparable_formula_conflicts_are_rejected_and_agreement_is_not_double_counted() { + let mut analyte = Analyte::from_sequence("PEPTC[UNIMOD:4]IDEK"); + for (carbon, valid) in [(45, true), (44, false)] { + analyte.formula = Property::Known(Formula { + elements: vec![(Element::C, None, carbon), (Element::S, None, 1)], + basis: FormulaBasis::NeutralMolecule, + }); + let result = library(std::slice::from_ref(&analyte), false); + if valid { + assert_eq!(envelope(&result.unwrap(), 0, 0), peptide_isotopes(45, 1)); + } else { + assert!(format!("{:?}", result.unwrap_err()).contains("disagrees")); + } + } + } + + #[test] + fn numerical_model_limits_are_checked() { + assert_eq!( + valid_counts((65536, 0)), + Err(UnavailableReason::InvalidCounts) + ); + assert_eq!(valid_counts((65535, 0)), Err(UnavailableReason::ModelRange)); + let mut builder = + TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); + builder.push_row(Row { + precursor_mz: 1e9, + charge: 2, + frags: &[(IonAnnot::try_from("y3").unwrap(), 300.0)], + ..Default::default() + }); + let result = ReferenceLibrary::try_from(TargetTable::Mzpaf { + geom: builder.seal(DecoyPolicy::Never).unwrap(), + frag_intens: Some(vec![1.0]), + }); + assert!( + matches!(result, Err(crate::errors::TargetReadingError::InvalidLibrary { message }) if message.contains("numerical range")) + ); + } + + #[test] + fn diann_and_mzspeclib_modification_annotations_reach_the_same_envelope() { + let diann = concat!( + "ModifiedPeptide\tStrippedPeptide\tPrecursorMz\tPrecursorCharge\tTr_recalibrated\tIonMobility\t", + "ProteinID\tDecoy\tFragmentMz\tFragmentType\tFragmentNumber\tFragmentCharge\tFragmentLossType\tRelativeIntensity\n", + "PEPTC(UniMod:4)IDEK\tPEPTCIDEK\t500\t2\t1\t1\tP1\t0\t300\ty\t3\t1\tnoloss\t1\n", + ); + let mzspeclib = concat!( + "\nMS:1003186|library format version=1.0\nMS:1003188|library name=isotopes\n", + "\nMS:1003061|library spectrum name=independent label\n", + "MS:1000744|selected ion m/z=500\nMS:1000041|charge state=2\n", + "MS:1003059|number of peaks=1\n\n", + "MS:1003270|proforma peptidoform ion notation=PEPTC[U:Carbamidomethyl]IDEK/2\n", + "\n300\t1\ty3/0.0\n", + ); + for (suffix, text) in [(".tsv", diann), (".mzspeclib.txt", mzspeclib)] { + let file = tempfile::Builder::new().suffix(suffix).tempfile().unwrap(); + std::fs::write(file.path(), text).unwrap(); + let lib = ReferenceLibrary::from_file(file.path(), Default::default()) + .unwrap_or_else(|e| panic!("{suffix}: {e:?}")); + assert_eq!( + lib.scoring_plan().isotopes().method, + IsotopeMethod::CompositionCs + ); + for query in lib.iter() { + let values: Vec<_> = query + .expected_precursor_envelope() + .iter() + .map(|(_, i)| *i) + .collect(); + assert_eq!(values, peptide_isotopes(45, 1), "{suffix}"); + } + } + } +} diff --git a/rust/timsseek/src/fragment_mass/mod.rs b/rust/timsseek/src/fragment_mass/mod.rs index 550b8ac2..7307793a 100644 --- a/rust/timsseek/src/fragment_mass/mod.rs +++ b/rust/timsseek/src/fragment_mass/mod.rs @@ -1,8 +1,6 @@ pub mod averagine; pub mod elution_group_converter; -pub use averagine::{ - IsotopeSource, - isotope_dist_from_mass, - isotope_dist_or_averagine, -}; +pub use averagine::isotope_dist_from_mass; + +pub mod isotope_plan; diff --git a/rust/timsseek/src/models/sequence.rs b/rust/timsseek/src/models/sequence.rs index caba84e0..3fbecabf 100644 --- a/rust/timsseek/src/models/sequence.rs +++ b/rust/timsseek/src/models/sequence.rs @@ -1,9 +1,10 @@ -//! Sequence-feature column names and the isotope helper's parsing fallback. +//! Sequence-feature column names. +#[cfg(test)] use timsquery::chemistry::ontologies; -/// Ontology-backed fallback for the existing isotope composition helper. -/// Search sequence features read stored analyte structure instead. +/// Test oracle for composition and modification masses. +#[cfg(test)] pub fn parse_proforma( sequence: &str, ) -> Result, String> { diff --git a/rust/timsseek/src/scoring/parquet_writer.rs b/rust/timsseek/src/scoring/parquet_writer.rs index 388ccc7f..b4656287 100644 --- a/rust/timsseek/src/scoring/parquet_writer.rs +++ b/rust/timsseek/src/scoring/parquet_writer.rs @@ -520,6 +520,9 @@ mod tests { .expect("plan metadata"); let plan: serde_json::Value = serde_json::from_str(plan.value.as_deref().unwrap()).unwrap(); assert_eq!(plan["rows"], 1); + assert_eq!(plan["isotopes"]["method"], "composition_cs"); + assert_eq!(plan["isotopes"]["composition_rows"], 1); + assert!(plan["isotopes"].get("envelopes").is_none()); assert_eq!(plan["unmodified_rows"], 1); assert_eq!(plan["operations"][0]["requirement"], "residue_sequence"); assert_eq!( diff --git a/rust/timsseek/src/scoring/plan.rs b/rust/timsseek/src/scoring/plan.rs index ba19f603..d1b7ddd9 100644 --- a/rust/timsseek/src/scoring/plan.rs +++ b/rust/timsseek/src/scoring/plan.rs @@ -19,6 +19,7 @@ use super::blocks::{ NameSink, ScoreBlock, }; +use crate::fragment_mass::isotope_plan::IsotopePlan; use serde::Serialize; use std::sync::Arc; use timsquery::chemistry::analyte::{ @@ -123,12 +124,14 @@ impl std::fmt::Display for OperationDecision { /// Owned by its reference library; callers cannot install a plan from another library. #[derive(Debug, Clone, Serialize)] pub struct ScoringPlan { + isotopes: IsotopePlan, rows: usize, unmodified_rows: usize, operations: Vec, } impl ScoringPlan { - pub(crate) fn resolve(geom: &TargetColumns) -> Self { + pub(crate) fn resolve(geom: &TargetColumns) -> Result { + let isotopes = IsotopePlan::resolve(geom)?; let mut unmodified_rows = 0; let mut residues = Coverage::default(); let mut modifications = Coverage::default(); @@ -172,21 +175,26 @@ impl ScoringPlan { } }) .collect(); - Self { + Ok(Self { + isotopes, rows, unmodified_rows, operations, - } + }) } /// 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 + "{} library entries; {} unmodified (known empty modification list); {}", + self.rows, self.unmodified_rows, self.isotopes ) } + pub fn isotopes(&self) -> &IsotopePlan { + &self.isotopes + } + pub fn unmodified_rows(&self) -> usize { self.unmodified_rows } @@ -268,7 +276,7 @@ mod tests { &"A".repeat(300), ] { let geom = arena(&Analyte::from_sequence(sequence)); - let plan = ScoringPlan::resolve(&geom); + let plan = ScoringPlan::resolve(&geom).unwrap(); assert!(plan.enabled(Operation::ResidueCounts), "{sequence}"); assert!(plan.enabled(Operation::ModificationCounts), "{sequence}"); assert_eq!(plan.width(), 22); @@ -284,7 +292,7 @@ mod tests { Analyte::from_sequence_fields("PEP[unresolved]TIDE", "PEPTIDE").unwrap(), ] { let geom = arena(&analyte); - let plan = ScoringPlan::resolve(&geom); + let plan = ScoringPlan::resolve(&geom).unwrap(); assert!(plan.enabled(Operation::ResidueCounts)); assert!(!plan.enabled(Operation::ModificationCounts)); assert_eq!(plan.width(), 21); @@ -314,7 +322,7 @@ mod tests { peptide, ..Default::default() }); - let plan = ScoringPlan::resolve(&geom); + let plan = ScoringPlan::resolve(&geom).unwrap(); assert_eq!(plan.width(), 0); assert_eq!(plan.unmodified_rows, 99); let coverage = &plan.operations()[0].coverage; @@ -344,7 +352,7 @@ mod tests { ..Default::default() }); let geom = builder.seal(DecoyPolicy::Never).unwrap(); - let plan = ScoringPlan::resolve(&geom); + let plan = ScoringPlan::resolve(&geom).unwrap(); assert_eq!(plan.width(), 0); assert_eq!(plan.operations()[0].coverage.missing, 1); } @@ -367,7 +375,7 @@ mod tests { }, ..Default::default() }); - let plan = ScoringPlan::resolve(&geom); + let plan = ScoringPlan::resolve(&geom).unwrap(); assert_eq!(plan.operations()[0].coverage.recovered, 1); assert!( plan.operations()[0] @@ -388,7 +396,7 @@ mod tests { #[should_panic(expected = "plan promised residues")] fn broken_plan_promise_is_an_invariant_failure() { let geom = arena(&Analyte::from_sequence("PEPTIDE")); - ScoringPlan::resolve(&geom).project( + ScoringPlan::resolve(&geom).unwrap().project( || AnalyteRef { peptide: PropertyRef::Missing, formula: PropertyRef::Missing, diff --git a/rust/timsseek_cli/src/predicted_library.rs b/rust/timsseek_cli/src/predicted_library.rs index 1caa1faf..de116b02 100644 --- a/rust/timsseek_cli/src/predicted_library.rs +++ b/rust/timsseek_cli/src/predicted_library.rs @@ -765,7 +765,7 @@ mod tests { /// The other route: msspeculator's own mzSpecLib writer to a file, then this /// project's reader back off it, which is what a `build-library` followed by /// a `search` does. - fn via_file(rows: &[SpectrumRow<'_>]) -> TargetColumns { + fn via_file_table(rows: &[SpectrumRow<'_>]) -> TargetTable { let dir = tempfile::tempdir().expect("temp dir"); // The name the sniffer dispatches on: it takes `.mzspeclib.` anywhere in // the file name, and reads plain text for anything not ending `.gz`. @@ -781,19 +781,48 @@ mod tests { } writer.finish().expect("file finishes"); } - let TargetTable::Mzpaf { geom, .. } = timsquery::serde::read_targets_with( + timsquery::serde::read_targets_with( &path, timsseek::LoadPolicy { decoys: DecoyPolicy::Never, ..Default::default() }, ) - .expect("the file this project writes reads back") else { + .expect("the file this project writes reads back") + } + + fn via_file(rows: &[SpectrumRow<'_>]) -> TargetColumns { + let TargetTable::Mzpaf { geom, .. } = via_file_table(rows) else { panic!("an mzSpecLib library is mzpaf-labelled"); }; geom } + #[test] + fn modified_prediction_and_reload_share_composition_envelopes() { + use timsseek::data_sources::reference_library::ExpectedIntensity; + use timsseek::fragment_mass::isotope_plan::IsotopeMethod; + let peptide = Fixture::new("PEPTC[UNIMOD:4]IDEK", "PEPTCIDEK"); + let rows = [peptide.row(2, false, None, peaks(4))]; + let predicted = build(&rows, DecoyPolicy::Never); + let reloaded = ReferenceLibrary::try_from(via_file_table(&rows)).unwrap(); + for lib in [&predicted.library, &reloaded] { + assert_eq!( + lib.scoring_plan().isotopes().method, + IsotopeMethod::CompositionCs + ); + let envelope: Vec<_> = lib + .iter() + .next() + .unwrap() + .expected_precursor_envelope() + .iter() + .map(|(_, value)| *value) + .collect(); + assert_eq!(envelope, timsseek::isotopes::peptide_isotopes(45, 1)); + } + } + /// This row's fragments as label to m/z, which is how the two routes can be /// compared at all: the writer sorts peaks by m/z and the sink keeps /// `(position, ion type)`, so the sequences differ where the mapping does not.