diff --git a/Cargo.lock b/Cargo.lock index 2eb4c5b7..865d7112 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4019,9 +4019,10 @@ dependencies = [ name = "micromzpaf" version = "0.33.0" dependencies = [ - "rustyms", "serde", + "serde_json", "thiserror 2.0.18", + "tracing", ] [[package]] diff --git a/Taskfile.yml b/Taskfile.yml index 676dc1d6..f4459c3d 100644 --- a/Taskfile.yml +++ b/Taskfile.yml @@ -30,6 +30,11 @@ tasks: # dashboard.rs only exists under the opt-in feature; CI only `cargo check`s it. - cargo clippy -p timsseek_cli --features dashboard --all-targets -- -D warnings + # Keep workspace documentation warning-free. + doc: + cmds: + - RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps + todos: cmds: - grep -nH -R TODO rust/tims*/src @@ -81,4 +86,3 @@ tasks: desc: Stop local Koina server cmds: - docker stop koina-local - diff --git a/rust/apex_sim/src/plots.rs b/rust/apex_sim/src/plots.rs index 79015208..936f2645 100644 --- a/rust/apex_sim/src/plots.rs +++ b/rust/apex_sim/src/plots.rs @@ -278,7 +278,7 @@ impl Default for TraceToggles { } /// Bottom panel: the apex-finder intermediate traces, min-max normalized to -/// [0,1] so heterogeneous scales overlay legibly. Apex markers overlaid. +/// `[0,1]` so heterogeneous scales overlay legibly. Apex markers overlaid. pub fn traces(ui: &mut egui::Ui, score: &ScoreResult, tog: &TraceToggles, true_apex: f32) { let t = &score.traces; let mut series: Vec<(&str, &[f32])> = Vec::new(); diff --git a/rust/micromzpaf/Cargo.toml b/rust/micromzpaf/Cargo.toml index 83c47f4c..41434f29 100644 --- a/rust/micromzpaf/Cargo.toml +++ b/rust/micromzpaf/Cargo.toml @@ -7,6 +7,9 @@ license.workspace = true [dependencies] serde = { workspace = true, features = ["derive"] } thiserror = { workspace = true } +# Used for the one-time warning when parsing drops a mass-error suffix. +tracing = { workspace = true } -# Workspace-inherited deps -rustyms = { workspace = true } +[dev-dependencies] +# Pins the string serialization format. +serde_json = { workspace = true } diff --git a/rust/micromzpaf/src/error.rs b/rust/micromzpaf/src/error.rs new file mode 100644 index 00000000..d8b1d241 --- /dev/null +++ b/rust/micromzpaf/src/error.rs @@ -0,0 +1,52 @@ +//! Why an annotation could not be parsed or represented. + +use thiserror::Error; + +use crate::{ + CHARGE_MAX, + CHARGE_MIN, + ISOTOPE_MAX, + ISOTOPE_MIN, +}; + +/// Why an annotation could not be parsed or represented. +#[derive(Debug, Error)] +pub enum IonParsingError { + #[error("Ordinal {ordinal} out of range for series '{series}'")] + OrdinalOutOfRange { ordinal: u8, series: char }, + #[error("Series '{series}' requires an ordinal")] + MissingOrdinal { series: char }, + #[error("Series '{series}' takes no ordinal, got {ordinal}")] + UnexpectedOrdinal { series: char, ordinal: u8 }, + #[error("Unsupported fragment type: '{fragment_type}'")] + UnsupportedFragmentType { fragment_type: char }, + #[error("Charge cannot be 0")] + ChargeCannotBeZero, + #[error("Charge {charge} outside the representable range {CHARGE_MIN}..={CHARGE_MAX}")] + ChargeOutOfRange { charge: i8 }, + #[error( + "Isotope offset {isotope} outside the representable range {ISOTOPE_MIN}..={ISOTOPE_MAX}" + )] + IsotopeOutOfRange { isotope: i8 }, + #[error("Neutral loss '{loss}' is not representable")] + UnsupportedNeutralLoss { loss: String }, + #[error("Immonium ions must be a bare uppercase residue, got '{annotation}'")] + UnsupportedModifiedImmonium { annotation: String }, + #[error("Ran out of distinct unknown-ion labels: the 8-bit ordinal is exhausted")] + UnknownIonsExhausted, + #[error("Could not parse '{annotation}': {context}")] + ParsingError { + annotation: String, + context: &'static str, + }, +} + +impl IonParsingError { + /// Create a parsing error with context. + pub(crate) fn parse(annotation: &str, context: &'static str) -> Self { + Self::ParsingError { + annotation: annotation.to_string(), + context, + } + } +} diff --git a/rust/micromzpaf/src/lib.rs b/rust/micromzpaf/src/lib.rs index 551e4b24..50b9af5a 100644 --- a/rust/micromzpaf/src/lib.rs +++ b/rust/micromzpaf/src/lib.rs @@ -1,72 +1,180 @@ //! Compact representation of fragment ion annotations for mass spectrometry. //! -//! This crate provides a minimal, memory-efficient encoding of fragment ion annotations -//! that pack ion series, charge state, ordinal position, and isotope offset into a -//! compact structure suitable for high-throughput proteomics workflows. +//! A packed `u32` representation of the fragment annotations used by spectral +//! libraries, including neutral losses, internal fragments and immonium ions. //! -//! # mzPAF Format Compliance +//! # Bit layout //! -//! This implementation supports a subset of the mzPAF annotation format: -//! - Common ion series (a, b, c, d, v, w, x, y, z, precursor) -//! - Charge states (^N notation) -//! - Positive isotope offsets (+Ni notation) +//! ```text +//! bit: 31 30 29 18 17 12 11 8 7 4 3 0 +//! ┌───────┬────────────────┬──────────┬────────┬────────┬───────┐ +//! │ spare │ payload │ loss │isotope │ charge │ kind │ +//! │ 2b │ 12b │ 6b │ 4b │ 4b zz │ 4b │ +//! └───────┴────────────────┴──────────┴────────┴────────┴───────┘ +//! ``` +//! +//! The widths and shifts are checked by `const` assertions. +//! +//! `payload` is reinterpreted per `kind` -- a tagged union inside the word: +//! +//! | kind | payload | +//! |---|---| +//! | backbone (a/b/c/d/v/w/x/y/z), unknown | ordinal, 8b (0..=255) | +//! | internal | start 6b │ end 6b (peptides to 63 residues) | +//! | immonium | residue index 5b | +//! | precursor | unused | +//! +//! `unknown` is discriminant 0 and `charge` is stored biased by one, so the +//! all-zero word is the valid annotation `?0` at charge 1. +//! +//! `charge` is zigzag-encoded to stay signed in 4 bits; `isotope` is unsigned. +//! Constructors validate both before packing. +//! +//! # The mzPAF subset //! -//! NOTABLY boes not support: -//! - Negative isotope offsets (not yet implemented) -//! - Complex neutral losses or modifications +//! Parses: the a/b/c/d/v/w/x/y/z series, precursor (`p`), unknown (`?`), +//! internal fragments (`m:`), bare immonium (`IA`), charge (`^N`), +//! positive isotopes (`+Ni`), neutral losses from [`NeutralLoss`], and the +//! mass-error suffix (`/-0.0003`, `/1.2ppm`). +//! +//! Rejects, rather than coercing onto a nearby representable ion: negative +//! isotope offsets, modified immonium (`IC[Carbamidomethyl]` carries an +//! arbitrary mod string), and losses outside the [`NeutralLoss`] table. //! //! # Examples //! //! ``` -//! use micromzpaf::IonAnnot; -//! -//! // Parse a simple b-ion annotation -//! let ion: IonAnnot = "b12".try_into().unwrap(); -//! assert_eq!(format!("{}", ion), "b12"); +//! use micromzpaf::{IonAnnot, NeutralLoss, split_mass_error}; //! -//! // Parse with charge and isotope +//! // A backbone ion, with charge and isotope suffixes. //! let ion: IonAnnot = "b12+i^3".try_into().unwrap(); -//! assert_eq!(ion.get_charge(), 3); +//! assert_eq!(ion.try_get_ordinal(), Some(12)); +//! assert_eq!((ion.get_charge(), ion.get_isotope()), (3, 1)); +//! +//! // A neutral loss. Two spellings of one chemical loss are one annotation, +//! // and render canonically. +//! let ion: IonAnnot = "y5-CH3SOH".try_into().unwrap(); +//! assert_eq!(ion.loss(), NeutralLoss::Methanesulfenic); +//! assert_eq!(ion.to_string(), "y5-CH4OS"); +//! +//! // An internal fragment spanning residues 2..=11, and a bare immonium ion. +//! // Neither sits on a ladder, so neither has an ordinal. +//! assert_eq!(IonAnnot::try_from("m2:11").unwrap().try_get_ordinal(), None); +//! assert_eq!(IonAnnot::try_from("IA").unwrap().to_string(), "IA"); +//! +//! // mzSpecLib peaks carry observed m/z; the suffix recovers theoretical. +//! let (ion, err) = split_mass_error("y1/-0.0005").unwrap(); +//! assert_eq!(ion, "y1"); +//! assert!((err.unwrap().theoretical_from_observed(175.1184) - 175.1189).abs() < 1e-4); +//! +//! // Unrepresentable annotations fail rather than losing the detail. +//! assert!(IonAnnot::try_from("y1-HCOOH").is_err()); +//! assert!(IonAnnot::try_from("IC[Carbamidomethyl]").is_err()); //! ``` -use rustyms::fragment::FragmentType; +pub mod error; +pub mod loss; +pub mod parse; +pub mod series; + +pub use error::IonParsingError; +pub use loss::NeutralLoss; +pub use parse::{ + MassError, + split_mass_error, +}; +pub use series::{ + INTERNAL_POS_MAX, + IonSeriesOrdinal, + Series, +}; + use serde::{ Deserialize, Serialize, }; +use series::{ + KIND_BITS, + PAYLOAD_BITS, +}; use std::fmt::Display; use std::hash::Hash; -use std::str::FromStr; -use thiserror::Error; + +const KIND_SHIFT: u32 = 0; +const CHARGE_SHIFT: u32 = 4; +const CHARGE_BITS: u32 = 4; +const ISOTOPE_SHIFT: u32 = 8; +const ISOTOPE_BITS: u32 = 4; +const LOSS_SHIFT: u32 = 12; +const LOSS_BITS: u32 = 6; +const PAYLOAD_SHIFT: u32 = 18; + +/// Charge bounds representable by the packed field. +pub const CHARGE_MIN: i8 = -7; +pub const CHARGE_MAX: i8 = 8; +/// Isotope offset bounds representable by the packed field. +pub const ISOTOPE_MIN: i8 = 0; +pub const ISOTOPE_MAX: i8 = mask(ISOTOPE_BITS) as i8; + +const _: () = assert!( + KIND_BITS + CHARGE_BITS + ISOTOPE_BITS + LOSS_BITS + PAYLOAD_BITS <= 32, + "the fields overflow the word" +); +const _: () = assert!(KIND_SHIFT == 0); +const _: () = assert!(CHARGE_SHIFT == KIND_SHIFT + KIND_BITS, "fields must abut"); +const _: () = assert!( + ISOTOPE_SHIFT == CHARGE_SHIFT + CHARGE_BITS, + "fields must abut" +); +const _: () = assert!( + LOSS_SHIFT == ISOTOPE_SHIFT + ISOTOPE_BITS, + "fields must abut" +); +const _: () = assert!(PAYLOAD_SHIFT == LOSS_SHIFT + LOSS_BITS, "fields must abut"); +const _: () = assert!( + zigzag_charge(CHARGE_MIN) <= mask(CHARGE_BITS) + && zigzag_charge(CHARGE_MAX) <= mask(CHARGE_BITS), + "the charge range does not survive zigzag inside CHARGE_BITS" +); +const _: () = assert!( + ISOTOPE_MIN >= 0 && ISOTOPE_MAX as u32 <= mask(ISOTOPE_BITS), + "the isotope range does not fit ISOTOPE_BITS" +); +const _: () = assert!( + NeutralLoss::COUNT as u32 <= mask(LOSS_BITS), + "the loss table has outgrown LOSS_BITS" +); + +#[inline] +pub(crate) const fn mask(bits: u32) -> u32 { + (1u32 << bits) - 1 +} + +/// Map a small signed value to an unsigned value without losing its sign. +#[inline] +const fn zigzag(v: i8) -> u32 { + (((v as i32) << 1) ^ ((v as i32) >> 31)) as u32 +} +#[inline] +const fn unzigzag(u: u32) -> i8 { + (((u >> 1) as i32) ^ -((u & 1) as i32)) as i8 +} + +/// Encode charge with a one-unit bias so zero decodes to charge 1. +#[inline] +const fn zigzag_charge(charge: i8) -> u32 { + zigzag(charge - 1) +} +#[inline] +const fn unzigzag_charge(u: u32) -> i8 { + unzigzag(u) + 1 +} /// Compact representation of fragment annotations. /// -/// This is a very compressed representation of a fragment -/// ion annotation. Essentially we are packing in 32 bytes -/// the ion series (b, y, ...), charge (+1 / -1 ...), -/// ordinal (12 in the ion series) and isotope. -/// -/// It is not meant to represent all possible ions but rather have -/// a very compact representation of the common ones. -/// -/// # Invariants -/// -/// - **charge**: Must be non-zero (±1 to ±127). Zero charge is invalid and rejected by constructors. -/// - **ordinal**: Limited to u8 range (1-255). Peptides with >255 residues cannot be represented. -/// - **isotope**: Isotope offset relative to monoisotopic peak (M+0), range -128 to +127. -/// -/// # Memory Layout -/// -/// The struct fits in 32 bytes with optimal packing: -/// - `IonSeriesOrdinal`: 2 bytes (enum discriminant + u8 ordinal) -/// - `charge`: 1 byte (i8) -/// - `isotope`: 1 byte (i8) -#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] -pub struct IonAnnot { - series_ordinal: IonSeriesOrdinal, - charge: i8, - isotope: i8, -} +/// A packed `u32`; see the crate docs for the bit layout. +#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, Default)] +pub struct IonAnnot(u32); impl Serialize for IonAnnot { fn serialize(&self, serializer: S) -> Result @@ -77,12 +185,10 @@ impl Serialize for IonAnnot { } } -/// Deserializes simple annotations for fragments. +/// Deserializes an annotation from its mzPAF string. /// -/// b12^3 -> b12 charge 3 (implicit 0 isotope) -/// b12+i^3 -> b12 charge 3 isotope 1 -/// b12+3i^3 -> b12 charge 3 isotope 2 -/// b13 -> b13 (implicit charge 1 and isotope 0) +/// The wire format is the string, not the packed word, so the bit layout can +/// change without breaking existing files. impl<'de> Deserialize<'de> for IonAnnot { fn deserialize(deserializer: D) -> Result where @@ -94,458 +200,383 @@ impl<'de> Deserialize<'de> for IonAnnot { } impl IonAnnot { - pub fn try_new( - ion_type: char, - ordinal: Option, + /// Build an annotation from its parts, validating the packed fields. + pub fn new( + series: IonSeriesOrdinal, + loss: NeutralLoss, charge: i8, isotope: i8, ) -> Result { - Ok(Self { - series_ordinal: IonSeriesOrdinal::try_new(ion_type, ordinal)?, - charge, - isotope, - }) + if charge == 0 { + return Err(IonParsingError::ChargeCannotBeZero); + } + if !(CHARGE_MIN..=CHARGE_MAX).contains(&charge) { + return Err(IonParsingError::ChargeOutOfRange { charge }); + } + if !(ISOTOPE_MIN..=ISOTOPE_MAX).contains(&isotope) { + return Err(IonParsingError::IsotopeOutOfRange { isotope }); + } + let (kind, payload) = series.to_parts(); + Ok(IonAnnot( + (kind << KIND_SHIFT) + | ((zigzag_charge(charge) & mask(CHARGE_BITS)) << CHARGE_SHIFT) + | (((isotope as u32) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT) + | ((loss as u32 & mask(LOSS_BITS)) << LOSS_SHIFT) + | ((payload & mask(PAYLOAD_BITS)) << PAYLOAD_SHIFT), + )) } - pub fn from_fragment( - frag: FragmentType, + /// Build a backbone, precursor or unknown annotation from its series letter. + pub fn try_new( + ion_type: char, + ordinal: Option, charge: i8, isotope: i8, ) -> Result { - Ok(Self { - series_ordinal: IonSeriesOrdinal::try_from(frag)?, + Self::new( + IonSeriesOrdinal::from_series_char(ion_type, ordinal)?, + NeutralLoss::None, charge, isotope, - }) + ) } - pub fn terminality(&self) -> IonSeriesTerminality { - self.series_ordinal.terminality() + #[inline] + fn payload(self) -> u32 { + (self.0 >> PAYLOAD_SHIFT) & mask(PAYLOAD_BITS) } - pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { - let new_isotope = - self.isotope - .checked_add(offset_neutrons) - .ok_or_else(|| IonParsingError::Custom { - error: format!( - "Isotope offset overflow: {} + {} exceeds i8 range", - self.isotope, offset_neutrons - ), - })?; - - Ok(Self { - series_ordinal: self.series_ordinal, - charge: self.charge, - isotope: new_isotope, - }) + #[inline] + pub fn get_charge(&self) -> i8 { + unzigzag_charge((self.0 >> CHARGE_SHIFT) & mask(CHARGE_BITS)) } - pub fn get_charge(&self) -> i8 { - self.charge + #[inline] + pub fn get_isotope(&self) -> i8 { + ((self.0 >> ISOTOPE_SHIFT) & mask(ISOTOPE_BITS)) as i8 + } + + /// The neutral loss this ion carries, [`NeutralLoss::None`] if it carries + /// none. + #[inline] + pub fn loss(&self) -> NeutralLoss { + NeutralLoss::from_discriminant(((self.0 >> LOSS_SHIFT) & mask(LOSS_BITS)) as u8) } + /// Shift the isotope by `offset_neutrons`. + /// + /// Errors when the result leaves [`ISOTOPE_MIN`]..=[`ISOTOPE_MAX`], naming + /// that range, since the caller cannot see the field width. + pub fn try_with_offset_neutrons(&self, offset_neutrons: i8) -> Result { + let new_isotope = self.get_isotope().saturating_add(offset_neutrons); + if !(ISOTOPE_MIN..=ISOTOPE_MAX).contains(&new_isotope) { + return Err(IonParsingError::IsotopeOutOfRange { + isotope: new_isotope, + }); + } + Ok(IonAnnot( + (self.0 & !(mask(ISOTOPE_BITS) << ISOTOPE_SHIFT)) + | (((new_isotope as u32) & mask(ISOTOPE_BITS)) << ISOTOPE_SHIFT), + )) + } + + /// The series ordinal, for the kinds that have one. + /// + /// `None` for precursor, unknown, internal and immonium: `?1` is a + /// uniqueness counter, not a position in a ladder. pub fn try_get_ordinal(&self) -> Option { - self.series_ordinal.try_get_ordinal() + use IonSeriesOrdinal as S; + match self.series_ordinal() { + S::backbone { ordinal, .. } => Some(ordinal), + S::unknown { .. } | S::precursor | S::internal { .. } | S::immonium { .. } => None, + } + } + + /// The backbone series this annotation belongs to, if any. + /// + /// `None` for precursor, unknown, internal and immonium ions, none of + /// which sit on a backbone ladder. + pub fn try_get_series(&self) -> Option { + match self.series_ordinal() { + IonSeriesOrdinal::backbone { series, .. } => Some(series), + _ => None, + } + } + + /// What kind of ion this is, and whatever that kind carries. + pub fn series_ordinal(&self) -> IonSeriesOrdinal { + IonSeriesOrdinal::from_parts((self.0 >> KIND_SHIFT) & mask(KIND_BITS), self.payload()) } } +/// Guards the one-shot warning in [`IonAnnot::try_from`]. +static MASS_ERROR_DISCARDED: std::sync::Once = std::sync::Once::new(); + impl TryFrom<&str> for IonAnnot { type Error = IonParsingError; + /// Parse an annotation, ignoring any mass-error suffix. + /// + /// Use [`split_mass_error`] when the suffix is needed for m/z correction. fn try_from(value: &str) -> Result { - let (rest, charge) = match value.split_once('^') { - Some((rest, charge)) => { - let charge = charge - .parse::() - .map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse the charge number"), - })?; - (rest, charge) - } - None => (value, 1), - }; - - // TODO: Implement 'negative isotopes' parsing ... right now I dont use them - // for serialization ... - // Note that this is not 100% compliant with mzPAF - let (rest, isotope) = match rest.split_once('+') { - Some((rest, adducts)) => { - // Make sure the adduct is only +{number}?i - let adducts = adducts - .strip_suffix('i') - .ok_or(IonParsingError::ParsingError { - error: adducts.to_string(), - context: Some("Unsupported adduct found"), - })?; - // If its empty its an implicit 1 isotope - // Since we stripped the 'i' from '+i' - let isotope = if adducts.is_empty() { - 1 - } else { - adducts - .parse::() - .map_err(|_| IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Unable to parse the isotope number"), - })? - }; - (rest, isotope) - } - None => (rest, 0), - }; - let series_ord = IonSeriesOrdinal::from_str(rest)?; - if charge == 0 { - return Err(IonParsingError::ParsingError { - error: value.to_string(), - context: Some("Charge cannot be 0"), + let (ion, mass_error) = split_mass_error(value)?; + if mass_error.is_some() { + MASS_ERROR_DISCARDED.call_once(|| { + tracing::warn!( + "annotation {value:?} carries a mass-error suffix, which \ + `IonAnnot` does not store; the m/z used is whatever the \ + caller supplied. Use `split_mass_error` to recover the \ + theoretical m/z. Warning once per process." + ); }); } - Ok(Self { - series_ordinal: series_ord, - charge, - isotope, - }) + let (series, suffixes) = parse::parse_ion(ion)?; + Self::new(series, suffixes.loss, suffixes.charge, suffixes.isotope) } } impl Display for IonAnnot { + /// Renders the canonical mzPAF spelling, the exact inverse of parsing except + /// that a non-canonical loss spelling (`-CH3SOH`) renders canonically + /// (`-CH4OS`). The mass-error suffix is not part of the annotation. fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "{}", self.series_ordinal)?; + write!(f, "{}", self.series_ordinal())?; + write!(f, "{}", self.loss())?; - match self.isotope { + match self.get_isotope() { 0 => {} 1 => write!(f, "+i")?, i => write!(f, "+{}i", i)?, } - if self.charge != 1 { - write!(f, "^{}", self.charge)?; + let charge = self.get_charge(); + if charge != 1 { + write!(f, "^{}", charge)?; } Ok(()) } } -#[derive(Debug, Error)] -pub enum IonParsingError { - #[error("Ordinal {ordinal} out of range for series '{series:?}'")] - OrdinalOutOfRange { ordinal: i32, series: Option }, - #[error("Unsupported fragment type: '{fragment_type}'")] - UnsupportedFragmentType { fragment_type: char }, - #[error("Parsing error: {error}{}", .context.map(|c| format!(" ({})", c)).unwrap_or_default())] - ParsingError { - error: String, - context: Option<&'static str>, - }, - #[error("{error}")] - Custom { error: String }, +/// Generates unique unknown-ion labels within a precursor. +#[derive(Debug, Default)] +pub struct UnknownIonCounter(u8); + +impl UnknownIonCounter { + /// The next unused unknown label at `charge`. + pub fn next_unknown(&mut self, charge: i8) -> Result { + let ordinal = self + .0 + .checked_add(1) + .ok_or(IonParsingError::UnknownIonsExhausted)?; + let annot = IonAnnot::try_new('?', Some(ordinal), charge, 0)?; + self.0 = ordinal; + Ok(annot) + } } -/// Refers to what terminus of the original peptide retains the -/// charge after a fragmentation. -#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] -pub enum IonSeriesTerminality { - NTerm, - CTerm, - None, -} +#[cfg(test)] +mod tests { + use super::*; -#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy, Default)] -#[allow(non_camel_case_types)] -pub enum IonSeriesOrdinal { - a { - ordinal: u8, - }, - b { - ordinal: u8, - }, - c { - ordinal: u8, - }, - d { - ordinal: u8, - }, - v { - ordinal: u8, - }, - w { - ordinal: u8, - }, - x { - ordinal: u8, - }, - y { - ordinal: u8, - }, - z { - ordinal: u8, - }, - unknown { - ordinal: u8, - }, - precursor, - - /// This variant should not be used directly ... its mainly added to satisfy trait constraints by TinyVec - #[default] - None, -} + fn ion(s: &str) -> IonAnnot { + IonAnnot::try_from(s).unwrap_or_else(|e| panic!("{s:?} must parse: {e}")) + } -impl IonSeriesOrdinal { - pub fn try_new(series: char, ordinal: Option) -> Result { - let tmp = match (series, ordinal) { - ('a', Some(ordinal)) => Self::a { ordinal }, - ('b', Some(ordinal)) => Self::b { ordinal }, - ('c', Some(ordinal)) => Self::c { ordinal }, - ('d', Some(ordinal)) => Self::d { ordinal }, - ('v', Some(ordinal)) => Self::v { ordinal }, - ('w', Some(ordinal)) => Self::w { ordinal }, - ('x', Some(ordinal)) => Self::x { ordinal }, - ('y', Some(ordinal)) => Self::y { ordinal }, - ('z', Some(ordinal)) => Self::z { ordinal }, - ('?', Some(ordinal)) => Self::unknown { ordinal }, - ('p', None) => Self::precursor, - ('p', Some(ordinal)) => { - return Err(IonParsingError::OrdinalOutOfRange { - ordinal: ordinal as i32, - series: Some(series), - }); - } - _ => { - return Err(IonParsingError::UnsupportedFragmentType { - fragment_type: series, - }); - } - }; + /// Keep the annotation and annotation/intensity pair compact. + #[test] + fn the_added_fields_cost_no_space() { + assert_eq!(size_of::(), 4); + assert_eq!(size_of::<(IonAnnot, f32)>(), 8); + + let loaded = IonAnnot::new( + IonSeriesOrdinal::try_internal(2, 11).expect("in range"), + NeutralLoss::PhosphoricAcidWater, + 3, + 2, + ) + .expect("every field at once"); + assert_eq!(size_of_val(&loaded), 4); + assert_eq!(loaded.loss(), NeutralLoss::PhosphoricAcidWater); + assert_eq!( + loaded.series_ordinal(), + IonSeriesOrdinal::internal { start: 2, end: 11 } + ); + assert_eq!((loaded.get_charge(), loaded.get_isotope()), (3, 2)); + } - Ok(tmp) + /// The default packed word is a valid annotation and round-trips via serde. + #[test] + fn the_default_annotation_is_a_real_annotation() { + let d = IonAnnot::default(); + assert_eq!(d.to_string(), "?0"); + assert_eq!(d.get_charge(), 1); + assert_eq!(serde_json::to_string(&d).unwrap(), "\"?0\""); + assert_eq!(IonAnnot::try_from("?0").expect("the default parses"), d); } - pub fn terminality(&self) -> IonSeriesTerminality { - match self { - IonSeriesOrdinal::a { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::b { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::c { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::d { ordinal: _ } => IonSeriesTerminality::NTerm, - IonSeriesOrdinal::v { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::w { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::x { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::y { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::z { ordinal: _ } => IonSeriesTerminality::CTerm, - IonSeriesOrdinal::unknown { ordinal: _ } => IonSeriesTerminality::None, - IonSeriesOrdinal::precursor => IonSeriesTerminality::None, - IonSeriesOrdinal::None => panic!("IonSeriesOrdinal::None should not be used directly"), + #[test] + fn test_deserialize() { + let cases = [ + ("b12", 'b', Some(12u8), 1i8, 0i8), + ("b12^3", 'b', Some(12), 3, 0), + ("y12^3", 'y', Some(12), 3, 0), + ("b12+i^3", 'b', Some(12), 3, 1), + ("b12+3i^3", 'b', Some(12), 3, 3), + ("b13", 'b', Some(13), 1, 0), + ("p^2", 'p', None, 2, 0), + ("p", 'p', None, 1, 0), + ("?12^2", '?', Some(12), 2, 0), + ]; + + for (input, series, ordinal, charge, isotope) in cases { + let annot = ion(input); + let expected = IonAnnot::try_new(series, ordinal, charge, isotope).unwrap(); + assert_eq!(annot, expected, "{input}"); + assert_eq!(annot.get_charge(), charge, "{input}"); + assert_eq!(annot.get_isotope(), isotope, "{input}"); + assert_eq!(format!("{}", annot), input); } } - pub fn try_get_ordinal(&self) -> Option { - match self { - IonSeriesOrdinal::a { ordinal } => Some(*ordinal), - IonSeriesOrdinal::b { ordinal } => Some(*ordinal), - IonSeriesOrdinal::c { ordinal } => Some(*ordinal), - IonSeriesOrdinal::d { ordinal } => Some(*ordinal), - IonSeriesOrdinal::v { ordinal } => Some(*ordinal), - IonSeriesOrdinal::w { ordinal } => Some(*ordinal), - IonSeriesOrdinal::x { ordinal } => Some(*ordinal), - IonSeriesOrdinal::y { ordinal } => Some(*ordinal), - IonSeriesOrdinal::z { ordinal } => Some(*ordinal), - IonSeriesOrdinal::unknown { .. } => None, - // ?1 does not mean its an ordinal, just a placeholder - IonSeriesOrdinal::precursor => None, - IonSeriesOrdinal::None => None, + /// Every field must survive the pack/unpack at its extremes. + #[test] + fn every_field_round_trips_at_its_extremes() { + for charge in CHARGE_MIN..=CHARGE_MAX { + if charge == 0 { + continue; + } + for isotope in ISOTOPE_MIN..=ISOTOPE_MAX { + for ordinal in [1u8, 2, 127, 255] { + for loss in [ + NeutralLoss::None, + NeutralLoss::Water, + NeutralLoss::PhosphoricAcidWater, + ] { + let series = IonSeriesOrdinal::backbone { + series: Series::y, + ordinal, + }; + let a = IonAnnot::new(series, loss, charge, isotope).expect("in range"); + assert_eq!(a.get_charge(), charge); + assert_eq!(a.get_isotope(), isotope); + assert_eq!(a.try_get_ordinal(), Some(ordinal)); + assert_eq!(a.loss(), loss); + } + } + } } } -} -impl Display for IonSeriesOrdinal { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - match self { - IonSeriesOrdinal::a { ordinal } => write!(f, "a{}", ordinal), - IonSeriesOrdinal::b { ordinal } => write!(f, "b{}", ordinal), - IonSeriesOrdinal::c { ordinal } => write!(f, "c{}", ordinal), - IonSeriesOrdinal::d { ordinal } => write!(f, "d{}", ordinal), - IonSeriesOrdinal::v { ordinal } => write!(f, "v{}", ordinal), - IonSeriesOrdinal::w { ordinal } => write!(f, "w{}", ordinal), - IonSeriesOrdinal::x { ordinal } => write!(f, "x{}", ordinal), - IonSeriesOrdinal::y { ordinal } => write!(f, "y{}", ordinal), - IonSeriesOrdinal::z { ordinal } => write!(f, "z{}", ordinal), - IonSeriesOrdinal::unknown { ordinal } => write!(f, "?{}", ordinal), - IonSeriesOrdinal::precursor => write!(f, "p"), - IonSeriesOrdinal::None => panic!("IonSeriesOrdinal::None should not be used directly"), - } + #[test] + fn out_of_range_is_rejected_not_truncated() { + assert!(matches!( + IonAnnot::try_new('y', Some(1), CHARGE_MAX + 1, 0), + Err(IonParsingError::ChargeOutOfRange { .. }) + )); + assert!(matches!( + IonAnnot::try_new('y', Some(1), 1, ISOTOPE_MAX + 1), + Err(IonParsingError::IsotopeOutOfRange { .. }) + )); + assert!(IonAnnot::try_new('y', Some(1), 0, 0).is_err()); + assert!(matches!( + IonAnnot::try_from("m64:1"), + Err(IonParsingError::OrdinalOutOfRange { .. }) + )); } -} -impl FromStr for IonSeriesOrdinal { - type Err = IonParsingError; - - fn from_str(s: &str) -> Result { - // "b12" split into "b" and "12" - let (series_chunk, ordinal_chunk) = s.split_at(1); - let series_id = series_chunk.chars().next(); - let series_ordinal = ordinal_chunk.parse::(); - - match (series_id, series_ordinal) { - (None, _) => Err(IonParsingError::ParsingError { - error: s.to_string(), - context: Some("Empty string"), - }), - (Some(x), Ok(y)) => IonSeriesOrdinal::try_new(x, Some(y)), - (Some('p'), Err(_)) => Ok(IonSeriesOrdinal::precursor), - (Some(_), Err(err)) => Err(IonParsingError::ParsingError { - error: format!("{ordinal_chunk} -> {err:?}"), - context: Some("Unable to parse the ordinal number"), - }), + #[test] + fn negative_isotopes_are_unrepresentable_not_just_unparsed() { + assert!(matches!( + IonAnnot::try_new('y', Some(5), 1, -1), + Err(IonParsingError::IsotopeOutOfRange { isotope: -1 }) + )); + assert!(ion("y5").try_with_offset_neutrons(-1).is_err()); + assert!(IonAnnot::try_from("y5-2i").is_err()); + assert!(IonAnnot::try_from("y5+-2i").is_err()); + for isotope in ISOTOPE_MIN..=ISOTOPE_MAX { + let a = IonAnnot::try_new('y', Some(5), 1, isotope).expect("in range"); + let text = a.to_string(); + assert!(!text.contains("+-"), "{text} is not mzPAF"); + assert_eq!(ion(&text), a, "{text}"); } } -} - -impl TryFrom for IonSeriesOrdinal { - type Error = IonParsingError; - fn try_from(value: FragmentType) -> Result { - fn try_convert_ordinal(ordinal: usize, series: char) -> Result { - ordinal - .try_into() - .map_err(|_| IonParsingError::OrdinalOutOfRange { - ordinal: ordinal as i32, - series: Some(series), - }) - } - let tmp = match value { - FragmentType::a(ordinal, _) => IonSeriesOrdinal::a { - ordinal: try_convert_ordinal(ordinal.series_number, 'a')?, - }, - FragmentType::b(ordinal, _) => IonSeriesOrdinal::b { - ordinal: try_convert_ordinal(ordinal.series_number, 'b')?, - }, - FragmentType::c(ordinal, _) => IonSeriesOrdinal::c { - ordinal: try_convert_ordinal(ordinal.series_number, 'c')?, - }, - FragmentType::d(ordinal, _, _, _, _) => IonSeriesOrdinal::d { - ordinal: try_convert_ordinal(ordinal.series_number, 'd')?, - }, - FragmentType::v(ordinal, _, _, _) => IonSeriesOrdinal::v { - ordinal: try_convert_ordinal(ordinal.series_number, 'v')?, - }, - FragmentType::w(ordinal, _, _, _, _) => IonSeriesOrdinal::w { - ordinal: try_convert_ordinal(ordinal.series_number, 'w')?, - }, - FragmentType::x(ordinal, _) => IonSeriesOrdinal::x { - ordinal: try_convert_ordinal(ordinal.series_number, 'x')?, - }, - FragmentType::y(ordinal, _) => IonSeriesOrdinal::y { - ordinal: try_convert_ordinal(ordinal.series_number, 'y')?, - }, - FragmentType::z(ordinal, _) => IonSeriesOrdinal::z { - ordinal: try_convert_ordinal(ordinal.series_number, 'z')?, - }, - FragmentType::Precursor => IonSeriesOrdinal::precursor, - _ => { - return Err(IonParsingError::Custom { - error: format!("Unsupported fragment type: {value:?}"), - }); - } - }; - Ok(tmp) + #[test] + fn isotope_offset_respects_the_field_bound() { + let a = ion("y5"); + assert_eq!(a.try_with_offset_neutrons(2).unwrap().get_isotope(), 2); + assert!(a.try_with_offset_neutrons(ISOTOPE_MAX + 1).is_err()); } -} -#[cfg(test)] -mod tests { - use super::*; + #[test] + fn parses_neutral_losses() { + let a = ion("y5-H2O"); + assert_eq!(a.loss(), NeutralLoss::Water); + assert_eq!(a.try_get_ordinal(), Some(5)); + assert_eq!(format!("{}", a), "y5-H2O"); + + assert_eq!(ion("y5-CH3SOH"), ion("y5-CH4OS")); + assert_eq!(format!("{}", ion("y5-CH3SOH")), "y5-CH4OS"); + + let b = ion("y10-NH3+i^2"); + assert_eq!(b.loss(), NeutralLoss::Ammonia); + assert_eq!(b.get_isotope(), 1); + assert_eq!(b.get_charge(), 2); + } #[test] - fn test_ion_series_ord_from_str() { - let ion: IonSeriesOrdinal = IonSeriesOrdinal::from_str("b12").unwrap(); - assert_eq!(ion, IonSeriesOrdinal::b { ordinal: 12 }); + fn parses_internal_fragments() { + let a = ion("m2:11"); + assert_eq!( + a.series_ordinal(), + IonSeriesOrdinal::internal { start: 2, end: 11 } + ); + assert_eq!(a.try_get_ordinal(), None); + assert_eq!(format!("{}", a), "m2:11"); + + let b = ion("m11:12-CO"); + assert_eq!(b.loss(), NeutralLoss::CarbonMonoxide); + assert_eq!(format!("{}", b), "m11:12-CO"); } #[test] - fn test_deserialize() { - let serde_pairs = vec![ - ( - "b12", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 1, - isotope: 0, - }, - ), - ( - "b12^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 0, - }, - ), - ( - "y12^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::y { ordinal: 12 }, - charge: 3, - isotope: 0, - }, - ), - ( - "b12+i^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 1, - }, - ), - ( - "b12+3i^3", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 12 }, - charge: 3, - isotope: 3, - }, - ), - ( - "b13", - IonAnnot { - series_ordinal: IonSeriesOrdinal::b { ordinal: 13 }, - charge: 1, - isotope: 0, - }, - ), - ( - "p^2", - IonAnnot { - series_ordinal: IonSeriesOrdinal::precursor, - charge: 2, - isotope: 0, - }, - ), - ( - "p", - IonAnnot { - series_ordinal: IonSeriesOrdinal::precursor, - charge: 1, - isotope: 0, - }, - ), - ( - "?12^2", - IonAnnot { - series_ordinal: IonSeriesOrdinal::unknown { ordinal: 12 }, - charge: 2, - isotope: 0, - }, - ), - ]; + fn parses_bare_immonium_and_rejects_modified() { + let a = ion("IA"); + assert_eq!( + a.series_ordinal(), + IonSeriesOrdinal::immonium { residue: 'A' } + ); + assert_eq!(format!("{}", a), "IA"); + + assert!(matches!( + IonAnnot::try_from("IC[Carbamidomethyl]"), + Err(IonParsingError::UnsupportedModifiedImmonium { .. }) + )); + } - for (input, expected) in serde_pairs { - let annot = IonAnnot::try_from(input).unwrap(); - assert_eq!(annot, expected); + #[test] + fn parses_and_applies_the_mass_error_suffix() { + let (rest, err) = split_mass_error("y1/-0.0005").unwrap(); + assert_eq!(ion(rest), ion("y1")); + assert_eq!(err, Some(MassError::Da(-0.0005))); + let theo = err.unwrap().theoretical_from_observed(175.1184); + assert!((theo - 175.1189).abs() < 1e-9, "got {theo}"); + + let (_, err) = split_mass_error("y6/1.2ppm").unwrap(); + assert_eq!(err, Some(MassError::Ppm(1.2))); + let theo = err.unwrap().theoretical_from_observed(700.0); + assert!((theo - 699.99916).abs() < 1e-4, "got {theo}"); + + assert_eq!(split_mass_error("y6").unwrap(), ("y6", None)); + assert_eq!(IonAnnot::try_from("y1/-0.0005").unwrap(), ion("y1")); + } - // Re-serialize and check that its the same - let serialized = format!("{}", annot); - assert_eq!(serialized, input); - } + #[test] + fn unrepresentable_loss_is_rejected_not_stripped() { + assert!(matches!( + IonAnnot::try_from("y1-HCOOH"), + Err(IonParsingError::UnsupportedNeutralLoss { .. }) + )); } } diff --git a/rust/micromzpaf/src/loss.rs b/rust/micromzpaf/src/loss.rs new file mode 100644 index 00000000..29289c84 --- /dev/null +++ b/rust/micromzpaf/src/loss.rs @@ -0,0 +1,385 @@ +//! Neutral losses, keyed by atomic composition rather than by spelling. +//! +//! Libraries may write the same chemical loss in different ways: +//! +//! | written | library | composition | +//! |---|---|---| +//! | `-CH3SOH` | NIST | C1H4O1S1 | +//! | `-CH4OS` | SpectraST | C1H4O1S1 | +//! | `-NH2-CO-CH2SH` | NIST | C2H5N1O1S1 | +//! | `-C2H5NOS` | SpectraST | C2H5N1O1S1 | +//! +//! Parsing goes `text -> composition -> discriminant`; the packed annotation +//! stores only the discriminant. + +use std::fmt::Display; + +use crate::IonParsingError; + +/// Atom counts for the elements that appear in peptide neutral losses. +/// +/// These losses use only C/H/N/O/S/P. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub(crate) struct Composition { + c: u8, + h: u8, + n: u8, + o: u8, + s: u8, + p: u8, +} + +/// A [`Composition`] naming only the elements it contains: `C!(h: 2, o: 1)`. +macro_rules! C { + ($($field:ident: $count:expr),+ $(,)?) => { + Composition { $($field: $count,)+ ..Composition::ZERO } + }; +} + +impl Composition { + /// All-zero, for [`TABLE`] rows to fill in only the elements they contain. + const ZERO: Self = Self { + c: 0, + h: 0, + n: 0, + o: 0, + s: 0, + p: 0, + }; + + /// The count for one element symbol, or `None` if this crate does not + /// represent that element. + fn count_mut(&mut self, symbol: u8) -> Option<&mut u8> { + Some(match symbol { + b'C' => &mut self.c, + b'H' => &mut self.h, + b'N' => &mut self.n, + b'O' => &mut self.o, + b'S' => &mut self.s, + b'P' => &mut self.p, + _ => return None, + }) + } + + /// The six counts, for tests that need to range over them. + #[cfg(test)] + fn counts(self) -> [u8; 6] { + [self.c, self.h, self.n, self.o, self.s, self.p] + } + + /// Combine element-wise, saturating each count. + fn zip(self, other: Self, f: impl Fn(u8, u8) -> u8) -> Self { + Self { + c: f(self.c, other.c), + h: f(self.h, other.h), + n: f(self.n, other.n), + o: f(self.o, other.o), + s: f(self.s, other.s), + p: f(self.p, other.p), + } + } + + /// Multiply every count. Used for the `2H2O` multiplier form. + fn scaled(self, k: u8) -> Self { + self.zip(Self::ZERO, |n, _| n.saturating_mul(k)) + } + + fn plus(self, other: Self) -> Self { + self.zip(other, u8::saturating_add) + } + + /// Parse a bare formula like `H2O`, `CH4OS`, `C2H5NOS`. + /// + /// Only single-letter C/H/N/O/S/P elements are recognized. + fn parse_formula(s: &str) -> Result { + if s.is_empty() { + return Err(IonParsingError::parse(s, "Empty neutral-loss formula")); + } + let mut out = Composition::default(); + let b = s.as_bytes(); + let mut i = 0; + while i < b.len() { + let symbol = b[i]; + i += 1; + let start = i; + while i < b.len() && b[i].is_ascii_digit() { + i += 1; + } + let count: u8 = if start == i { + 1 + } else { + s[start..i].parse().map_err(|_| { + IonParsingError::parse(s, "Neutral-loss atom count out of range") + })? + }; + let slot = out + .count_mut(symbol) + .ok_or_else(|| IonParsingError::parse(s, "Unsupported element in neutral loss"))?; + *slot = slot.saturating_add(count); + } + Ok(out) + } + + /// Parse a full loss expression: `-` separated terms, each optionally + /// prefixed by a repeat count. `2H2O`, `H2O-NH3`, `NH2-CO-CH2SH`. + /// + /// Terms are summed, so ordering and multiplier spelling are normalized. + pub(crate) fn parse_expression(s: &str) -> Result { + let mut total = Composition::default(); + for term in s.split('-') { + let term = term.trim(); + if term.is_empty() { + return Err(IonParsingError::parse( + s, + "Empty term in neutral-loss expression", + )); + } + // Leading digits are a repeat count for the whole term. + let digits = term + .find(|c: char| !c.is_ascii_digit()) + .unwrap_or(term.len()); + let (mult, formula) = if digits > 0 { + let m: u8 = term[..digits].parse().map_err(|_| { + IonParsingError::parse(s, "Neutral-loss repeat count out of range") + })?; + (m, &term[digits..]) + } else { + (1, term) + }; + total = total.plus(Composition::parse_formula(formula)?.scaled(mult)); + } + Ok(total) + } +} + +/// The neutral losses this crate can represent, as a packed discriminant. +/// +/// Inputs outside this table are reported as unrepresentable. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)] +#[repr(u8)] +pub enum NeutralLoss { + #[default] + None = 0, + /// H2O + Water = 1, + /// NH3 + Ammonia = 2, + /// CO + CarbonMonoxide = 3, + /// CO2 + CarbonDioxide = 4, + /// 2 H2O + WaterX2 = 5, + /// 2 NH3 + AmmoniaX2 = 6, + /// H2O + NH3 + WaterAmmonia = 7, + /// CH4OS -- methanesulfenic acid, off oxidized Met. Also spelled `CH3SOH`. + Methanesulfenic = 8, + /// C2H5NOS -- also spelled `NH2-CO-CH2SH`. + Carbamidomethylthiol = 9, + /// H3PO4 -- phospho-Ser/Thr. + PhosphoricAcid = 10, + /// HPO3 -- phospho-Tyr, and phospho-Ser/Thr. + Metaphosphoric = 11, + /// H3PO4 + H2O + PhosphoricAcidWater = 12, +} + +/// `(composition, discriminant, canonical spelling)`, indexed by discriminant. +/// Entries are ordered by discriminant; canonical spelling is emitted by +/// [`Display`]. +const TABLE: &[(Composition, NeutralLoss, &str)] = &[ + (C!(h: 2, o: 1), NeutralLoss::Water, "H2O"), + (C!(h: 3, n: 1), NeutralLoss::Ammonia, "NH3"), + (C!(c: 1, o: 1), NeutralLoss::CarbonMonoxide, "CO"), + (C!(c: 1, o: 2), NeutralLoss::CarbonDioxide, "CO2"), + (C!(h: 4, o: 2), NeutralLoss::WaterX2, "2H2O"), + (C!(h: 6, n: 2), NeutralLoss::AmmoniaX2, "2NH3"), + (C!(h: 5, n: 1, o: 1), NeutralLoss::WaterAmmonia, "H2O-NH3"), + ( + C!(c: 1, h: 4, o: 1, s: 1), + NeutralLoss::Methanesulfenic, + "CH4OS", + ), + ( + C!(c: 2, h: 5, n: 1, o: 1, s: 1), + NeutralLoss::Carbamidomethylthiol, + "C2H5NOS", + ), + (C!(h: 3, o: 4, p: 1), NeutralLoss::PhosphoricAcid, "H3PO4"), + (C!(h: 1, o: 3, p: 1), NeutralLoss::Metaphosphoric, "HPO3"), + ( + C!(h: 5, o: 5, p: 1), + NeutralLoss::PhosphoricAcidWater, + "H3PO4-H2O", + ), +]; + +impl NeutralLoss { + /// Number of loss entries in the table. + pub(crate) const COUNT: u8 = TABLE.len() as u8; + + /// Decode a packed loss discriminant. + pub(crate) fn from_discriminant(d: u8) -> Self { + Self::at(d).map_or(Self::None, |(_, loss, _)| *loss) + } + + /// The [`TABLE`] row a discriminant names, or `None` for [`Self::None`] and + /// for the reserved values above the table. + fn at(d: u8) -> Option<&'static (Composition, NeutralLoss, &'static str)> { + TABLE.get((d as usize).checked_sub(1)?) + } + + /// Resolve a loss expression (without the leading `-`) to a discriminant. + /// + /// `None` means the composition is valid but not in the supported table. + pub(crate) fn from_expression(s: &str) -> Result, IonParsingError> { + let comp = Composition::parse_expression(s)?; + Ok(TABLE + .iter() + .find(|(c, _, _)| *c == comp) + .map(|(_, l, _)| *l)) + } + + /// Canonical spelling, without the leading `-`. Empty for [`Self::None`]. + pub(crate) fn canonical(self) -> &'static str { + Self::at(self as u8).map_or("", |(_, _, spelling)| *spelling) + } +} + +impl Display for NeutralLoss { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + if *self == NeutralLoss::None { + return Ok(()); + } + write!(f, "-{}", self.canonical()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn atom_counts_parse_and_reject_out_of_range() { + assert_eq!( + Composition::parse_expression("H2O").unwrap(), + C!(h: 2, o: 1), + "an element with no digits is one atom" + ); + assert_eq!( + Composition::parse_expression("C10H12").unwrap(), + C!(c: 10, h: 12), + "counts are multi-digit, not one digit per element" + ); + assert_eq!( + Composition::parse_expression("C255").unwrap(), + C!(c: 255), + "the u8 slot is full at 255" + ); + assert!( + Composition::parse_expression("C256").is_err(), + "one past the slot must be an error, not a wrap to 0" + ); + } + + #[test] + fn different_spellings_resolve_to_one_loss() { + assert_eq!( + NeutralLoss::from_expression("CH3SOH").unwrap(), + Some(NeutralLoss::Methanesulfenic) + ); + assert_eq!( + NeutralLoss::from_expression("CH4OS").unwrap(), + Some(NeutralLoss::Methanesulfenic) + ); + assert_eq!( + NeutralLoss::from_expression("NH2-CO-CH2SH").unwrap(), + Some(NeutralLoss::Carbamidomethylthiol) + ); + assert_eq!( + NeutralLoss::from_expression("C2H5NOS").unwrap(), + Some(NeutralLoss::Carbamidomethylthiol) + ); + } + + #[test] + fn ordering_and_multipliers_normalize() { + assert_eq!( + NeutralLoss::from_expression("H2O-NH3").unwrap(), + NeutralLoss::from_expression("NH3-H2O").unwrap() + ); + assert_eq!( + NeutralLoss::from_expression("2H2O").unwrap(), + NeutralLoss::from_expression("H2O-H2O").unwrap() + ); + assert_eq!( + NeutralLoss::from_expression("2H2O").unwrap(), + Some(NeutralLoss::WaterX2) + ); + } + + #[test] + fn unrepresentable_is_distinct_from_malformed() { + assert_eq!(NeutralLoss::from_expression("HCOOH").unwrap(), None); + assert!(NeutralLoss::from_expression("Xe2").is_err()); + assert!(NeutralLoss::from_expression("").is_err()); + assert!(NeutralLoss::from_expression("H2O-").is_err()); + } + + #[test] + fn table_is_indexed_by_discriminant() { + for (i, (_comp, loss, canon)) in TABLE.iter().enumerate() { + assert_eq!( + *loss as usize, + i + 1, + "{canon} sits at offset {i} but its discriminant is {}", + *loss as u8 + ); + } + // Nothing outside the table decodes to a loss, in either direction. + assert_eq!(NeutralLoss::from_discriminant(0), NeutralLoss::None); + assert_eq!( + NeutralLoss::from_discriminant(TABLE.len() as u8 + 1), + NeutralLoss::None + ); + assert_eq!(NeutralLoss::None.canonical(), ""); + } + + #[test] + fn table_round_trips_through_canonical_spelling() { + for (_comp, loss, canon) in TABLE { + assert_eq!( + NeutralLoss::from_expression(canon).unwrap(), + Some(*loss), + "canonical spelling {canon} must resolve to its own loss" + ); + assert_eq!(loss.canonical(), *canon); + assert_eq!( + NeutralLoss::from_discriminant(*loss as u8), + *loss, + "{canon} does not survive the discriminant round trip" + ); + } + } + + #[test] + fn table_compositions_are_unique() { + for (i, (a, _, sa)) in TABLE.iter().enumerate() { + for (b, _, sb) in TABLE.iter().skip(i + 1) { + assert_ne!(a, b, "{sa} and {sb} share a composition"); + } + assert!( + a.counts().iter().all(|&n| n < u8::MAX), + "{sa} holds a saturated atom count" + ); + } + } + + #[test] + fn display_adds_the_prefix_and_none_renders_empty() { + assert_eq!(NeutralLoss::None.to_string(), ""); + assert_eq!(NeutralLoss::Water.to_string(), "-H2O"); + } +} diff --git a/rust/micromzpaf/src/parse.rs b/rust/micromzpaf/src/parse.rs new file mode 100644 index 00000000..875078f2 --- /dev/null +++ b/rust/micromzpaf/src/parse.rs @@ -0,0 +1,144 @@ +//! Peeling an mzPAF annotation apart, and the mass-error suffix. +//! +//! The suffixes come off right to left, and that order is not arbitrary: each +//! delimiter can only appear to the right of the previous one, so peeling in +//! this sequence is unambiguous while any other order is not. +//! +//! ```text +//! y5 - H2O + 2i ^ 3 / -0.0003 +//! └core┘└loss┘└iso┘└ch┘└error┘ +//! 4 3 2 1 0 <- order peeled +//! ``` +//! +//! `0` is [`split_mass_error`], `1..=3` are `peel_suffixes`, and `4` is +//! `IonSeriesOrdinal::parse`, which owns the spelling of the ion itself. + +use crate::series::IonSeriesOrdinal; +use crate::{ + IonParsingError, + NeutralLoss, +}; + +/// A mass-error suffix on an mzPAF annotation: observed minus theoretical. +/// +/// mzSpecLib peak lists carry *observed* m/z, so the theoretical value a +/// library reader wants is `observed - error`. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum MassError { + /// Absolute, in daltons. + Da(f64), + /// Relative, in parts per million. + Ppm(f64), +} + +impl MassError { + /// Recover the theoretical m/z from the observed one. + pub fn theoretical_from_observed(&self, observed: f64) -> f64 { + match self { + MassError::Da(d) => observed - d, + MassError::Ppm(p) => observed / (1.0 + p * 1e-6), + } + } +} + +/// Split the trailing `/[ppm]` off an annotation, if present. +/// +/// Public because a library reader needs the error even when the ion itself is +/// unrepresentable: the error is what recovers theoretical m/z, so a peak that +/// ends up with an unknown label still gets a correct mass. +pub fn split_mass_error(s: &str) -> Result<(&str, Option), IonParsingError> { + let Some((head, tail)) = s.rsplit_once('/') else { + return Ok((s, None)); + }; + let (num, is_ppm) = match tail.strip_suffix("ppm") { + Some(n) => (n, true), + None => (tail, false), + }; + let v: f64 = num + .parse() + .map_err(|_| IonParsingError::parse(s, "mass-error suffix is not a number"))?; + Ok(( + head, + Some(if is_ppm { + MassError::Ppm(v) + } else { + MassError::Da(v) + }), + )) +} + +/// The three suffixes an annotation can carry after its ion kind. +pub(crate) struct Suffixes { + pub(crate) charge: i8, + pub(crate) isotope: i8, + pub(crate) loss: NeutralLoss, +} + +/// Peel `^N`, `+Ni` and `-` off `value`, returning the ion core and the +/// suffixes. `value` must already have had any mass-error suffix removed. +/// +/// Defaults when a suffix is absent: charge 1, isotope 0, no loss. +pub(crate) fn peel_suffixes(value: &str) -> Result<(&str, Suffixes), IonParsingError> { + // Charge: `^N`, rightmost of the three. + let (rest, charge) = match value.split_once('^') { + Some((rest, charge)) => ( + rest, + charge + .parse::() + .map_err(|_| IonParsingError::parse(value, "charge is not a number"))?, + ), + None => (value, 1), + }; + + // Isotope: `+Ni`, or bare `+i` for one. + let (rest, isotope) = match rest.split_once('+') { + Some((rest, adducts)) => { + let digits = adducts.strip_suffix('i').ok_or(IonParsingError::parse( + value, + "only the isotope adduct '+Ni' is supported", + ))?; + let isotope = if digits.is_empty() { + 1 + } else { + digits + .parse::() + .map_err(|_| IonParsingError::parse(value, "isotope count is not a number"))? + }; + (rest, isotope) + } + None => (rest, 0), + }; + + // Loss: everything from the first `-`. Internal fragments spell themselves + // `m:` and so never contain a `-` left of the loss. + let (core, loss) = match rest.split_once('-') { + Some((core, expr)) => { + let loss = NeutralLoss::from_expression(expr)?.ok_or_else(|| { + IonParsingError::UnsupportedNeutralLoss { + loss: expr.to_string(), + } + })?; + (core, loss) + } + None => (rest, NeutralLoss::None), + }; + + Ok(( + core, + Suffixes { + charge, + isotope, + loss, + }, + )) +} + +/// Parse a whole annotation, minus any mass-error suffix, into its parts. +/// +/// A comma-separated list of alternatives is NOT handled here -- that is +/// ambiguity, and resolving it needs the caller's policy. Split on `,` and +/// parse each alternative. +pub(crate) fn parse_ion(value: &str) -> Result<(IonSeriesOrdinal, Suffixes), IonParsingError> { + let (core, suffixes) = peel_suffixes(value)?; + Ok((IonSeriesOrdinal::parse(core)?, suffixes)) +} diff --git a/rust/micromzpaf/src/series.rs b/rust/micromzpaf/src/series.rs new file mode 100644 index 00000000..025b1ae6 --- /dev/null +++ b/rust/micromzpaf/src/series.rs @@ -0,0 +1,422 @@ +//! What kind of ion an annotation names, and how that fits in a packed payload. +//! +//! This module owns three related details: +//! +//! 1. the mzPAF **spelling** of each kind (`Display` and `IonSeriesOrdinal::parse`, +//! which are exact inverses), +//! 2. the **payload layout** each kind packs into (`IonSeriesOrdinal::to_parts` +//! and `from_parts`, also inverses), +//! 3. the **bounds** that layout imposes, enforced by the constructors. +//! +//! [`IonAnnot`](crate::IonAnnot) owns the surrounding word -- charge, isotope, +//! loss, and their positions. + +use std::fmt::Display; + +use crate::{ + IonParsingError, + mask, +}; + +/// Width of the `kind` discriminant in the packed word. +pub(crate) const KIND_BITS: u32 = 4; +/// Width of the `payload` this module encodes into. +pub(crate) const PAYLOAD_BITS: u32 = 12; +/// Width of each internal-fragment endpoint inside `payload`. Two share it. +const INTERNAL_POS_BITS: u32 = 6; +/// Width of the immonium residue index inside `payload`. +const IMMONIUM_BITS: u32 = 5; + +/// Widest residue position an internal fragment endpoint holds. +/// +/// Narrower than a backbone ordinal because two endpoints share one payload. +/// Internal fragments are bounded by peptide length, so 63 is well past tryptic. +pub const INTERNAL_POS_MAX: u8 = mask(INTERNAL_POS_BITS) as u8; + +const _: () = assert!( + 2 * INTERNAL_POS_BITS <= PAYLOAD_BITS, + "two internal-fragment endpoints must fit in one payload" +); +const _: () = assert!( + IMMONIUM_BITS <= PAYLOAD_BITS && mask(IMMONIUM_BITS) >= (b'Z' - b'A') as u32, + "the immonium field must fit the payload and hold every uppercase residue" +); + +/// `kind` discriminants. The numbering lives only here; [`IonSeriesOrdinal::to_parts`] +/// and `from_parts` are its only readers. +/// +/// Backbone ions take 1..=9 from [`Series`] itself, so a series can never drift +/// out of step with its letter. `unknown` is 0 so that the all-zero word decodes +/// to a real annotation. 13..=15 are unassigned and decode as `unknown`. +const KIND_UNKNOWN: u32 = 0; +const KIND_PRECURSOR: u32 = 10; +const KIND_INTERNAL: u32 = 11; +const KIND_IMMONIUM: u32 = 12; + +const _: () = assert!( + KIND_IMMONIUM <= mask(KIND_BITS), + "a kind discriminant overflows KIND_BITS" +); + +/// A backbone fragment ion series. +/// +/// The nine mzPAF backbone series differ only by their letter, so the +/// letter-to-discriminant pairing lives in exactly one private table and every +/// match over them is a single arm. +#[derive(Clone, Copy, Eq, PartialEq, Ord, PartialOrd, Hash, Debug)] +#[allow(non_camel_case_types)] +#[repr(u8)] +pub enum Series { + a = 1, + b, + c, + d, + v, + w, + x, + y, + z, +} + +impl Series { + /// Every series, in discriminant order. + pub const ALL: [Self; 9] = [ + Self::a, + Self::b, + Self::c, + Self::d, + Self::v, + Self::w, + Self::x, + Self::y, + Self::z, + ]; + /// The mzPAF letters, in discriminant order. + const CHARS: &'static [u8; 9] = b"abcdvwxyz"; + + /// The mzPAF letter for this series. + pub const fn as_char(self) -> char { + Self::CHARS[self as usize - 1] as char + } + + /// The series for an mzPAF letter, or `None` if it names no backbone series. + fn from_char(c: char) -> Option { + let idx = Self::CHARS.iter().position(|&b| b as char == c)?; + Some(Self::ALL[idx]) + } +} + +impl Display for Series { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.as_char()) + } +} + +/// What kind of ion an annotation names, with whatever that kind carries. +/// +/// A *view*: [`IonAnnot`](crate::IonAnnot) stores a packed word and reconstructs +/// this on demand, so building one allocates no annotation. Fields are public +/// because reading them is the point; the validating constructors +/// ([`Self::try_internal`], [`Self::try_immonium`]) are how you get one that is known to +/// fit the payload. +#[derive(Clone, Eq, PartialEq, Ord, PartialOrd, Hash, Debug, Copy)] +#[allow(non_camel_case_types)] +pub enum IonSeriesOrdinal { + /// One of the nine backbone series, at position `ordinal` in the ladder. + backbone { + series: Series, + ordinal: u8, + }, + /// An unannotated peak. `ordinal` is a uniqueness counter, not a position. + unknown { + ordinal: u8, + }, + precursor, + /// An internal fragment spanning residues `start..=end`. + internal { + start: u8, + end: u8, + }, + /// A bare immonium ion for an uppercase residue code. + immonium { + residue: char, + }, +} + +impl IonSeriesOrdinal { + /// An internal fragment spanning residues `start..=end`. + /// + /// Endpoints are bounded by [`INTERNAL_POS_MAX`], which is this module's + /// payload width and is checked here rather than by the caller. + pub fn try_internal(start: u8, end: u8) -> Result { + if start > INTERNAL_POS_MAX || end > INTERNAL_POS_MAX { + return Err(IonParsingError::OrdinalOutOfRange { + ordinal: start.max(end), + series: 'm', + }); + } + Ok(Self::internal { start, end }) + } + + /// A bare immonium ion for an uppercase residue code. + /// + /// Modified immonium (`IC[Carbamidomethyl]`) carries an arbitrary mod string + /// that no fixed-width field holds, so only a bare uppercase residue is + /// accepted. + pub fn try_immonium(residue: char) -> Result { + if !residue.is_ascii_uppercase() { + return Err(IonParsingError::UnsupportedFragmentType { + fragment_type: residue, + }); + } + Ok(Self::immonium { residue }) + } + + /// A backbone / precursor / unknown ion from its mzPAF letter. + /// + /// `p` is the only letter that takes no ordinal; every other one requires + /// one. Internal fragments and immonium ions are spelled differently and + /// have their own constructors. + pub fn from_series_char(c: char, ordinal: Option) -> Result { + if c == 'p' { + return match ordinal { + None => Ok(Self::precursor), + Some(ordinal) => Err(IonParsingError::UnexpectedOrdinal { series: c, ordinal }), + }; + } + let ordinal = ordinal.ok_or(IonParsingError::MissingOrdinal { series: c })?; + if c == '?' { + return Ok(Self::unknown { ordinal }); + } + Series::from_char(c) + .map(|series| Self::backbone { series, ordinal }) + .ok_or(IonParsingError::UnsupportedFragmentType { fragment_type: c }) + } + + /// Parse the ion-kind part of an annotation: everything left of the charge, + /// isotope and loss suffixes. + /// + /// The inverse of this type's [`Display`]. + pub(crate) fn parse(core: &str) -> Result { + // Internal fragment: `m:`. Checked before the single-letter + // forms because `m` is not a backbone letter, so there is no ambiguity. + if let Some(spans) = core.strip_prefix('m') + && let Some((a, b)) = spans.split_once(':') + { + let start = a + .parse::() + .map_err(|_| IonParsingError::parse(core, "internal-fragment start is not a u8"))?; + let end = b + .parse::() + .map_err(|_| IonParsingError::parse(core, "internal-fragment end is not a u8"))?; + return Self::try_internal(start, end); + } + + // Immonium: `I`, bare only. + if let Some(res) = core.strip_prefix('I') { + let mut ch = res.chars(); + return match (ch.next(), ch.next()) { + (Some(r), None) => Self::try_immonium(r), + (None, _) => Err(IonParsingError::parse( + core, + "immonium ion names no residue", + )), + // Carries a mod string no fixed-width field can hold. + _ => Err(IonParsingError::UnsupportedModifiedImmonium { + annotation: core.to_string(), + }), + }; + } + + // Backbone / precursor / unknown: a series letter then an ordinal. + let mut chars = core.chars(); + let series = chars + .next() + .ok_or(IonParsingError::parse(core, "empty annotation"))?; + let rest = chars.as_str(); + let ordinal = + if rest.is_empty() { + None + } else { + Some(rest.parse::().map_err(|_| { + IonParsingError::parse(core, "ordinal is not a number in 0..=255") + })?) + }; + Self::from_series_char(series, ordinal) + } + + /// Split into the `kind` discriminant and its `payload`, the two fields + /// [`IonAnnot`](crate::IonAnnot) packs. + /// + /// Split into the fields used by the packed word. + pub(crate) const fn to_parts(self) -> (u32, u32) { + match self { + Self::backbone { series, ordinal } => (series as u32, ordinal as u32), + Self::unknown { ordinal } => (KIND_UNKNOWN, ordinal as u32), + Self::precursor => (KIND_PRECURSOR, 0), + Self::internal { start, end } => ( + KIND_INTERNAL, + ((start as u32) & mask(INTERNAL_POS_BITS)) + | (((end as u32) & mask(INTERNAL_POS_BITS)) << INTERNAL_POS_BITS), + ), + // `saturating_sub` rather than `-`: the constructor rejects a + // non-uppercase residue, but the variant is public and this must not + // underflow for one built directly. + Self::immonium { residue } => ( + KIND_IMMONIUM, + (residue as u32).saturating_sub(b'A' as u32) & mask(IMMONIUM_BITS), + ), + } + } + + /// Inverse of [`Self::to_parts`]. Unknown discriminants decode as unknown + /// ions. + pub(crate) const fn from_parts(kind: u32, payload: u32) -> Self { + let ordinal = payload as u8; + match kind { + 1..=9 => Self::backbone { + // `kind` is in range, so this is the inverse of `series as u32`. + series: Series::ALL[kind as usize - 1], + ordinal, + }, + KIND_PRECURSOR => Self::precursor, + KIND_INTERNAL => Self::internal { + start: (payload & mask(INTERNAL_POS_BITS)) as u8, + end: ((payload >> INTERNAL_POS_BITS) & mask(INTERNAL_POS_BITS)) as u8, + }, + KIND_IMMONIUM => Self::immonium { + residue: (b'A' + (payload & mask(IMMONIUM_BITS)) as u8) as char, + }, + _ => Self::unknown { ordinal }, + } + } +} + +impl Display for IonSeriesOrdinal { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::backbone { series, ordinal } => write!(f, "{}{}", series.as_char(), ordinal), + Self::unknown { ordinal } => write!(f, "?{}", ordinal), + Self::precursor => write!(f, "p"), + Self::internal { start, end } => write!(f, "m{}:{}", start, end), + Self::immonium { residue } => write!(f, "I{}", residue), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn all_series() -> Vec { + let mut out: Vec = Series::ALL + .iter() + .enumerate() + .map(|(i, &series)| IonSeriesOrdinal::backbone { + series, + ordinal: i as u8 + 1, + }) + .collect(); + out.extend([ + IonSeriesOrdinal::precursor, + IonSeriesOrdinal::unknown { ordinal: 10 }, + IonSeriesOrdinal::internal { start: 2, end: 11 }, + IonSeriesOrdinal::internal { + start: INTERNAL_POS_MAX, + end: INTERNAL_POS_MAX, + }, + IonSeriesOrdinal::immonium { residue: 'A' }, + IonSeriesOrdinal::immonium { residue: 'Z' }, + ]); + out + } + + #[test] + fn every_variant_round_trips_through_the_packed_parts() { + for series in all_series() { + let (kind, payload) = series.to_parts(); + assert_eq!(IonSeriesOrdinal::from_parts(kind, payload), series); + assert!(kind <= mask(KIND_BITS), "{series:?} kind overflows"); + assert!( + payload <= mask(PAYLOAD_BITS), + "{series:?} payload overflows" + ); + } + + let mut kinds: Vec = Series::ALL + .iter() + .map(|&series| { + IonSeriesOrdinal::backbone { series, ordinal: 1 } + .to_parts() + .0 + }) + .collect(); + kinds.sort_unstable(); + kinds.dedup(); + assert_eq!(kinds.len(), Series::ALL.len(), "two series share a kind"); + } + + #[test] + fn every_kind_bit_pattern_decodes_without_panicking() { + for kind in 0..=mask(KIND_BITS) { + for payload in [0, 1, mask(PAYLOAD_BITS)] { + let _ = IonSeriesOrdinal::from_parts(kind, payload).to_string(); + } + } + } + + #[test] + fn every_variant_round_trips_through_its_mzpaf_spelling() { + for series in all_series() { + let text = series.to_string(); + assert_eq!( + IonSeriesOrdinal::parse(&text).expect("own spelling parses"), + series, + "{text}" + ); + } + } + + #[test] + fn series_letters_and_discriminants_agree() { + assert_eq!(Series::ALL.len(), Series::CHARS.len()); + for &series in &Series::ALL { + assert_eq!(Series::from_char(series.as_char()), Some(series)); + } + for c in ['p', '?', 'm', 'I', 'q', 'A'] { + assert_eq!(Series::from_char(c), None, "{c} is not a backbone series"); + } + } + + #[test] + fn constructors_reject_what_the_payload_cannot_hold() { + assert!(IonSeriesOrdinal::try_internal(INTERNAL_POS_MAX, INTERNAL_POS_MAX).is_ok()); + assert!(matches!( + IonSeriesOrdinal::try_internal(INTERNAL_POS_MAX + 1, 1), + Err(IonParsingError::OrdinalOutOfRange { .. }) + )); + assert!(IonSeriesOrdinal::try_immonium('A').is_ok()); + assert!(IonSeriesOrdinal::try_immonium('Z').is_ok()); + assert!(IonSeriesOrdinal::try_immonium('a').is_err()); + + // A bare `I` names no residue; a modified one carries a mod string. + assert!(IonSeriesOrdinal::parse("I").is_err()); + assert!(matches!( + IonSeriesOrdinal::parse("IC[Carbamidomethyl]"), + Err(IonParsingError::UnsupportedModifiedImmonium { .. }) + )); + // `p` is the only letter that takes no ordinal, and requires none. + assert!(matches!( + IonSeriesOrdinal::parse("p1"), + Err(IonParsingError::UnexpectedOrdinal { .. }) + )); + assert!(matches!( + IonSeriesOrdinal::parse("y"), + Err(IonParsingError::MissingOrdinal { .. }) + )); + assert!(matches!( + IonSeriesOrdinal::parse("q1"), + Err(IonParsingError::UnsupportedFragmentType { .. }) + )); + } +} diff --git a/rust/rescore_dash/src/app.rs b/rust/rescore_dash/src/app.rs index 85468857..26e3a27d 100644 --- a/rust/rescore_dash/src/app.rs +++ b/rust/rescore_dash/src/app.rs @@ -294,7 +294,7 @@ fn available() -> bool { /// would unwind past the caller's warn-only `if let Err` -- after the results /// have already been written to disk. /// -/// The terminal is restored on every returning path. [`catch_panics`] covers +/// The terminal is restored on every returning path. `catch_panics` covers /// the event loop, but only in a debug build; see its doc. pub fn run(dash: Dashboard) -> std::io::Result<()> { if !available() { diff --git a/rust/rescore_dash/src/lib.rs b/rust/rescore_dash/src/lib.rs index e90fed8c..a747f7c4 100644 --- a/rust/rescore_dash/src/lib.rs +++ b/rust/rescore_dash/src/lib.rs @@ -2,7 +2,8 @@ //! target/decoy separation, FDR curve and decoy calibration. //! //! Two steps, deliberately separate. [`Dashboard::build`] materializes -//! everything on screen from a [`RescoreView`] -- see [`precompute`] for what is +//! everything on screen from a [`RescoreView`] -- see the `precompute` module +//! for what is //! exact and what is sampled -- and [`run`] opens the TUI over the result. //! Splitting them lets the caller drop the feature matrix, gigabytes at a //! realistic library size, before the TUI blocks for as long as the user leaves diff --git a/rust/speclib_build_cli/src/config.rs b/rust/speclib_build_cli/src/config.rs index b399c92c..7d301154 100644 --- a/rust/speclib_build_cli/src/config.rs +++ b/rust/speclib_build_cli/src/config.rs @@ -91,9 +91,11 @@ impl Default for DigestionConfig { #[derive(Debug, Deserialize)] #[serde(default, deny_unknown_fields)] pub struct ModificationsConfig { - /// Fixed modifications applied to every matching residue, e.g. ["Carbamidomethyl@C"]. + /// Fixed modifications applied to every matching residue, e.g. + /// `["Carbamidomethyl@C"]`. pub fixed: Vec, - /// Variable modifications considered during peptide generation, e.g. ["Oxidation@M"]. + /// Variable modifications considered during peptide generation, e.g. + /// `["Oxidation@M"]`. pub variable: Vec, /// Maximum number of variable modifications per peptide. pub max_variable: usize, diff --git a/rust/speclib_build_cli/src/entry.rs b/rust/speclib_build_cli/src/entry.rs index 68a7668e..3d253596 100644 --- a/rust/speclib_build_cli/src/entry.rs +++ b/rust/speclib_build_cli/src/entry.rs @@ -30,7 +30,7 @@ pub struct EntryFilters { /// Strip bracket-enclosed modifications from a sequence. /// -/// "PEPTC[U:4]IDEK" → "PEPTCIDEK" +/// `"PEPTC[U:4]IDEK"` -> `"PEPTCIDEK"` pub fn strip_mods(seq: &str) -> String { let mut out = String::with_capacity(seq.len()); let mut depth = 0usize; diff --git a/rust/speclib_build_cli/src/mods.rs b/rust/speclib_build_cli/src/mods.rs index a1c7d7cb..f7edfb69 100644 --- a/rust/speclib_build_cli/src/mods.rs +++ b/rust/speclib_build_cli/src/mods.rs @@ -1,6 +1,6 @@ /// Proforma-like modification parsing and application for speclib_build. /// -/// Supports notations like "C[U:4]", "M[U:35]", "S[U:21]", "M[+15.995]". +/// Supports notations like `C[U:4]`, `M[U:35]`, `S[U:21]`, `M[+15.995]`. /// Fixed mods are inserted after every matching residue in the sequence. /// Variable mods generate all combinations up to `max_mods` sites. @@ -10,7 +10,7 @@ pub struct Modification { /// Uppercase single-letter amino acid code. pub target_residue: char, - /// Bracket notation including brackets, e.g. "[U:4]" or "[+15.995]". + /// Bracket notation including brackets, e.g. `[U:4]` or `[+15.995]`. pub notation: String, } diff --git a/rust/timsquery/src/lib.rs b/rust/timsquery/src/lib.rs index e0f2c8d9..65aa71aa 100644 --- a/rust/timsquery/src/lib.rs +++ b/rust/timsquery/src/lib.rs @@ -58,7 +58,9 @@ pub mod ion { IonAnnot, IonParsingError, IonSeriesOrdinal, - IonSeriesTerminality, + NeutralLoss, + Series, + UnknownIonCounter, }; } diff --git a/rust/timsquery/src/serde/diann_io.rs b/rust/timsquery/src/serde/diann_io.rs index 9d18d66b..02251d65 100644 --- a/rust/timsquery/src/serde/diann_io.rs +++ b/rust/timsquery/src/serde/diann_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use crate::models::OwnedSourceId; use arrow::array::{ @@ -26,7 +27,7 @@ use tracing::{ pub enum DiannReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(DiannPrecursorParsingError), Parquet(String), Arrow(String), } @@ -36,8 +37,8 @@ impl std::fmt::Display for DiannReadingError { match self { DiannReadingError::Io => write!(f, "IO error"), DiannReadingError::Csv => write!(f, "CSV parsing error"), - DiannReadingError::PrecursorParsing => { - write!(f, "DIA-NN precursor parsing error") + DiannReadingError::PrecursorParsing(err) => { + write!(f, "DIA-NN precursor parsing error: {}", err) } DiannReadingError::Parquet(msg) => write!(f, "Parquet error: {}", msg), DiannReadingError::Arrow(msg) => write!(f, "Arrow error: {}", msg), @@ -49,7 +50,11 @@ impl std::error::Error for DiannReadingError {} #[derive(Debug)] pub enum DiannPrecursorParsingError { - IonParsingError, + /// Ion parsing failed at this fragment row. + IonParsing { + row: usize, + source: IonParsingError, + }, /// A library that names its other precursors left this one blank. See /// [`Naming`]. UnnamedPrecursor, @@ -58,16 +63,32 @@ pub enum DiannPrecursorParsingError { Other, } -impl From for DiannPrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - DiannPrecursorParsingError::IonParsingError +impl std::fmt::Display for DiannPrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => { + write!(f, "fragment row {}: {}", row, source) + } + Self::UnnamedPrecursor => write!( + f, + "a library that names its other precursors left this one blank" + ), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::EmptyIonString => write!(f, "empty FragmentType"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl DiannPrecursorParsingError { + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for DiannReadingError { - fn from(_err: DiannPrecursorParsingError) -> Self { - DiannReadingError::PrecursorParsing + fn from(err: DiannPrecursorParsingError) -> Self { + DiannReadingError::PrecursorParsing(err) } } @@ -85,7 +106,7 @@ impl From for DiannReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct DiannPrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -413,7 +434,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(rows.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -433,8 +454,9 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -457,7 +479,8 @@ fn parse_precursor_group( DiannPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); @@ -743,7 +766,7 @@ fn parse_precursor_group_from_parquet( let mut fragment_mzs = Vec::with_capacity(indices.len()); buffers.fragment_labels.clear(); let mut rel_intensities = Vec::with_capacity(indices.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::default(); for (i, &idx) in indices.iter().enumerate() { let fragment_mz = columns.product_mzs[idx] as f64; @@ -765,8 +788,9 @@ fn parse_precursor_group_from_parquet( columns.fragment_loss_types[idx], i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); rel_intensities.push((ion_annot, rel_intensity)); @@ -789,7 +813,8 @@ fn parse_precursor_group_from_parquet( DiannPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(DiannPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); @@ -868,7 +893,9 @@ AAAAAAALQAK\tAAAAAAALQAK\t478.7\t2\t11.0\t0.9\tP2\t0\t300.0\ty\t3\t1\tnoloss\t1. assert!( matches!( read_targets(blank.path()), - Err(DiannReadingError::PrecursorParsing) + Err(DiannReadingError::PrecursorParsing( + DiannPrecursorParsingError::UnnamedPrecursor + )) ), "a blank name in a naming library must not load" ); @@ -966,29 +993,28 @@ AAAAAAALQAK\tAAAAAAALQAK\t478.7\t2\t11.0\t0.9\tP2\t0\t300.0\ty\t3\t1\tnoloss\t1. let mgrysgk = &elution_groups[1].0; let hgdtgrr = &elution_groups[0].0; - let mut mgrysgk_expected_labels = vec!["y6", "b6", "b3", "b5", "y4"] - .into_iter() - .map(|s| IonAnnot::try_from(s).unwrap()) - .collect::>(); - let mut actual_labels: Vec = mgrysgk - .iter_fragments() - .map(|(label, _mz)| *label) - .collect(); - mgrysgk_expected_labels.sort(); - actual_labels.sort(); - assert_eq!(actual_labels, mgrysgk_expected_labels); - - let mut hgdtgrr_expected_labels = vec!["y6", "b3", "y4", "y5"] - .into_iter() - .map(|s| IonAnnot::try_from(s).unwrap()) - .collect::>(); - let mut actual_labels: Vec = hgdtgrr - .iter_fragments() - .map(|(label, _mz)| *label) - .collect(); - hgdtgrr_expected_labels.sort(); - actual_labels.sort(); - assert_eq!(actual_labels, hgdtgrr_expected_labels); + // Labels as a set: the reader's fragment order is not part of the + // contract being tested here. Keyed on the mzPAF spelling because + // `IonAnnot` is deliberately not `Ord`. + fn label_set(eg: &Target) -> Vec { + let mut out: Vec = eg + .iter_fragments() + .map(|(label, _mz)| label.to_string()) + .collect(); + out.sort(); + out + } + fn expected_set(labels: &[&str]) -> Vec { + let mut out: Vec = labels.iter().map(|s| s.to_string()).collect(); + out.sort(); + out + } + + assert_eq!( + label_set(mgrysgk), + expected_set(&["y6", "b6", "b3", "b5", "y4"]) + ); + assert_eq!(label_set(hgdtgrr), expected_set(&["y6", "b3", "y4", "y5"])); } #[test] diff --git a/rust/timsquery/src/serde/diann_speclib_io.rs b/rust/timsquery/src/serde/diann_speclib_io.rs index 6c36efa3..4659857b 100644 --- a/rust/timsquery/src/serde/diann_speclib_io.rs +++ b/rust/timsquery/src/serde/diann_speclib_io.rs @@ -9,11 +9,11 @@ //! order, so there is no random access across sections. //! //! The reader is three layers: -//! 1. Decode ([`Cursor`], [`Fragment`], [`Peptide`]): typed, zero-copy views +//! 1. Decode (`Cursor`, `Fragment`, `Peptide`): typed, zero-copy views //! over byte ranges of the in-memory file. No domain knowledge. -//! 2. Emit ([`SpecLib`] + [`EntryIter`]): a pull parser that walks the -//! variable-length envelope and yields one [`EntryView`] per precursor. -//! 3. Map ([`map_entry`]): [`EntryView`] -> rows pushed directly into a +//! 2. Emit (`SpecLib` + `EntryIter`): a pull parser that walks the +//! variable-length envelope and yields one `EntryView` per precursor. +//! 3. Map (`map_entry`): `EntryView` -> rows pushed directly into a //! columnar `TargetColumns` (plus a parallel reference-intensity //! sidecar). This is where the ion-type table, dedup, and drop stats live. //! @@ -29,7 +29,11 @@ use super::library_file::{ TargetReadingError, TargetTable, }; -use crate::ion::IonAnnot; +use crate::ion::{ + IonAnnot, + IonSeriesOrdinal, + NeutralLoss, +}; use crate::models::{ Row, TargetCapabilities, @@ -569,9 +573,11 @@ fn bounded_capacity(count: usize, remaining_bytes: usize, min_record_bytes: usiz #[derive(Debug, Default, Clone, Copy)] pub struct SpeclibDecodeStats { /// Fragments flagged `ExcludeFromAssay` (`type & 0x80`). Kept (see - /// [`map_entry`]); counted for reporting only. + /// `map_entry`); counted for reporting only. pub exclude_flagged: usize, - /// Fragments with a neutral loss (`loss != 0`) -- `IonAnnot` cannot hold loss. + /// Fragments whose neutral-loss code `loss_from_code` does not map. Water + /// and ammonia are mapped and kept; anything else is dropped rather than + /// guessed, and counted so the cost stays visible. pub loss_dropped: usize, /// Fragments whose ion-type code or recovered series was unusable. pub unknown_ion_dropped: usize, @@ -700,6 +706,17 @@ fn append_arena(dst: &mut TargetColumns, mut src: TargetColumns Option { + match code { + 0 => Some(NeutralLoss::None), + 1 => Some(NeutralLoss::Water), + 2 => Some(NeutralLoss::Ammonia), + _ => None, + } +} + /// Map a decoded target `EntryView` directly into the columnar arena, folding /// drop counters into `stats` and appending each kept fragment's reference /// intensity to `frag_intens` IN THE SAME ORDER as the pushed fragment labels @@ -764,11 +781,18 @@ fn map_entry( if f.typ() & 0x80 != 0 { stats.exclude_flagged += 1; } - // IonAnnot cannot represent neutral loss; drop lossy fragments. - if f.loss() != 0 { - stats.loss_dropped += 1; - continue; - } + let loss = match loss_from_code(f.loss()) { + Some(loss) => loss, + None => { + warn!( + "speclib entry {:?}: unknown neutral-loss code {}; dropping fragment", + name, + f.loss() + ); + stats.loss_dropped += 1; + continue; + } + }; let type_char = match f.typ() & 0x7F { 1 => 'b', @@ -807,11 +831,27 @@ fn map_entry( continue; } - let ion = match IonAnnot::try_new(type_char, Some(series as u8), f.charge() as i8, 0) { + // The wire value is unsigned; reject values that do not fit IonAnnot. + let charge = match i8::try_from(f.charge()) { + Ok(charge) => charge, + Err(_) => { + warn!( + "speclib entry {:?}: fragment charge byte {} is not a charge; dropping", + name, + f.charge() + ); + stats.unknown_ion_dropped += 1; + continue; + } + }; + + let ion = match IonSeriesOrdinal::from_series_char(type_char, Some(series as u8)) + .and_then(|series| IonAnnot::new(series, loss, charge, 0)) + { Ok(ion) => ion, Err(e) => { warn!( - "speclib entry {:?}: failed to build IonAnnot ({:?}); dropping fragment", + "speclib entry {:?}: failed to build IonAnnot ({}); dropping fragment", name, e ); stats.unknown_ion_dropped += 1; @@ -872,7 +912,7 @@ type ParsedSpeclib = (TargetColumns, Vec, SpeclibDecodeStats, boo /// structural desync. /// /// Reads the file, parses the header, then parses+maps entries in parallel -/// ([`SpecLib::open`] + [`SpecLib::parse_parallel`]). +/// (`SpecLib::open` + `SpecLib::parse_parallel`). pub fn parse_speclib_reader(reader: R) -> Result { SpecLib::open(reader)?.parse_parallel() } @@ -1134,11 +1174,15 @@ mod tests { let file = std::fs::File::open(fixture_path()).unwrap(); let (geom, _intens, stats, _eof) = parse_speclib_reader(std::io::BufReader::new(file)).unwrap(); - // Reference parser: 8384 ExcludeFromAssay-flagged (kept, counted only), - // 152 neutral-loss dropped, no dc != 0 entries. assert_eq!(stats.exclude_flagged, 8384); - assert_eq!(stats.loss_dropped, 152); assert_eq!(stats.decoys_dropped, 0); + assert_eq!(stats.loss_dropped, 0); + let with_loss = geom + .frag_labels + .iter() + .filter(|l| l.loss() != NeutralLoss::None) + .count(); + assert_eq!(with_loss, 152, "every dropped loss is now a labelled one"); // Exclude-flagged fragments are present in the output, not dropped: // entry[0] keeps its flagged y9. assert!( @@ -1148,4 +1192,77 @@ mod tests { "flagged y9 must be kept" ); } + + #[test] + fn loss_codes_are_water_and_ammonia() { + use std::collections::BTreeMap; + let mut buf = Vec::new(); + std::fs::File::open(fixture_path()) + .unwrap() + .read_to_end(&mut buf) + .unwrap(); + let lib = SpecLib::open(std::io::Cursor::new(&buf[..])).unwrap(); + + let mut hist: BTreeMap = BTreeMap::new(); + let mut deltas: BTreeMap> = BTreeMap::new(); + + let mut c = Cursor::new(lib.data.bytes()); + c.pos = lib.entries_start; + for _ in 0..lib.n_entries { + let entry = read_entry(&mut c, lib.version, &lib.pg_ids).unwrap(); + let frags: Vec<(u8, u8, u8, u8, f64)> = entry + .peptide + .fragments() + .map(|f| { + ( + f.typ() & 0x7F, + f.index(), + f.charge(), + f.loss(), + f.mz() as f64, + ) + }) + .collect(); + for &(_, _, _, loss, _) in &frags { + *hist.entry(loss).or_default() += 1; + } + for &(t, i, c, loss, mz) in &frags { + if loss == 0 || c == 0 { + continue; + } + if let Some(&(_, _, _, _, base)) = frags + .iter() + .find(|&&(t2, i2, c2, l2, _)| (t2, i2, c2, l2) == (t, i, c, 0)) + { + deltas.entry(loss).or_default().push((base - mz) * c as f64); + } + } + } + for &code in hist.keys() { + assert!( + loss_from_code(code).is_some(), + "unmapped loss code {code} in the fixture: measure it before mapping it" + ); + } + let expected = [ + (1u8, NeutralLoss::Water, 18.0106_f64), + (2, NeutralLoss::Ammonia, 17.0265), + ]; + assert_eq!( + deltas.keys().copied().collect::>(), + expected.iter().map(|e| e.0).collect::>(), + "the fixture's set of lossy codes changed" + ); + for (code, loss, mass) in expected { + assert_eq!(loss_from_code(code), Some(loss)); + let observed = &deltas[&code]; + assert!(!observed.is_empty()); + for &d in observed { + assert!( + (d - mass).abs() < 3e-3, + "code {code} measures {d:.4}, not {mass:.4} -- it is not {loss:?}" + ); + } + } + } } diff --git a/rust/timsquery/src/serde/elution_group_inputs.rs b/rust/timsquery/src/serde/elution_group_inputs.rs index 97b6533c..53f5eee8 100644 --- a/rust/timsquery/src/serde/elution_group_inputs.rs +++ b/rust/timsquery/src/serde/elution_group_inputs.rs @@ -1,4 +1,7 @@ -use crate::ion::IonAnnot; +use crate::ion::{ + IonAnnot, + UnknownIonCounter, +}; use crate::tinyvec::{ TinyVec, tiny_vec, @@ -10,10 +13,19 @@ use crate::{ #[derive(Debug)] pub enum ElutionGroupInputError { - MismatchedFragmentLabelsLength { expected: usize, found: usize }, + MismatchedFragmentLabelsLength { + expected: usize, + found: usize, + }, AlreadyHasFragmentLabels, - IonConversionError { inner: String }, + IonConversionError { + inner: String, + }, MissingFragmentLabels, + /// Too many fragments for the unknown-ion label space. + TooManyFragmentsToLabel { + found: usize, + }, } /// User-friendly format for specifying elution groups in an input file @@ -44,6 +56,11 @@ impl ElutionGroupInput { if self.fragment_labels.is_some() { return Err(ElutionGroupInputError::AlreadyHasFragmentLabels); } + if num_fragments > u8::MAX as usize { + return Err(ElutionGroupInputError::TooManyFragmentsToLabel { + found: num_fragments, + }); + } let fragment_labels: Vec = (0..num_fragments).map(|i| i as u8).collect(); Ok(ElutionGroupInput { @@ -58,24 +75,30 @@ impl ElutionGroupInput { }) } + /// Fill in unique unknown-ion labels for an input with no fragment labels. pub fn try_fill_labels_annot( self, ) -> Result, ElutionGroupInputError> { - let tmp = self.try_fill_labels_u8()?; - let new_frags = tmp - .fragment_labels - .unwrap() - .into_iter() - .map(|lbl| IonAnnot::try_new('?', Some(lbl), 1, 1).unwrap()) - .collect(); + if self.fragment_labels.is_some() { + return Err(ElutionGroupInputError::AlreadyHasFragmentLabels); + } + let mut unknown_ions = UnknownIonCounter::default(); + let new_frags = self + .fragments + .iter() + .map(|_| unknown_ions.next_unknown(1)) + .collect::, _>>() + .map_err(|e| ElutionGroupInputError::IonConversionError { + inner: e.to_string(), + })?; Ok(ElutionGroupInput { - id: tmp.id, - mobility: tmp.mobility, - rt_seconds: tmp.rt_seconds, - precursor: tmp.precursor, - precursor_charge: tmp.precursor_charge, - precursor_isotopes: tmp.precursor_isotopes, - fragments: tmp.fragments, + id: self.id, + mobility: self.mobility, + rt_seconds: self.rt_seconds, + precursor: self.precursor, + precursor_charge: self.precursor_charge, + precursor_isotopes: self.precursor_isotopes, + fragments: self.fragments, fragment_labels: Some(new_frags), }) } diff --git a/rust/timsquery/src/serde/skyline_io.rs b/rust/timsquery/src/serde/skyline_io.rs index f91c4203..78683fc3 100644 --- a/rust/timsquery/src/serde/skyline_io.rs +++ b/rust/timsquery/src/serde/skyline_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use serde::{ Deserialize, @@ -20,7 +21,7 @@ use tracing::{ pub enum SkylineReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(SkylinePrecursorParsingError), } #[derive(Debug)] @@ -36,8 +37,8 @@ impl std::fmt::Display for SkylineReadingError { match self { SkylineReadingError::Io => write!(f, "IO error"), SkylineReadingError::Csv => write!(f, "CSV parsing error"), - SkylineReadingError::PrecursorParsing => { - write!(f, "Skyline precursor parsing error") + SkylineReadingError::PrecursorParsing(err) => { + write!(f, "Skyline precursor parsing error: {}", err) } } } @@ -47,21 +48,34 @@ impl std::error::Error for SkylineReadingError {} #[derive(Debug)] pub enum SkylinePrecursorParsingError { - IonParsingError, + /// Ion parsing failed at this fragment row. + IonParsing { + row: usize, + source: IonParsingError, + }, IonOverCapacity, Other, } -impl From for SkylinePrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - SkylinePrecursorParsingError::IonParsingError +impl std::fmt::Display for SkylinePrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => write!(f, "fragment row {}: {}", row, source), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl SkylinePrecursorParsingError { + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for SkylineReadingError { - fn from(_err: SkylinePrecursorParsingError) -> Self { - SkylineReadingError::PrecursorParsing + fn from(err: SkylinePrecursorParsingError) -> Self { + SkylineReadingError::PrecursorParsing(err) } } @@ -79,7 +93,7 @@ impl From for SkylineReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct SkylinePrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -309,7 +323,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(fragment_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(fragment_rows.len()); - let mut num_unknown_losses = 0u8; + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in fragment_rows.iter().enumerate() { let fragment_mz = row.product_mz; @@ -340,7 +354,8 @@ fn parse_precursor_group( ); SkylinePrecursorParsingError::IonOverCapacity })?; - IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)? + IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(SkylinePrecursorParsingError::ion(i))? } None => { warn!( @@ -348,8 +363,9 @@ fn parse_precursor_group( falling back to unknown ion", frag_type, i ); - num_unknown_losses = num_unknown_losses.saturating_add(1); - IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)? + unknown_ions + .next_unknown(frag_charge as i8) + .map_err(SkylinePrecursorParsingError::ion(i))? } }; diff --git a/rust/timsquery/src/serde/spectronaut_io.rs b/rust/timsquery/src/serde/spectronaut_io.rs index c850ce94..4b81e9ac 100644 --- a/rust/timsquery/src/serde/spectronaut_io.rs +++ b/rust/timsquery/src/serde/spectronaut_io.rs @@ -2,6 +2,7 @@ use crate::Target; use crate::ion::{ IonAnnot, IonParsingError, + UnknownIonCounter, }; use serde::Deserialize; use std::path::Path; @@ -16,7 +17,7 @@ use tracing::{ pub enum SpectronautReadingError { Io, Csv, - PrecursorParsing, + PrecursorParsing(SpectronautPrecursorParsingError), } /// Error type for library format detection (sniffing) @@ -35,8 +36,8 @@ impl std::fmt::Display for SpectronautReadingError { match self { SpectronautReadingError::Io => write!(f, "IO error"), SpectronautReadingError::Csv => write!(f, "CSV parsing error"), - SpectronautReadingError::PrecursorParsing => { - write!(f, "Spectronaut precursor parsing error") + SpectronautReadingError::PrecursorParsing(err) => { + write!(f, "Spectronaut precursor parsing error: {}", err) } } } @@ -46,22 +47,36 @@ impl std::error::Error for SpectronautReadingError {} #[derive(Debug)] pub enum SpectronautPrecursorParsingError { - IonParsingError, + /// Ion parsing failed at this fragment row. + IonParsing { + row: usize, + source: IonParsingError, + }, IonOverCapacity, EmptyIonString, Other, } -impl From for SpectronautPrecursorParsingError { - fn from(err: IonParsingError) -> Self { - error!("Ion parsing error: {:?}", err); - SpectronautPrecursorParsingError::IonParsingError +impl std::fmt::Display for SpectronautPrecursorParsingError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::IonParsing { row, source } => write!(f, "fragment row {}: {}", row, source), + Self::IonOverCapacity => write!(f, "fragment ordinal or charge out of range"), + Self::EmptyIonString => write!(f, "empty FragmentType"), + Self::Other => write!(f, "malformed precursor group"), + } + } +} + +impl SpectronautPrecursorParsingError { + fn ion(row: usize) -> impl Fn(IonParsingError) -> Self { + move |source| Self::IonParsing { row, source } } } impl From for SpectronautReadingError { - fn from(_err: SpectronautPrecursorParsingError) -> Self { - SpectronautReadingError::PrecursorParsing + fn from(err: SpectronautPrecursorParsingError) -> Self { + SpectronautReadingError::PrecursorParsing(err) } } @@ -79,7 +94,7 @@ impl From for SpectronautReadingError { } } -#[derive(Debug, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Clone, PartialEq)] pub struct SpectronautPrecursorExtras { pub modified_peptide: String, pub stripped_peptide: String, @@ -283,7 +298,7 @@ fn parse_precursor_group( let mut fragment_mzs = Vec::with_capacity(included_rows.len()); buffers.fragment_labels.clear(); let mut relative_intensities = Vec::with_capacity(included_rows.len()); - let mut num_unknown_losses = 0; + let mut unknown_ions = UnknownIonCounter::default(); for (i, row) in included_rows.iter().enumerate() { let fragment_mz = row.fragment_mz; @@ -303,8 +318,9 @@ fn parse_precursor_group( row.fragment_loss_type, i ); - num_unknown_losses += 1; - let ion_annot = IonAnnot::try_new('?', Some(num_unknown_losses), frag_charge as i8, 0)?; + let ion_annot = unknown_ions + .next_unknown(frag_charge as i8) + .map_err(SpectronautPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); relative_intensities.push((ion_annot, rel_intensity)); @@ -327,7 +343,8 @@ fn parse_precursor_group( SpectronautPrecursorParsingError::IonOverCapacity })?; - let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0)?; + let ion_annot = IonAnnot::try_new(frag_char, Some(frag_num), frag_charge as i8, 0) + .map_err(SpectronautPrecursorParsingError::ion(i))?; buffers.fragment_labels.push(ion_annot); fragment_mzs.push(fragment_mz); diff --git a/rust/timsquery/src/traits/fragment_label.rs b/rust/timsquery/src/traits/fragment_label.rs index e415a698..335e08be 100644 --- a/rust/timsquery/src/traits/fragment_label.rs +++ b/rust/timsquery/src/traits/fragment_label.rs @@ -1,5 +1,6 @@ use crate::IonAnnot; use crate::traits::KeyLike; +use micromzpaf::Series; use std::sync::Arc; /// Capability to apply a decoy m/z shift to a fragment label's m/z. @@ -20,11 +21,12 @@ pub trait DecoyShift { /// `FragmentLabel`. NOTE: decoy m/z shifting is gated on `DecoyShift`, not this /// trait -- a decoy variant computes its shift through `DecoyShift::decoy_shift_mz`, /// which every label implements. -// TODO(fragment-features): add series() accessor once a standalone series -// enum exists (today series and ordinal are fused in `IonSeriesOrdinal`), and -// wire this bound into the score that consumes it. +// TODO(fragment-features): wire this bound into the score that consumes it. pub trait FragmentLabel: KeyLike + DecoyShift { fn try_get_ordinal(&self) -> Option; + /// The backbone series, or `None` for ions that sit on no ladder + /// (precursor, unknown, internal, immonium). + fn try_get_series(&self) -> Option; fn get_charge(&self) -> i8; } @@ -45,6 +47,10 @@ impl FragmentLabel for IonAnnot { IonAnnot::try_get_ordinal(self) } + fn try_get_series(&self) -> Option { + IonAnnot::try_get_series(self) + } + fn get_charge(&self) -> i8 { IonAnnot::get_charge(self) } diff --git a/rust/timsquery/src/utils/sorting.rs b/rust/timsquery/src/utils/sorting.rs index 4bac6dea..d9ca8e04 100644 --- a/rust/timsquery/src/utils/sorting.rs +++ b/rust/timsquery/src/utils/sorting.rs @@ -1,48 +1,3 @@ -/// Macro that sorts an arbitrary number of vecs by a the values -/// first one. -/// -/// NOTE: This macro creates a new ordered vec for each one. -/// In theory its possible to have this happen in-place (see commit history) -/// but it seems ineficient for the most part when I benchmarked it. -/// -/// TODO: Make a variant that sort in parallel. -/// TODO: Add a parameter to specify ascending or descending. -/// -/// # Example -/// ``` -/// use timsquery::sort_vecs_by_first; -/// -/// let va = vec![9, 8, 7]; -/// let vb = vec![1, 2, 3]; -/// let vc = vec!['a', 'b', 'c']; -/// let out = sort_vecs_by_first!(&va, &vb, &vc); -/// -/// assert_eq!(out.0, vec![7, 8, 9]); -/// assert_eq!(out.1, vec![3, 2, 1]); -/// assert_eq!(out.2, vec!['c', 'b', 'a']); -/// ``` -#[macro_export] -macro_rules! sort_vecs_by_first { - ($first:expr $(,$rest:expr)*) => {{ - let first_vec = $first; - let len = first_vec.len(); - - // Create and sort indices - let mut indices: Vec<_> = (0..len).collect(); - indices.sort_unstable_by_key(|&i| &first_vec[i]); - - // Reorder first vector - let sorted_first: Vec<_> = indices.iter().map(|&i| first_vec[i]).collect(); - - // Reorder all other vectors - (sorted_first, $( { - let other_vec = $rest; - assert_eq!(other_vec.len(), len, "All vectors must have the same length"); - indices.iter().map(|&i| other_vec[i]).collect::>() - }, )*) - }}; -} - /// Returns the top n elements of a slice. /// /// The indices of the elements are also returned. @@ -94,45 +49,6 @@ pub fn top_n(slice: &[T], n: usize) -> (Vec, Vec mod tests { use super::*; - #[test] - fn test_sort_two_vecs() { - let v1 = vec![3, 1, 4, 1, 5]; - let v2 = vec!['a', 'b', 'c', 'd', 'e']; - - let (sorted_v1, sorted_v2) = sort_vecs_by_first!(v1, v2); - - assert_eq!(sorted_v1, vec![1, 1, 3, 4, 5]); - assert_eq!(sorted_v2, vec!['b', 'd', 'a', 'c', 'e']); - } - - #[test] - fn test_sort_three_vecs() { - let v1 = vec![3, 1, 4]; - let v2 = vec!['x', 'y', 'z']; - let v3 = vec![true, false, true]; - - let (sorted_v1, sorted_v2, sorted_v3) = sort_vecs_by_first!(v1, v2, v3); - - assert_eq!(sorted_v1, vec![1, 3, 4]); - assert_eq!(sorted_v2, vec!['y', 'x', 'z']); - assert_eq!(sorted_v3, vec![false, true, true]); - } - - #[test] - fn test_sort_four_vecs() { - let v1 = vec![3, 1, 4]; - let v2 = vec!['x', 'y', 'z']; - let v3 = vec![true, false, true]; - let v4 = vec![1.0, 2.0, 3.0]; - - let (sorted_v1, sorted_v2, sorted_v3, sorted_v4) = sort_vecs_by_first!(v1, v2, v3, v4); - - assert_eq!(sorted_v1, vec![1, 3, 4]); - assert_eq!(sorted_v2, vec!['y', 'x', 'z']); - assert_eq!(sorted_v3, vec![false, true, true]); - assert_eq!(sorted_v4, vec![2.0, 1.0, 3.0]); - } - #[test] fn test_top_n() { let v = vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; diff --git a/rust/timsquery_viewer/src/calibration.rs b/rust/timsquery_viewer/src/calibration.rs index 392e9c2e..729500c6 100644 --- a/rust/timsquery_viewer/src/calibration.rs +++ b/rust/timsquery_viewer/src/calibration.rs @@ -653,7 +653,7 @@ impl ViewerCalibrationState { /// /// `indexed_data` and `elution_groups` are needed to enable the Start /// button (we need both loaded). `tolerance` is written when the user - /// clicks [Apply]. + /// clicks `Apply`. pub fn render_panel( &mut self, ui: &mut egui::Ui, diff --git a/rust/timsseek/src/data_sources/speclib.rs b/rust/timsseek/src/data_sources/speclib.rs index 1905ea8d..45909066 100644 --- a/rust/timsseek/src/data_sources/speclib.rs +++ b/rust/timsseek/src/data_sources/speclib.rs @@ -176,7 +176,7 @@ fn strip_mods(s: &str) -> String { out } -/// Summary of a [`finalize_reference_library`] call, for load-time logging. +/// Summary of a `finalize_reference_library` call, for load-time logging. #[derive(Debug, Clone, Copy)] pub struct LoadReport { /// Physical stored rows (pre decoy expansion), i.e. `TargetColumns::n_rows`. diff --git a/rust/timsseek/src/lib.rs b/rust/timsseek/src/lib.rs index 2ea374f6..88908d39 100644 --- a/rust/timsseek/src/lib.rs +++ b/rust/timsseek/src/lib.rs @@ -33,6 +33,5 @@ pub use scoring::{ pub use timsquery::ion::{ IonAnnot, IonParsingError, - IonSeriesTerminality, }; pub use traits::ScorerQueriable; diff --git a/rust/timsseek/src/ml/cv.rs b/rust/timsseek/src/ml/cv.rs index 8ea27621..480d6ef4 100644 --- a/rust/timsseek/src/ml/cv.rs +++ b/rust/timsseek/src/ml/cv.rs @@ -354,24 +354,15 @@ impl FoldDataset for StreamingDataset<'_, T> { /// A model [`CrossValidatedScorer`] can cross-fit: fit on a row-index slice, /// score another, and report per-column importance. /// -/// `fit` receives BOTH the training rows and an early-stopping (`val`) slice. -/// A model without early stopping is expected to ignore `val`; only the -/// scorer's guarantee matters, namely that neither slice is ever scored by the -/// model fitted from it. +/// `fit` receives training rows and an optional early-stopping slice. Neither +/// slice is scored by the model fitted from it. pub(crate) trait FoldModel: Sized { type Config; type Error; /// Fit one fold's model. /// - /// `fold` is the identity of the fold being fitted, and it is a PARAMETER - /// rather than config state because the caller is the only thing that knows - /// it: a config is shared across every fold, and deriving the fold from - /// `train[0]` breaks on an empty slice and is meaningless for a partition - /// that trains on several folds at once (`qvalues::crossfit_lda` trains on - /// all folds but one). Models with a stochastic initialization mix it into - /// their seed, so that folds differ from each other while a rerun of the - /// same fold does not. Deterministic models ignore it, the same way a model - /// without early stopping ignores `val`. + /// `fold` identifies the fold being fitted. Models may use it to derive a + /// fold-specific random seed. fn fit( cfg: &Self::Config, data: &D, @@ -418,10 +409,6 @@ pub(crate) trait FoldModel: Sized { /// Dataset over an already-materialized row-major matrix. `get_values` is a /// `copy_from_slice` out of the existing slab. -/// -/// `pub(crate)`, matching [`PrecomputedFeatures`]: the only constructor is -/// `pub(crate)` and there is no other public inherent method, so `pub` made the -/// name visible outside the crate while leaving it unconstructable there. pub(crate) struct RowMajorDataset { features: PrecomputedFeatures, names: Vec>, @@ -544,12 +531,8 @@ impl DataBuffer { /// Gather the FEATURES of `rows` (in the given order) out of `data` into /// feature-major layout: `fold_buffer[feature_idx * nrows + sample_idx]`. /// - /// Deliberately does NOT read [`FoldDataset::is_decoy`]. The scoring path has - /// no use for the labels, and this module is where leak-freedom is enforced, - /// so "the labels are not even read here" is worth being structurally true - /// rather than true by inspection of a discarded binding. `response_buffer` - /// is cleared, so [`Self::as_matrix`] cannot hand out a stale label slice - /// from a previous gather; use [`Self::features_as_matrix`] after this. + /// Does not read labels. Clears the response buffer so a features-only gather + /// cannot expose labels from a previous gather. fn fill_features_from(&mut self, data: &D, rows: &[usize], ncols: usize) { self.fold_buffer.clear(); self.response_buffer.clear(); @@ -583,15 +566,7 @@ impl DataBuffer { /// Features + labels. Only valid after [`Self::fill_from`]. /// - /// The assertion catches a features-only gather OF A NON-EMPTY ROW SET -- - /// [`Self::fill_features_from`] clears `response_buffer`, so the labels are - /// missing while `nrows` is not, and this panics instead of handing out an - /// empty response slice. It does NOT catch the zero-row case: with - /// `nrows == 0` both sides are `0` and the two fills are indistinguishable - /// here, so an empty gather returns an empty view either way. That is the - /// shape `crossfit` produces (it passes `val = &[]`), and it is harmless - /// because that path uses [`Self::fill_from`] anyway -- but the assertion is - /// not what makes it safe. + /// Panics unless the label buffer matches the row count. fn as_matrix(&self) -> FoldView<'_> { let mat = self.features_as_matrix(); assert_eq!(self.response_buffer.len(), self.nrows); @@ -610,10 +585,8 @@ pub(crate) fn fold_weights(data: &D, rows: &[usize]) -> Vec .collect() } -/// [`FoldModel`] adapter for `forust`'s [`GradientBooster`]. A newtype only -/// because `GradientBooster` is a foreign type; it also carries the lane width -/// so [`FoldModel::importance`] can return a full-width, lane-indexed vector -/// (forust reports only the columns it split on). +/// [`FoldModel`] adapter for `forust`'s [`GradientBooster`]. It carries the lane +/// width so importance can be returned in full feature order. pub(crate) struct GbmFoldModel { booster: GradientBooster, ncols: usize, @@ -623,10 +596,7 @@ impl FoldModel for GbmFoldModel { type Config = GBMConfig; type Error = ForustError; - /// `fold` is DELIBERATELY IGNORED: forust seeds itself from - /// `GBMConfig::seed`, so mixing the fold in here would fight the one seed - /// the booster already has. Folds still differ, because they are fitted on - /// different rows. + /// The booster owns its seed; fold identity is not mixed into it. fn fit( cfg: &GBMConfig, data: &D, @@ -667,7 +637,7 @@ impl FoldModel for GbmFoldModel { /// forust reports only the columns it actually split on. Every other lane /// is `NAN` -- "not reported" -- rather than `0.0`, because this model has no /// gain measurement for a column it never used. See the - /// [`FoldModel::importance`] contract. + /// `FoldModel::importance` contract. fn importance(&self) -> Vec { let raw = self .booster @@ -720,10 +690,8 @@ impl PrecomputedFeatures { } /// Build from an already-materialized row-major matrix - /// (`features[i*ncols + j]`) + responses, instead of walking - /// Rows MUST align with the `data` the scorer is constructed from (same - /// order). This is how the lane-matrix consumer trains GBM on a prebuilt - /// lane feature set. + /// (`features[i*ncols + j]`) and responses. Rows must align with the + /// scorer's data order. pub(crate) fn from_row_major(features: Vec, ncols: usize, responses: Vec) -> Self { assert_eq!( features.len(), @@ -752,18 +720,8 @@ pub struct FeatureStat { pub nan_ratio: f32, } -/// Per-fold feature statistics. -/// -/// `feature_stats` is in the dataset's own column order -/// ([`FoldDataset::column_names`], i.e. the matrix's column order), one entry per -/// column. `feature_importance` is sorted by importance DESCENDING (top features -/// first), and carries only the columns the model reported a finite value for -- -/// see [`FoldModel::importance`] -- so it is generally shorter than -/// `feature_stats` and in a different order. -/// -/// Descending because the two TSV sidecars `timsseek_cli` writes -/// (`results.feature_stats.tsv` / `results.feature_importance.tsv`) emit these -/// vectors in order, one row per entry, and the importance one is read top-down. +/// Per-fold feature statistics in dataset column order. Feature importance is +/// sorted descending and includes only finite reported values. #[derive(Debug, Serialize)] pub struct FoldStats { pub fold: u8, @@ -773,20 +731,9 @@ pub struct FoldStats { pub type RescoreFeatureStats = Vec; -/// Per-fold feature means/NaN ratios + model importance, for ANY fold -/// partition. -/// -/// Partition-agnostic on purpose: `fold_rows[f]` is simply "the rows summarized -/// under fold `f`" and `models[f]` is "the model whose importance is reported -/// there". The two rescoring partitions in this crate disagree on both -- the -/// [`CrossValidatedScorer`] fits on fold `f` and scores the others, while -/// `qvalues::crossfit_lda` fits on everything BUT fold `f` and scores only fold -/// `f` -- and both are leak-free. Keeping this function ignorant of which one it -/// is handed is what lets the sidecar have one implementation without the two -/// partitions being forced to converge. -/// -/// Column names (and therefore the row width) come from `data`, so the stats -/// align with the matrix the model saw by construction. +/// Per-fold feature means, NaN ratios and model importance. `fold_rows[f]` +/// contains the rows summarized for fold `f`, and `models[f]` supplies that +/// fold's importance values. pub(crate) fn fold_feature_stats( data: &D, fold_rows: &[Vec], @@ -796,16 +743,7 @@ pub(crate) fn fold_feature_stats( let mut out: RescoreFeatureStats = Vec::with_capacity(fold_rows.len()); let mut row_buf = vec![0.0f64; names.len()]; for (fold, rows) in fold_rows.iter().enumerate() { - // --- Importance, back to the (name, gain) sidecar shape --- - // UNREPORTED (`NAN`) columns are dropped; every REPORTED column is - // emitted, `0.0` included. See the [`FoldModel::importance`] contract - // for why those are different things. The drop matters because the - // sidecar and the dashboard's fold-averaged gain both treat a feature - // as "reported by this fold" simply by being present, so a model's - // unmeasured columns would otherwise pad the averaging divisor with - // values it never produced. - // - // Retain finite zeroes: they are reported measurements. + // Keep finite importances, including zero; NaN means unreported. let importance: Vec<(Arc, f32)> = match models.get(fold).copied().flatten() { Some(model) => { let raw_imp = model.importance(); @@ -1472,17 +1410,7 @@ mod test { } } - /// THE [`FoldModel::importance`] sentinel contract, at the sidecar boundary: - /// `NAN` means "this model reports nothing for this column" and is dropped; - /// every FINITE value reaches the sidecar, `0.0` included. - /// - /// Regression guard. The boundary used to filter `!= 0.0`, which was right - /// for a tree model (where 0.0 did mean "never split on") and wrong for - /// every other model: an LDA's `|coef|` of exactly 0.0 is a measurement of a - /// dead or constant column, and silently deleting those rows removes - /// exactly what an operator reads the sidecar to find. Both halves are - /// asserted -- dropping the NAN alone would also pass if 0.0 were dropped - /// too, so the surviving zero is the load-bearing assertion. + /// NaN importances are omitted; finite values, including zero, are retained. #[test] fn importance_nan_is_unreported_but_zero_is_a_value() { let dataset = RowMajorDataset::new( @@ -1514,9 +1442,7 @@ mod test { assert_eq!(stats[0].feature_stats.len(), 3); } - /// The GBM side of the same contract: forust reports only the columns it - /// split on, and the ones it did not must come back as `NAN` (absent from - /// the sidecar), NOT as a 0.0 gain it never measured. + /// Forust reports only columns it split on; other columns are NaN. /// /// `feature_4` is constant, so no tree can split on it. The other four are /// the usual separable draws, so their presence keeps the absence diff --git a/rust/timsseek/src/ml/qvalues.rs b/rust/timsseek/src/ml/qvalues.rs index 0fefd4a9..01d81c54 100644 --- a/rust/timsseek/src/ml/qvalues.rs +++ b/rust/timsseek/src/ml/qvalues.rs @@ -230,7 +230,7 @@ fn finalize( /// A cross-fit (leak-free) model over the canonical rescore fold partition. /// -/// See [`crossfit`] for the partition contract. +/// See `crossfit` for the partition contract. struct CrossFit { /// Held-out score per row, row-aligned with the input matrix. scores: Vec, @@ -259,48 +259,13 @@ impl CrossFit { /// Cross-fit a [`FoldModel`] over the canonical rescore fold partition and /// return each row's HELD-OUT score. /// -/// The shared statement of leak-freedom for [`crossfit_lda`], `rescore_lda`, -/// and `rescore_hybrid`. -/// /// Generic over the model because the partition is independent of the fitted /// model. `what` names the model in failure logs. /// -/// # Why a row may never be scored by a model that saw it -/// These models are label-aware: they fit on the target/decoy labels. One -/// in-sample fit over all rows, scoring those same rows, lets every row's -/// discriminant peek at its own label; the target/decoy separation then looks -/// better than it is, and `assign_qval` derives the q-values from exactly that -/// separation. The result is an FDR that is wrong in the flattering direction -/// with nothing downstream to catch it. -/// -/// The hybrid needs this even more sharply: there the held-out score is not the -/// final score but one `lda_score` column fed to a -/// cross-validated GBM. An in-sample column would smuggle a row's own label -/// into a feature the GBM reads while that row is held out of its GBM fold -- so -/// the GBM's own CV cannot notice, and the leak surfaces only as an optimistic -/// FDR. -/// /// # Partition -/// Fold membership comes from [`FoldDataset::get_fold`] -- for the -/// production datasets, that is `i % N_RESCORE_FOLDS`. For -/// fold `f` the model is fitted on every row with `get_fold(i) != f` and then -/// scores only the rows with `get_fold(i) == f`, so no row contributes to the -/// model scoring it. -/// -/// This is NOT `CrossValidatedScorer`'s partition and must not be unified with -/// it: that one fits on fold `f` alone, early-stops on `f + 1`, and scores the -/// remaining `n_folds - 2` folds, so a row is scored by several models that -/// each saw a fraction of the data. Here a row is scored exactly once, by a -/// model that saw everything else. Both satisfy leak-freedom, which is the only -/// property either has to satisfy; how many rows a model sees and how many -/// models score a row are free to differ, and do. -/// -/// What the two DO have to agree on is the fold ASSIGNMENT `get_fold`, or a -/// hybrid row's cross-fit column can come from a model trained on rows the GBM -/// is holding out -- leak restored, silently. Both sides read that assignment -/// from [`FoldDataset::get_fold`], so it has one definition; the -/// `crossfit_holds_out_exactly_the_rows_the_gbm_scorer_trains_on` test pins the -/// two partitions against each other rather than against a re-typed modulo. +/// For fold `f`, the model is fitted on rows where `get_fold(i) != f` and +/// scores only rows where `get_fold(i) == f`. Fold assignment is shared with +/// the GBM scorer when this output is used as a hybrid feature. /// /// Returns an error if any fold cannot fit or score all of its held-out rows. /// Callers therefore cannot accidentally consume a partially filled score @@ -324,18 +289,10 @@ where ); for f in 0..n_folds { - // TRAIN rows = every row NOT in fold f. HELD = exactly fold f. - // Ascending in both cases, so the row order the fit reduces over is the - // dataset's own order. + // Train on every row outside fold f and score fold f. let train: Vec = (0..nrows).filter(|&i| data.get_fold(i) != f).collect(); let held: Vec = (0..nrows).filter(|&i| data.get_fold(i) == f).collect(); - // `val` is EMPTY and the partition is why: this walk trains on - // all-but-fold-`f` and scores fold `f`, so every row is already spoken - // for and there is no third slice to hand over. The LDA is closed-form - // and ignores it. `MlpFoldModel` DOES early-stop, and handles the empty - // slice by carving a deterministic inner validation set out of `train` - // -- see its `fit` for the rule; nothing about it reaches fold `f`. let model = M::fit(cfg, data, f, &train, &[]).map_err(|e| RescoreError::CrossFit { model: what, fold: f, @@ -343,7 +300,6 @@ where stage: "fit", reason: e.to_string(), })?; - // Score ONLY the held-out fold, with a model that never saw it. let preds = model .predict(data, &held) .map_err(|e| RescoreError::CrossFit { @@ -380,7 +336,7 @@ where }) } -/// Cross-fit an [`LdaModel`] on [`LdaConfig::default`] -- see [`crossfit`] for +/// Cross-fit an [`LdaModel`] on [`LdaConfig::default`] -- see `crossfit` for /// the partition and the failure policy. /// /// `N_RESCORE_FOLDS` fits instead of 1 is affordable precisely because the LDA @@ -447,14 +403,9 @@ pub fn rescore(mut data: Vec) -> RescoreResult { /// is generally much cheaper and less sensitive than GBM, with neither gap /// constant across candidate counts. /// -/// The FDR machinery (`assign_qval`, target-decoy competition) is untouched -- -/// only the discriminant score source changes. -/// -/// CROSS-FIT, not a single in-sample fit: every row's score comes from an LDA -/// fitted without that row, via [`crossfit_lda`] -- see [`crossfit`] for the -/// partition and why it is mandatory. +/// Scores each row with an LDA model fitted without that row. /// -/// Returns PER-FOLD [`FoldStats`] (one per fold, like the GBM path): feature +/// Returns per-fold `FoldStats`: feature /// means/NaN ratios over each fold's held-out rows, `|coef|` importance from /// that fold's model. /// @@ -499,10 +450,8 @@ pub fn rescore_lda(mut data: Vec) -> RescoreResult { /// /// `data` must already be through [`canonicalize_and_shuffle`]. /// -/// The fold count is `N_RESCORE_FOLDS`, i.e. the SAME -/// `get_fold` the GBM's [`CrossValidatedScorer`] derives its partition from -/// below. That shared definition is what makes "the same fold assignment on both -/// sides" structural rather than a comment; see [`crossfit`]. +/// The fold count and assignment match the GBM's +/// [`CrossValidatedScorer`] partition. fn hybrid_linear_dataset(data: &[CompetedCandidate]) -> StreamingDataset<'_, CompetedCandidate> { StreamingDataset::new( data, @@ -555,16 +504,8 @@ fn hybrid_frame( ) } -/// Hybrid rescorer: cross-fit an LDA on the LINEAR lane, push its (leak-free) -/// `lda_score` as one extra column into the NONLINEAR lane, then train the GBM -/// CV on `nonlinear + lda_score` instead of the full feature frame. -/// -/// LEAK-FREEDOM: `lda_score` is cross-fit via [`crossfit`] -- see there for the -/// partition, why a label-aware feature fed to a CV'd GBM in particular must be -/// leak-free, and why the fold ASSIGNMENT has to match the one -/// `CrossValidatedScorer` derives its own partition from. ASSIGNMENT, not -/// partition: the two train/score splits differ deliberately and both are -/// leak-free, so "unifying" them is the refactor to not make. +/// Hybrid rescorer: cross-fit an LDA on the LINEAR lane, append its held-out +/// score to the NONLINEAR lane, then train the GBM on both lanes. /// /// Selected via the `rescore_model` config field / `--rescore-model` CLI flag /// ([`crate::ml::RescoreModel::Hybrid`]). @@ -645,7 +586,7 @@ fn rescore_mlp_with(mut data: Vec, config: MlpConfig) -> Resc /// [`rescore`] gives the GBM, so the two are directly comparable. The one /// [`crate::ml::RescoreModel::Mlp`] selects. /// -/// See [`rescore_mlp_with`] for the cross-fit and determinism contracts. +/// See `rescore_mlp_with` for the cross-fit and determinism contracts. /// /// Runtime and sensitivity comparisons are not constant across candidate counts. pub fn rescore_mlp(data: Vec) -> RescoreResult { @@ -774,8 +715,7 @@ fn write_competed_linear_row(candidate: &CompetedCandidate, out: &mut [f64]) { sink.finish(); } -/// The LINEAR-lane matrix for `data` in its CURRENT order (call AFTER any -/// shuffle, so row `i` aligns with `data[i]`). `LINEAR_NCOLS` wide. +/// 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); @@ -799,12 +739,8 @@ fn build_nonlinear_matrix(data: &[CompetedCandidate]) -> Vec { /// The ALL-lane matrix (linear then nonlinear, per row) -- the GBM feature set, /// `ALL_NCOLS` wide, matching [`all_feature_name_set`]'s order. /// -/// ONE pass over `rows` with ONE `Derived::compute` per row: the two lanes are -/// adjacent within a row, so there is nothing to gain from walking twice. -/// -/// Takes `(scoring, meta)` pairs rather than a row type, because both sides of -/// rescoring feed it and they agree on nothing else. See [`competed_rows`] for -/// the pre-rescore side and [`feature_frame`] for the post-rescore one. +/// Each row is projected once, with its linear features followed by its +/// nonlinear features. fn build_all_matrix<'a>( rows: impl ExactSizeIterator, ) -> Vec { @@ -817,7 +753,7 @@ fn build_all_matrix<'a>( out } -/// Competed candidates in the shape [`build_all_matrix`] consumes. +/// Competed candidates in the shape `build_all_matrix` consumes. fn competed_rows( data: &[CompetedCandidate], ) -> impl ExactSizeIterator { @@ -833,7 +769,7 @@ pub fn feature_frame(data: &[FinalResult]) -> (Vec>, Vec) { (all_feature_name_set(), build_all_matrix(rows)) } -/// LINEAR-lane feature names (LDA), in [`project_linear_row`]'s order. +/// LINEAR-lane feature names (LDA), in `project_linear_row`'s order. pub fn linear_feature_name_set() -> Vec> { let mut n = NameSink::new(); ::linear_feature_names(&mut n); @@ -842,7 +778,7 @@ pub fn linear_feature_name_set() -> Vec> { n.into_names() } -/// NONLINEAR-lane feature names, in [`project_nonlinear_row`]'s order. The +/// NONLINEAR-lane feature names, in `project_nonlinear_row`'s order. The /// `sequence_counts` names are unconditional -- a peptide with no parsed /// sequence contributes NaN values under them, not a shorter row. pub fn nonlinear_feature_name_set() -> Vec> { @@ -855,7 +791,7 @@ pub fn nonlinear_feature_name_set() -> Vec> { } /// The ALL-lane feature names (GBM) = linear ++ nonlinear, matching -/// [`build_all_matrix`]'s column order. +/// `build_all_matrix`'s column order. pub fn all_feature_name_set() -> Vec> { let mut v = linear_feature_name_set(); v.extend(nonlinear_feature_name_set()); diff --git a/rust/timsseek/src/utils/elution_group_ops.rs b/rust/timsseek/src/utils/elution_group_ops.rs index dbd8f3e2..44fea8e7 100644 --- a/rust/timsseek/src/utils/elution_group_ops.rs +++ b/rust/timsseek/src/utils/elution_group_ops.rs @@ -16,9 +16,12 @@ pub fn apply_isotope_offset_fragments_into( ) { dst.reset_from(src); for (k, v) in dst.iter_fragments_refs_mut() { - let new_ions = k.try_with_offset_neutrons(offset).expect( - "Isotope offset overflow - this should never happen with realistic isotope offsets", - ); + // The error names the representable range, so a library carrying an + // isotope offset too close to the ceiling says what the limit is rather + // than asserting the situation is impossible. + let new_ions = k.try_with_offset_neutrons(offset).unwrap_or_else(|e| { + panic!("fragment {k} cannot take a {offset:+} isotope offset: {e}") + }); let mz_offset = (C13_C12_MASS_DIFF / k.get_charge() as f64) * offset as f64; *v += mz_offset; *k = new_ions;