diff --git a/Cargo.lock b/Cargo.lock index 06310e6..adc01ad 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1516,6 +1516,7 @@ dependencies = [ "log", "mzdata", "osprey-core", + "quick-xml", "rusqlite", "serde", "tempfile", diff --git a/Cargo.toml b/Cargo.toml index fd2f4de..d1bfd05 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -38,6 +38,12 @@ parquet = "58" # File I/O mzdata = "0.63" +# Used by osprey-io::mzml to pre-parse isolation-window cvParams as +# f64, bypassing mzdata 0.63's `param.to_f32()` quantization which +# produces ~3e-5 m/z drift on isolation bounds that land between two +# f32-representable values. mzdata transitively pins quick-xml 0.30, +# so this version matches and does not bloat the dep graph. +quick-xml = "0.30" rusqlite = { version = "0.39", features = ["bundled"] } csv = "1.3" # Use flate2's `zlib-default` feature (vendored stock zlib via libz-sys) diff --git a/crates/osprey-chromatography/src/calibration/rt.rs b/crates/osprey-chromatography/src/calibration/rt.rs index e77b544..5a05fae 100644 --- a/crates/osprey-chromatography/src/calibration/rt.rs +++ b/crates/osprey-chromatography/src/calibration/rt.rs @@ -114,13 +114,18 @@ impl RTCalibrator { ))); } - // Sort by library RT + // Sort by library RT, breaking ties on measured RT so the order is + // deterministic for duplicate-x inputs. Without the secondary key, + // pairs with equal library_rt land in input-order on Rust (stable + // sort_by) but may land in a different order cross-impl, producing + // cosmetic swaps in abs_residuals[i] downstream that the cross-impl + // calibration.json diff catches. let mut pairs: Vec<(f64, f64)> = library_rts .iter() .zip(measured_rts.iter()) .map(|(&x, &y)| (x, y)) .collect(); - pairs.sort_by(|a, b| a.0.total_cmp(&b.0)); + pairs.sort_by(|a, b| a.0.total_cmp(&b.0).then_with(|| a.1.total_cmp(&b.1))); let x: Vec = pairs.iter().map(|(x, _)| *x).collect(); let y: Vec = pairs.iter().map(|(_, y)| *y).collect(); diff --git a/crates/osprey-core/src/config.rs b/crates/osprey-core/src/config.rs index 503c42a..52c1ba1 100644 --- a/crates/osprey-core/src/config.rs +++ b/crates/osprey-core/src/config.rs @@ -580,8 +580,35 @@ n_threads: 0 # 0 = auto-detect } /// Compute SHA-256 hash of parameters that affect reconciliation. - /// Includes the search hash (if search changed, reconciliation is also invalid). + /// Includes the search hash (if search changed, reconciliation is also + /// invalid). The file set affects the consensus RTs that drive + /// reconciliation actions, so it is part of the hash. Uses + /// `self.input_files` stems by default; the per-file rescore worker + /// (`--join-at-pass=1 --no-join`) has only one input parquet in + /// `input_files` but participated in a multi-file join, so it must + /// instead compute the hash with the full file_stems list from its + /// reconciliation.json via [`Self::reconciliation_parameter_hash_for_stems`]. pub fn reconciliation_parameter_hash(&self) -> String { + let stems: Vec = self + .input_files + .iter() + .filter_map(|p| { + p.file_stem() + .and_then(|s| s.to_str()) + .map(|s| s.to_string()) + }) + .collect(); + self.reconciliation_parameter_hash_for_stems(&stems) + } + + /// Variant of [`Self::reconciliation_parameter_hash`] that accepts the + /// file_stems list explicitly. The per-file rescore worker uses this + /// with the file_stems carried in `reconciliation.json` so the hash + /// it writes into the reconciled parquet matches the hash the + /// downstream `--join-at-pass=2` merge node will compute (which is + /// based on the full set of `--input-scores` it receives, not the + /// single parquet the worker rescored). + pub fn reconciliation_parameter_hash_for_stems(&self, file_stems: &[String]) -> String { let mut hasher = Sha256::new(); hasher.update(self.search_parameter_hash().as_bytes()); hasher @@ -594,17 +621,11 @@ n_threads: 0 # 0 = auto-detect .as_bytes(), ); hasher.update(format!("run_fdr:{}\n", self.run_fdr).as_bytes()); - // File set affects consensus RTs - let mut stems: Vec = self - .input_files - .iter() - .filter_map(|p| { - p.file_stem() - .and_then(|s| s.to_str()) - .map(|s| s.to_string()) - }) - .collect(); + // File set affects consensus RTs. Sort + dedup so the same logical + // set produces the same hash regardless of caller iteration order. + let mut stems = file_stems.to_vec(); stems.sort(); + stems.dedup(); hasher.update(format!("file_stems:{:?}\n", stems).as_bytes()); format!("{:x}", hasher.finalize()) } diff --git a/crates/osprey-fdr/src/percolator.rs b/crates/osprey-fdr/src/percolator.rs index ef76050..dac427e 100644 --- a/crates/osprey-fdr/src/percolator.rs +++ b/crates/osprey-fdr/src/percolator.rs @@ -199,22 +199,83 @@ pub fn run_percolator( // two sides, cascading through every downstream computation. dump_stage5_standardizer(&standardizer, config.feature_names.as_deref()); - // 3. Subsample by peptide groups if needed (before fold splitting, per PMC5059416) - // Keeps target-decoy pairs and charge states together. - // The subsampled set is used for fold assignment + SVM training. - // ALL entries are scored with the trained model. - let train_subset: Option> = if config.max_train_size > 0 && n > config.max_train_size - { - Some(subsample_by_peptide_group( - &labels, - &entry_ids, - &peptides, - config.max_train_size, - config.seed, - )) - } else { - None - }; + // One-shot diagnostic for 2nd-pass divergence localization. Dumps + // per-entry raw feature vectors so the cross-impl compare can see + // exactly which rows differ. Sorted by (entry_id, native_position) + // to align with the C# dump. + dump_stage5_perc_input(entries, config.feature_names.as_deref()); + + // 3a. Best-per-precursor: pick the single best-scoring observation per + // (base_id, is_decoy) tuple across all files. With N files per peptide, + // this avoids the SVM seeing the same precursor's target/decoy pair + // N times, which would inflate apparent target/decoy separation and + // cause the SVM to learn file-specific noise rather than peptide + // discriminating features. The streaming path applies this dedup + // inline (pipeline.rs::run_percolator_fdr); applying it here makes + // the direct path statistically consistent. Without this, on multi- + // file inputs below the streaming threshold (Stellar 3-file at 393K + // entries) the SVM trained on N-times-redundant precursor pairs and + // produced inflated target/decoy separation diverging from the + // deduped C# port. + // + // Dedup key is features[0] (fragment_coelution_sum, the first PIN + // feature) which matches the streaming path's `coelution_sum` + // dedup field exactly — same per-entry value. + let mut best_target: HashMap = HashMap::new(); + let mut best_decoy: HashMap = HashMap::new(); + for i in 0..n { + let base_id = entries[i].entry_id & 0x7FFF_FFFF; + let score = entries[i].features[0]; + let map = if labels[i] { + &mut best_decoy + } else { + &mut best_target + }; + match map.get(&base_id) { + Some(&existing) if entries[existing].features[0] >= score => {} + _ => { + map.insert(base_id, i); + } + } + } + let mut dedup_indices: Vec = best_target + .values() + .chain(best_decoy.values()) + .copied() + .collect(); + dedup_indices.sort(); + log::debug!( + " Best-per-precursor: {} entries ({} targets, {} decoys) from {} total", + dedup_indices.len(), + best_target.len(), + best_decoy.len(), + n + ); + + // 3b. Subsample by peptide groups if dedup'd count still exceeds the + // training cap (before fold splitting, per PMC5059416). Keeps target- + // decoy pairs and charge states together. The subsampled set is used + // for fold assignment + SVM training; ALL entries are scored with + // the trained model regardless. + let train_subset: Option> = + if config.max_train_size > 0 && dedup_indices.len() > config.max_train_size { + // Build dedup-local arrays for the subsample call. + let dedup_labels: Vec = dedup_indices.iter().map(|&i| labels[i]).collect(); + let dedup_entry_ids: Vec = dedup_indices.iter().map(|&i| entry_ids[i]).collect(); + let dedup_peptides: Vec = + dedup_indices.iter().map(|&i| peptides[i].clone()).collect(); + let local = subsample_by_peptide_group( + &dedup_labels, + &dedup_entry_ids, + &dedup_peptides, + config.max_train_size, + config.seed, + ); + // Remap local indices back into the original entry index space. + Some(local.into_iter().map(|li| dedup_indices[li]).collect()) + } else { + Some(dedup_indices) + }; let sub_n = train_subset.as_ref().map_or(n, |s| s.len()); @@ -1305,7 +1366,21 @@ fn compute_per_run_peptide_qvalues( qvalues } -/// Compute experiment-level precursor q-values (across all files) +/// Compute experiment-level precursor q-values (across all files). +/// +/// Propagates the winner's q-value to every observation sharing the +/// same `base_id` (target and decoy). The streaming path +/// (`run_percolator_fdr_streaming` in `crates/osprey/src/pipeline.rs`) +/// already does this via its `base_id_exp_prec_q` map; the direct +/// path used to assign the q-value only to the single winner, +/// leaving every non-winning per-file observation at q=1.0. That +/// undercount silently broke downstream stages that gate on +/// `experiment_precursor_qvalue` (Stage 6 consensus selection / +/// calibration refit, Stage 7 protein FDR) on multi-file inputs +/// using the direct path (Stellar 3-file, ~393K entries below the +/// 600K streaming threshold). The OspreySharp port matched the +/// streaming path's "propagate to all base_id observations" +/// semantics, so the direct path was the divergent side. fn compute_experiment_precursor_qvalues( scores: &[f64], labels: &[bool], @@ -1320,8 +1395,20 @@ fn compute_experiment_precursor_qvalues( let ws: Vec = winner_indices.iter().map(|&i| scores[i]).collect(); compute_conservative_qvalues(&ws, &winner_is_decoy, &mut q); + // Build base_id -> winner q-value map (one entry per base_id, since + // compete_all picks a single winner per base_id). + let mut base_id_q: HashMap = HashMap::with_capacity(winner_indices.len()); for (rank, &idx) in winner_indices.iter().enumerate() { - qvalues[idx] = q[rank]; + let base_id = entry_ids[idx] & 0x7FFF_FFFF; + base_id_q.insert(base_id, q[rank]); + } + + // Propagate winner's q-value to every observation of that base_id. + for i in 0..n { + let base_id = entry_ids[i] & 0x7FFF_FFFF; + if let Some(&qv) = base_id_q.get(&base_id) { + qvalues[i] = qv; + } } qvalues @@ -1672,6 +1759,68 @@ fn dump_stage5_standardizer(standardizer: &FeatureStandardizer, feature_names: O exit_if_only("OSPREY_STANDARDIZER_ONLY", "Stage 5 standardizer dump"); } +/// One-shot diagnostic dump of the per-entry raw feature vectors that +/// feed `FeatureStandardizer::fit_transform`. Gated by +/// `OSPREY_DUMP_PERC_INPUT=1`. Writes `rust_stage5_perc_input.tsv` with +/// columns `native_position, entry_id, is_decoy, ..` +/// sorted by (entry_id, native_position) for stable cross-impl compare. +/// +/// This is a localizer for cross-impl standardizer divergence: when +/// `rust_stage5_standardizer.tsv` differs but `2nd-pass entries[]` +/// positions match, the divergence is in feature values themselves. +fn dump_stage5_perc_input(entries: &[PercolatorEntry], feature_names: Option<&[String]>) { + if !is_dump_enabled("OSPREY_DUMP_PERC_INPUT") { + return; + } + let path = "rust_stage5_perc_input.tsv"; + let Ok(mut f) = + std::fs::File::create(path).map(|file| std::io::BufWriter::with_capacity(8 << 20, file)) + else { + log::warn!("Could not create {}", path); + return; + }; + write!(f, "native_position\tentry_id\tis_decoy").ok(); + let n_features = if entries.is_empty() { + 0 + } else { + entries[0].features.len() + }; + for i in 0..n_features { + let name = feature_names + .and_then(|n| n.get(i)) + .map(|s| s.as_str()) + .unwrap_or("unknown"); + write!(f, "\t{}", name).ok(); + } + writeln!(f).ok(); + + let mut order: Vec = (0..entries.len()).collect(); + order.sort_by_key(|&i| (entries[i].entry_id, i)); + for i in order { + let e = &entries[i]; + write!( + f, + "{}\t{}\t{}", + i, + e.entry_id, + if e.is_decoy { "true" } else { "false" } + ) + .ok(); + for v in &e.features { + write!(f, "\t{}", format_f64_roundtrip(*v)).ok(); + } + writeln!(f).ok(); + } + let _ = f.flush(); + drop(f); + log::info!( + "Wrote Stage 5 Percolator input dump: {} ({} rows)", + path, + entries.len() + ); + exit_if_only("OSPREY_PERC_INPUT_ONLY", "Stage 5 Percolator input dump"); +} + /// Subsample entries by peptide group, keeping target-decoy pairs and charge states together. /// /// Groups entries by target peptide (via base_id), then randomly selects groups until diff --git a/crates/osprey-io/Cargo.toml b/crates/osprey-io/Cargo.toml index dfd2a89..e1f001b 100644 --- a/crates/osprey-io/Cargo.toml +++ b/crates/osprey-io/Cargo.toml @@ -11,6 +11,7 @@ repository.workspace = true [dependencies] osprey-core = { workspace = true } mzdata = { workspace = true } +quick-xml = { workspace = true } rusqlite = { workspace = true } csv = { workspace = true } flate2 = { workspace = true } diff --git a/crates/osprey-io/src/mzml/parser.rs b/crates/osprey-io/src/mzml/parser.rs index 0768a79..3daba95 100644 --- a/crates/osprey-io/src/mzml/parser.rs +++ b/crates/osprey-io/src/mzml/parser.rs @@ -14,10 +14,168 @@ use mzdata::io::mzml::MzMLReader; use mzdata::prelude::*; use mzdata::spectrum::RefPeakDataLevel; use osprey_core::{IsolationWindow, MS1Spectrum, OspreyError, Result, Spectrum, SpectrumSource}; +use quick_xml::events::Event; +use quick_xml::reader::Reader; +use std::collections::HashMap; use std::fs::File; use std::io::BufReader; use std::path::{Path, PathBuf}; +/// Raw f64 isolation window cvParam values lifted directly from the +/// mzML XML, indexed by spectrum index (matching mzdata's +/// `MultiLayerSpectrum.description().index`). Used to override +/// mzdata 0.63's f32-quantized `IsolationWindow.lower_bound` / +/// `upper_bound` values; see [`read_isolation_cvparams_f64`] for the +/// rationale. +#[derive(Default, Clone, Debug)] +struct IsolationCvParams { + /// MS:1000827 isolation window target m/z (f64 from XML text). + target_mz: Option, + /// MS:1000828 isolation window lower offset (positive magnitude). + lower_offset: Option, + /// MS:1000829 isolation window upper offset (positive magnitude). + upper_offset: Option, +} + +/// One-pass streaming scan of an mzML file to extract every +/// `` block's cvParam values at f64 precision, +/// directly from the XML text. mzdata 0.63's reader pipes these +/// values through `param.to_f32()` (see mzdata-0.63 +/// `src/io/mzml/reader.rs:281`), which quantizes through f32 and +/// produces ~3e-5 m/z drift on isolation-window edges that land +/// between two f32-representable values. The C# port +/// (OspreySharp) parses the same cvParams at f64, so the f32 +/// quantization on the Rust side is the sole source of a 1,732-row +/// cross-impl `iso_upper` mismatch on Stellar 3-file. This +/// function reads the same XML text mzdata sees but preserves full +/// f64 precision. +/// +/// The override is wired in via [`make_isolation_window`] at the two +/// MS2 parsing sites in this module. When upstream mzdata moves to +/// f64 storage for `IsolationWindow` bounds, this pre-pass + override +/// can be deleted in one commit. See workspace `Cargo.toml` +/// `quick-xml = "0.30"` line for the lifetime expectation. +fn read_isolation_cvparams_f64(path: &Path) -> Result> { + let mut reader = Reader::from_file(path).map_err(|e| { + OspreyError::MzmlParseError(format!("quick-xml open '{}': {}", path.display(), e)) + })?; + let mut buf = Vec::new(); + let mut result: HashMap = HashMap::new(); + let mut current_index: Option = None; + let mut depth_in_isolation_window: i32 = 0; + let mut current_params = IsolationCvParams::default(); + + loop { + let event = reader.read_event_into(&mut buf).map_err(|e| { + OspreyError::MzmlParseError(format!( + "quick-xml read at pos {}: {}", + reader.buffer_position(), + e + )) + })?; + match event { + Event::Start(ref e) if e.name().as_ref() == b"spectrum" => { + current_index = None; + for attr in e.attributes().flatten() { + if attr.key.as_ref() == b"index" { + if let Ok(s) = std::str::from_utf8(attr.value.as_ref()) { + current_index = s.parse().ok(); + } + } + } + current_params = IsolationCvParams::default(); + } + Event::Start(ref e) if e.name().as_ref() == b"isolationWindow" => { + depth_in_isolation_window += 1; + } + Event::End(ref e) if e.name().as_ref() == b"isolationWindow" => { + depth_in_isolation_window -= 1; + if depth_in_isolation_window == 0 { + if let Some(idx) = current_index { + result.insert(idx, current_params.clone()); + } + } + } + // cvParam elements inside isolationWindow are typically + // self-closing (`Event::Empty`); handle the start variant + // too for defensive parsing. + Event::Empty(ref e) | Event::Start(ref e) + if depth_in_isolation_window > 0 && e.name().as_ref() == b"cvParam" => + { + let mut accession: Option<&[u8]> = None; + let mut value: Option<&[u8]> = None; + let attrs: Vec<_> = e.attributes().flatten().collect(); + for attr in &attrs { + match attr.key.as_ref() { + b"accession" => accession = Some(attr.value.as_ref()), + b"value" => value = Some(attr.value.as_ref()), + _ => {} + } + } + if let (Some(acc), Some(v)) = (accession, value) { + if let (Ok(acc_str), Ok(v_str)) = + (std::str::from_utf8(acc), std::str::from_utf8(v)) + { + if let Ok(f) = v_str.parse::() { + match acc_str { + "MS:1000827" => current_params.target_mz = Some(f), + "MS:1000828" => current_params.lower_offset = Some(f), + "MS:1000829" => current_params.upper_offset = Some(f), + _ => {} + } + } + } + } + } + Event::Eof => break, + _ => {} + } + buf.clear(); + } + Ok(result) +} + +/// Build an [`IsolationWindow`] for a spectrum, preferring f64 cvParam +/// values from the pre-parsed map if available. Falls back to +/// mzdata's f32-quantized bounds when (a) no map is provided, (b) the +/// map has no entry for this scan index, or (c) the entry is missing +/// the requested cvParams (e.g., older mzML converters that omit +/// MS:1000827 / 828 / 829). +fn make_isolation_window( + precursor: &mzdata::spectrum::Precursor, + scan_number: u32, + iso_cv_params: Option<&HashMap>, +) -> Option { + let ion = precursor.ion()?; + let center_fallback = ion.mz; + + if let Some(map) = iso_cv_params { + if let Some(cv) = map.get(&scan_number) { + if cv.target_mz.is_some() || cv.lower_offset.is_some() || cv.upper_offset.is_some() { + let center = cv.target_mz.unwrap_or(center_fallback); + let lower_offset = cv.lower_offset.unwrap_or(12.5); + let upper_offset = cv.upper_offset.unwrap_or(12.5); + return Some(IsolationWindow::new(center, lower_offset, upper_offset)); + } + } + } + + // Fallback: mzdata-quantized values (f32 → f64 cast). + let isolation = &precursor.isolation_window; + let center = center_fallback; + let lower_offset = if isolation.lower_bound > 0.0 { + center - isolation.lower_bound as f64 + } else { + 12.5 + }; + let upper_offset = if isolation.upper_bound > 0.0 { + isolation.upper_bound as f64 - center + } else { + 12.5 + }; + Some(IsolationWindow::new(center, lower_offset, upper_offset)) +} + /// Some mzML producers emit peaks that are not strictly ascending in m/z /// (observed in a HeLa Astral 3 mz DIA file: ~1 row in 1.7M had a single /// inverted pair of consecutive centroids). Downstream fragment matching @@ -75,12 +233,20 @@ pub struct MzmlReader { reader: MzMLReader>, total_spectra: Option, current_index: usize, + /// Pre-parsed f64 isolation window cvParams, keyed by spectrum + /// index. Overrides mzdata 0.63's f32-quantized bounds in + /// `convert_spectrum`. See [`read_isolation_cvparams_f64`]. + iso_cv_params: HashMap, } impl MzmlReader { /// Open an mzML file for reading pub fn open>(path: P) -> Result { let path = path.as_ref().to_path_buf(); + // First pass: scan the XML for isolation-window cvParams as + // f64. Cheap streaming pass (~5% of bytes parsed) that lets us + // override mzdata 0.63's f32 quantization downstream. + let iso_cv_params = read_isolation_cvparams_f64(&path)?; let file = File::open(&path).map_err(|e| { OspreyError::MzmlParseError(format!("Failed to open file '{}': {}", path.display(), e)) })?; @@ -93,6 +259,7 @@ impl MzmlReader { reader: mzml_reader, total_spectra: None, // Will be determined on first iteration current_index: 0, + iso_cv_params, }) } @@ -116,31 +283,18 @@ impl MzmlReader { scan.start_time // mzdata returns time in minutes }); - // Get isolation window from precursor (now a Vec in mzdata 0.63) - let isolation_window = if let Some(precursor) = desc.precursor.first() { - let ion = match precursor.ion() { - Some(i) => i, - None => return Ok(None), - }; - let isolation = &precursor.isolation_window; - - let center = ion.mz; - // In mzdata 0.63, bounds are f32 directly - let lower_offset = if isolation.lower_bound > 0.0 { - center - isolation.lower_bound as f64 - } else { - 12.5 - }; - let upper_offset = if isolation.upper_bound > 0.0 { - isolation.upper_bound as f64 - center - } else { - 12.5 - }; - - IsolationWindow::new(center, lower_offset, upper_offset) - } else { - // No precursor info - skip this spectrum - return Ok(None); + // Get isolation window from precursor (now a Vec in mzdata 0.63). + // Routes through `make_isolation_window` so we use f64 cvParam + // values from the pre-parsed XML, bypassing mzdata 0.63's f32 + // quantization. See `read_isolation_cvparams_f64`. + let isolation_window = match desc.precursor.first() { + Some(precursor) => { + match make_isolation_window(precursor, scan_number, Some(&self.iso_cv_params)) { + Some(w) => w, + None => return Ok(None), // No selected ion + } + } + None => return Ok(None), // No precursor info }; // Get peaks - use the peaks() method which returns RefPeakDataLevel @@ -333,6 +487,8 @@ impl MS1Index { /// Returns: (MS2 spectra, MS1 index) pub fn load_all_spectra>(path: P) -> Result<(Vec, MS1Index)> { let path = path.as_ref(); + // First pass: f64 isolation-window cvParams from raw XML. + let iso_cv_params = read_isolation_cvparams_f64(path)?; let file = File::open(path).map_err(|e| { OspreyError::MzmlParseError(format!("Failed to open file '{}': {}", path.display(), e)) })?; @@ -400,30 +556,18 @@ pub fn load_all_spectra>(path: P) -> Result<(Vec, MS1In }); } 2 => { - // Process MS2 spectrum - // Get isolation window from precursor - let isolation_window = if let Some(precursor) = desc.precursor.first() { - let ion = match precursor.ion() { - Some(i) => i, - None => continue, - }; - let isolation = &precursor.isolation_window; - - let center = ion.mz; - let lower_offset = if isolation.lower_bound > 0.0 { - center - isolation.lower_bound as f64 - } else { - 12.5 - }; - let upper_offset = if isolation.upper_bound > 0.0 { - isolation.upper_bound as f64 - center - } else { - 12.5 - }; - - IsolationWindow::new(center, lower_offset, upper_offset) - } else { - continue; + // Process MS2 spectrum. + // Route through `make_isolation_window` for the f64 + // cvParam override path; see `convert_spectrum` and + // `read_isolation_cvparams_f64` for the rationale. + let isolation_window = match desc.precursor.first() { + Some(precursor) => { + match make_isolation_window(precursor, scan_number, Some(&iso_cv_params)) { + Some(w) => w, + None => continue, // No selected ion + } + } + None => continue, }; // Get peaks diff --git a/crates/osprey-scoring/src/diagnostics.rs b/crates/osprey-scoring/src/diagnostics.rs index 695bad9..10035e8 100644 --- a/crates/osprey-scoring/src/diagnostics.rs +++ b/crates/osprey-scoring/src/diagnostics.rs @@ -53,11 +53,16 @@ pub fn dump_cal_match(library: &[LibraryEntry], results: &[CalibrationMatch]) { for entry in entries { if let Some(m) = by_id.get(&entry.id) { - // Use :.10 everywhere so we don't hit banker's vs round- - // half-up rounding differences between Rust and C#. + // Use :.17 (17 fractional digits) so f64 values round-trip + // exactly. Earlier :.10 was chosen "to avoid banker's vs + // round-half-up rounding differences between Rust and C#", + // but evidence shows the formatters still disagree at the + // 10th decimal when f64 values land near a rounding + // boundary. :.17 sidesteps the disagreement entirely: + // the underlying f64 bits are reproduced exactly. writeln!( f, - "{}\t{}\t{}\t1\t{}\t{:.10}\t{:.10}\t{:.10}\t{}\t{:.10}\t{:.10}", + "{}\t{}\t{}\t1\t{}\t{:.17}\t{:.17}\t{:.17}\t{}\t{:.17}\t{:.17}", entry.id, if entry.is_decoy { 1 } else { 0 }, entry.charge, @@ -97,8 +102,11 @@ pub fn dump_cal_match(library: &[LibraryEntry], results: &[CalibrationMatch]) { /// Dump per-entry LDA discriminant + q-value to `rust_lda_scores.txt`, /// sorted by entry_id for a stable diff against `cs_lda_scores.txt`. /// -/// Uses `{:.10}` formatting to avoid banker's-vs-half-up text rounding -/// mismatches with C#. +/// Uses `{:.17}` formatting so the underlying f64 round-trips exactly. +/// Earlier `{:.10}` was chosen "to avoid banker's vs half-up rounding +/// mismatches with C#"; in practice Rust `{:.10}` (RHE) and .NET +/// Framework `F10` (HAFZ) still disagree by 1 in the last digit on +/// boundary-case f64s. See the cal_match dump for the same fix. /// /// Gated by `OSPREY_DUMP_LDA_SCORES=1`. When `OSPREY_LDA_SCORES_ONLY=1` /// is also set, exits the process after writing. @@ -115,7 +123,7 @@ pub fn dump_lda_scores(matches: &[CalibrationMatch]) { let m = &matches[i]; writeln!( f, - "{}\t{}\t{:.10}\t{:.10}", + "{}\t{}\t{:.17}\t{:.17}", m.entry_id, if m.is_decoy { 1 } else { 0 }, m.discriminant_score, diff --git a/crates/osprey-scoring/src/lib.rs b/crates/osprey-scoring/src/lib.rs index fc2c184..2cfeead 100644 --- a/crates/osprey-scoring/src/lib.rs +++ b/crates/osprey-scoring/src/lib.rs @@ -430,20 +430,25 @@ pub fn compute_cosine_at_scan( return 0.0; } - // L2 normalize - let lib_norm = lib_preprocessed.iter().map(|x| x * x).sum::().sqrt(); - let obs_norm = obs_preprocessed.iter().map(|x| x * x).sum::().sqrt(); + // Single-pass accumulation of squared norms + dot product, then divide once + // at the end. Matches `cosine_angle` in this crate and the C# port's + // ComputeCosineAtScan, ensuring cross-impl bit-equality for sg_weighted_cosine. + let mut lib_sq = 0.0; + let mut obs_sq = 0.0; + let mut dot = 0.0; + for (a, b) in lib_preprocessed.iter().zip(obs_preprocessed.iter()) { + lib_sq += a * a; + obs_sq += b * b; + dot += a * b; + } + let lib_norm = lib_sq.sqrt(); + let obs_norm = obs_sq.sqrt(); if lib_norm < 1e-10 || obs_norm < 1e-10 { return 0.0; } - // Cosine = dot(lib_norm, obs_norm) - lib_preprocessed - .iter() - .zip(obs_preprocessed.iter()) - .map(|(a, b)| (a / lib_norm) * (b / obs_norm)) - .sum() + dot / (lib_norm * obs_norm) } /// Compute per-fragment mass accuracy statistics at the apex scan. @@ -2081,43 +2086,14 @@ impl SpectralScorer { // First get LibCosine for additional metrics let lib_cosine_score = self.lib_cosine(observed, library); - // Bin observed spectrum using Comet BIN macro - let mut obs_binned = vec![0.0f32; self.bin_config.n_bins]; - for (&mz, &intensity) in observed.mzs.iter().zip(observed.intensities.iter()) { - if let Some(bin) = self.bin_config.mz_to_bin(mz) { - // Apply sqrt transformation to experimental spectrum - obs_binned[bin] += intensity.sqrt(); - } - } - - // Apply windowing normalization (Comet's MakeCorrData) - let windowed = self.apply_windowing_normalization(&obs_binned); - - // Apply sliding window subtraction (fast XCorr preprocessing) - let xcorr_preprocessed = self.apply_sliding_window(&windowed); - - // XCorr = sum of preprocessed experimental values at UNIQUE fragment - // bin positions. When two library fragments fall into the same bin, - // the bin's contribution must count once, not twice -- the Comet - // theoretical spectrum uses unit intensity per bin (see - // preprocess_library_for_xcorr, which sets binned[bin] = 1.0 per - // unique bin). Summing preprocessed[bin] once per fragment instead - // of once per unique bin double-counts collisions and over-scores - // dense fragment lists. - let n_bins = xcorr_preprocessed.len(); - let mut visited = vec![false; n_bins]; - let mut xcorr_raw: f64 = 0.0; - for frag in &library.fragments { - if let Some(bin) = self.bin_config.mz_to_bin(frag.mz) { - if !visited[bin] { - visited[bin] = true; - xcorr_raw += xcorr_preprocessed[bin] as f64; - } - } - } - - // Scale XCorr (pyXcorrDIA uses 0.005 for spectrum-centric) - let xcorr_scaled = xcorr_raw * 0.005; + // XCorr via the same cache-narrowing path C# OspreySharp's + // calibration uses (XcorrFromPreprocessed(float[] cache, entry)): + // f64-internal preprocessing narrows to an f32 cache, then the + // sparse sum widens each cached value back to f64 on read and + // accumulates in f64. Bit-equal cross-impl as long as the f32 + // cache is bit-equal (Option 3 guarantees this) and both impls + // promote-on-read into f64 accumulators. + let xcorr_scaled = self.xcorr_at_scan(observed, library); SpectralScore { xcorr: xcorr_scaled, @@ -2281,20 +2257,12 @@ impl SpectralScorer { matches } - /// Apply Comet-style windowing normalization - /// - /// Divides spectrum into 10 windows and normalizes each to max=50.0 - fn apply_windowing_normalization(&self, spectrum: &[f32]) -> Vec { - let mut result = vec![0.0f32; spectrum.len()]; - self.apply_windowing_normalization_into(spectrum, &mut result); - result - } - - /// In-place variant of [`apply_windowing_normalization`] that writes - /// into a caller-provided buffer. The buffer must be zero on entry - /// because below-threshold and empty-window positions retain their - /// initial value. - fn apply_windowing_normalization_into(&self, spectrum: &[f32], result: &mut [f32]) { + /// In-place Comet-style windowing normalization (f64 internal, + /// dual-purpose helper for both calibration and the HRAM per-window + /// cache path). The buffer must be zero on entry because + /// below-threshold and empty-window positions retain their initial + /// value. Matches C# `ApplyWindowingNormalizationD` bit-for-bit. + fn apply_windowing_normalization_into(&self, spectrum: &[f64], result: &mut [f64]) { debug_assert_eq!( spectrum.len(), result.len(), @@ -2303,23 +2271,20 @@ impl SpectralScorer { let num_windows = 10; let window_size = (spectrum.len() / num_windows) + 1; - // Find global max for threshold - let global_max = spectrum.iter().cloned().fold(0.0f32, f32::max); + let global_max = spectrum.iter().cloned().fold(0.0f64, f64::max); let threshold = global_max * 0.05; for window_idx in 0..num_windows { let start = window_idx * window_size; let end = ((window_idx + 1) * window_size).min(spectrum.len()); - // Find max in this window - let mut window_max = 0.0f32; + let mut window_max = 0.0f64; for &val in &spectrum[start..end] { if val > window_max { window_max = val; } } - // Normalize this window to 50.0 if window_max > 0.0 { let norm_factor = 50.0 / window_max; for i in start..end { @@ -2331,34 +2296,24 @@ impl SpectralScorer { } } - /// Apply sliding window subtraction for fast XCorr (Comet-style) - /// - /// Uses prefix sum for O(n) performance instead of O(n × window). - /// Comet divides by (2*offset) = 150 regardless of boundary effects. - /// offset=75, matching Comet's iXcorrProcessingOffset default. - fn apply_sliding_window(&self, spectrum: &[f32]) -> Vec { - let n = spectrum.len(); - let mut prefix = vec![0.0f32; n + 1]; - let mut result = vec![0.0f32; n]; - self.apply_sliding_window_into(spectrum, &mut prefix, &mut result); - result - } - - /// In-place variant of [`apply_sliding_window`] that writes into - /// caller-provided prefix and result buffers. Every position of - /// `result` is overwritten; `prefix[0]` is assumed to be 0 on entry - /// and positions 1..=n are overwritten. - fn apply_sliding_window_into(&self, spectrum: &[f32], prefix: &mut [f32], result: &mut [f32]) { + /// In-place sliding-window subtraction (Comet fast XCorr) with f64 + /// internal math narrowing to an f32 result buffer at the final store. + /// `spectrum` and `prefix` are f64 scratch; `result` is the per-spectrum + /// f32 cache. Every position of `result` is overwritten; `prefix[0]` is + /// assumed to be 0 on entry and positions 1..=n are overwritten. + /// Matches C# `ApplySlidingWindowD` for the f64 cascade, with a single + /// final cast to f32 on store. Used by the HRAM per-window cache path. + fn apply_sliding_window_into(&self, spectrum: &[f64], prefix: &mut [f64], result: &mut [f32]) { let n = spectrum.len(); let offset: usize = 75; - // Comet uses (window_size - 1) = 2*offset = 150 as divisor - let norm_factor = 1.0f32 / (2 * offset) as f32; + // Comet uses (window_size - 1) = 2*offset = 150 as divisor. + // f64 literal so the final subtraction stays in f64. + let norm_factor = 1.0f64 / (2 * offset) as f64; debug_assert_eq!(prefix.len(), n + 1, "prefix buffer must be n+1"); debug_assert_eq!(result.len(), n, "result buffer must be n"); debug_assert_eq!(prefix[0], 0.0, "prefix[0] must be 0 on entry"); - // Build prefix sum for O(n) window sums for i in 0..n { prefix[i + 1] = prefix[i] + spectrum[i]; } @@ -2366,12 +2321,12 @@ impl SpectralScorer { for i in 0..n { let left = i.saturating_sub(offset); let right = if i + offset < n { i + offset + 1 } else { n }; - // Window sum including center let window_sum = prefix[right] - prefix[left]; - // Subtract center to get sum excluding center let sum_excluding_center = window_sum - spectrum[i]; - // Subtract local average from center value - result[i] = spectrum[i] - sum_excluding_center * norm_factor; + let centered = spectrum[i] - sum_excluding_center * norm_factor; + // f64 -> f32 narrowing at the final store: single deterministic + // rounding, identical bits cross-impl when the f64 inputs agree. + result[i] = centered as f32; } } @@ -2393,11 +2348,17 @@ impl SpectralScorer { } /// In-place XCorr preprocess using a pool-rented scratch and a - /// caller-provided output buffer. The output buffer is fully - /// overwritten. The scratch's accumulator fields (`binned`, - /// `windowed`) are zeroed at the start of the call, so callers can - /// reuse a single rented scratch across many preprocess calls - /// without interleaving recycle/rent. + /// caller-provided f32 output buffer (the per-spectrum HRAM cache). + /// The scratch fields are f64 so the windowing / sliding-window + /// cascade runs in f64 precision; only the final cache write narrows + /// to f32. Matches C# OspreySharp's PreprocessSpectrumForXcorrF32 with + /// f64-internal math: same f32 cache size, but values carry single- + /// cast precision rather than f32-cascade noise. + /// + /// The scratch's accumulator fields (`binned`, `windowed`) are zeroed + /// at the start of the call, so callers can reuse a single rented + /// scratch across many preprocess calls without interleaving + /// recycle/rent. pub fn preprocess_spectrum_for_xcorr_into( &self, spectrum: &Spectrum, @@ -2417,10 +2378,13 @@ impl SpectralScorer { scratch.binned.fill(0.0); scratch.windowed.fill(0.0); - // Bin observed spectrum with sqrt transformation using Comet BIN macro + // Bin observed spectrum with sqrt transformation using Comet BIN + // macro. Widen f32 intensity to f64 BEFORE sqrt to match C#'s + // `Math.Sqrt((double)float)` bit-for-bit; f32::sqrt and + // f64::sqrt(f32 as f64) can differ by 1 ULP at f32 magnitude. for (&mz, &intensity) in spectrum.mzs.iter().zip(spectrum.intensities.iter()) { if let Some(bin) = self.bin_config.mz_to_bin(mz) { - scratch.binned[bin] += intensity.sqrt(); + scratch.binned[bin] += (intensity as f64).sqrt(); } } @@ -2467,11 +2431,15 @@ impl SpectralScorer { /// zero*something. Iterating fragments directly is O(n_frags) and /// allocates only a small stack of bin indices for deduplication. /// - /// Bit-parity with the dense path: the final scale and cast follow the - /// exact pattern `(sum_f32 * 0.005_f32) as f64` so both return the same - /// numerical value; the only difference is the summation ORDER (linear - /// left-to-right vs SIMD reduction), which typically differs by less - /// than 1e-7 on xcorr values around 1.0. + /// Cross-impl alignment: each f32 cache value is widened to f64 on read + /// and accumulated in f64; final scale uses an f64 literal. Matches + /// C# `SpectralScorer.XcorrFromPreprocessed(float[], LibraryEntry)` + /// which sums `xcorrRaw += preprocessed[bin]` into a `double` accumulator + /// (implicit float-to-double promotion on read). The iteration order + /// differs (Rust sorts bins ascending; C# iterates fragments in input + /// order with a visited[bin] dedup) but at f64 precision the resulting + /// sum-order difference is bounded by ~f64_eps × intermediate magnitude + /// (~1e-13 for typical xcorr) — well below f64 epsilon noise. #[inline] pub fn xcorr_sparse(&self, spectrum_preprocessed: &[f32], entry: &LibraryEntry) -> f64 { // Small stack-friendly accumulator for unique fragment bins. At HRAM @@ -2488,17 +2456,17 @@ impl SpectralScorer { bins.dedup(); let n_bins = spectrum_preprocessed.len(); - let sum: f32 = bins + let sum: f64 = bins .iter() .filter_map(|&b| { if b < n_bins { - Some(spectrum_preprocessed[b]) + Some(spectrum_preprocessed[b] as f64) } else { None } }) .sum(); - (sum * 0.005f32) as f64 + sum * 0.005_f64 } /// Preprocess a library entry for XCorr (Comet-style) diff --git a/crates/osprey-scoring/src/xcorr_pool.rs b/crates/osprey-scoring/src/xcorr_pool.rs index 0e62949..e0f61c3 100644 --- a/crates/osprey-scoring/src/xcorr_pool.rs +++ b/crates/osprey-scoring/src/xcorr_pool.rs @@ -1,16 +1,22 @@ //! Buffer pool for XCorr preprocessing scratch space. //! -//! XCorr preprocessing allocates three `NBins`-sized `f32` buffers per call -//! (binned accumulator, windowed result, sliding-window prefix sum). On -//! HRAM binning (NBins ~= 100K) each buffer is ~400 KB, so per-window -//! scoring across ~1000 spectra churns through ~1.2 GB of transient -//! `Vec` allocations plus the ~400 MB per-window preprocessed cache -//! itself. All of this memory is dropped at the end of each window and -//! re-allocated by the next window -- pure churn. +//! XCorr preprocessing allocates three `NBins`-sized `f64` scratch buffers +//! per call (binned accumulator, windowed result, sliding-window prefix +//! sum) and writes the final preprocessed values into a per-spectrum `f32` +//! cache. The f64 scratch is dual-purpose: it keeps the windowing / +//! sliding-window cascade math in f64 precision (avoiding the ~1e-5 +//! prefix-sum error f32 would accumulate over ~100K bins), and only the +//! final write to the cache narrows to f32. The per-spectrum cache stays +//! f32 to keep HRAM memory budget at ~400 KB per spectrum. //! -//! This pool hands out reusable scratch bundles and single output buffers -//! so that after the first few windows reach the high-water mark the hot -//! path never allocates. Mirrors the C# OspreySharp `XcorrScratchPool` in +//! On HRAM binning (NBins ~= 100K) each scratch buffer is ~800 KB, so per- +//! window scoring across ~1000 spectra would otherwise churn through +//! ~2.4 GB of transient `Vec` allocations plus the ~400 MB per-window +//! preprocessed cache. This pool hands out reusable scratch bundles and +//! single output buffers so that after the first few windows reach the +//! high-water mark the hot path never allocates. +//! +//! Mirrors the C# OspreySharp `XcorrScratchPool` in //! `pwiz_tools/OspreySharp/OspreySharp.Scoring/XcorrScratchPool.cs`. //! //! Thread-safety: the pool uses `Mutex>` bags. Contention is low @@ -24,16 +30,16 @@ use std::sync::Mutex; /// Rented by value, returned by value via [`XcorrScratchPool::recycle`]. pub struct XcorrScratch { /// `NBins` accumulator for binned-and-sqrt-transformed spectrum - /// intensities. Written via `+=` so must start zeroed; zeroed on - /// recycle. - pub binned: Vec, - /// `NBins` windowing-normalization output. Filled by + /// intensities, in f64 so the windowing pipeline preserves precision. + /// Written via `+=` so must start zeroed; zeroed on recycle. + pub binned: Vec, + /// `NBins` windowing-normalization output, f64. Filled by /// `apply_windowing_normalization_into`; zeroed on recycle so /// below-threshold positions retain 0. - pub windowed: Vec, - /// `NBins + 1` sliding-window prefix sum. `prefix[0]` is always 0; + pub windowed: Vec, + /// `NBins + 1` sliding-window prefix sum, f64. `prefix[0]` is always 0; /// positions `1..=NBins` are fully overwritten each call. - pub prefix: Vec, + pub prefix: Vec, } impl XcorrScratch { diff --git a/crates/osprey/src/pipeline.rs b/crates/osprey/src/pipeline.rs index da8a1f8..f08c800 100644 --- a/crates/osprey/src/pipeline.rs +++ b/crates/osprey/src/pipeline.rs @@ -62,6 +62,22 @@ use osprey_scoring::{ DecoyGenerator, DecoyMethod, Enzyme, SpectralScorer, }; +/// Resolve the precursor tolerance to a ppm value usable for MS1 isotope +/// envelope extraction. If the configured precursor tolerance is already +/// in ppm, return it as-is. If it's a Da/Mz tolerance (typical for unit- +/// resolution data like Stellar), treating it as ppm would produce an +/// absurdly tight window (e.g. 1 Da -> 0.0005 mDa at 500 m/z), making +/// `envelope.has_m0()` fail on ~99.8% of matches. Fall back to a 10 ppm +/// default in that case, matching C# OspreySharp's PerFileScoringTask +/// ScoreCalibrationEntry behavior: +/// `unit == Ppm ? tolerance : 10.0`. +fn ms1_envelope_tolerance_ppm(precursor_tolerance: &osprey_core::FragmentToleranceConfig) -> f64 { + match precursor_tolerance.unit { + osprey_core::ToleranceUnit::Ppm => precursor_tolerance.tolerance, + osprey_core::ToleranceUnit::Mz => 10.0, + } +} + /// Wrapper to implement MS1SpectrumLookup for MS1Index struct MS1IndexWrapper<'a>(&'a MS1Index); @@ -342,7 +358,24 @@ fn build_scores_metadata(config: &OspreyConfig) -> Vec Vec { +/// Build the reconciled-parquet metadata block. +/// +/// `file_stems_override`: when present, used to compute the +/// reconciliation hash via +/// [`OspreyConfig::reconciliation_parameter_hash_for_stems`]. The +/// per-file rescore worker (`--join-at-pass=1 --no-join`) passes the +/// file_stems list from reconciliation.json here so the hash it writes +/// matches the multi-file hash the downstream `--join-at-pass=2` merge +/// node will compute. Pass `None` from the in-process pipeline path — +/// `config.input_files` already reflects the full set. +fn build_reconciled_metadata( + config: &OspreyConfig, + file_stems_override: Option<&[String]>, +) -> Vec { + let recon_hash = match file_stems_override { + Some(stems) if !stems.is_empty() => config.reconciliation_parameter_hash_for_stems(stems), + _ => config.reconciliation_parameter_hash(), + }; vec![ parquet::file::metadata::KeyValue { key: META_OSPREY_VERSION.to_string(), @@ -362,7 +395,7 @@ fn build_reconciled_metadata(config: &OspreyConfig) -> Vec= rt_stats.r_squared * 0.99 { rt_calibration = rt_cal_refined; rt_stats = rt_stats_refined; + // Update metadata to reflect the refined fit that's + // actually being used (was pass 1's count). C# parity: + // PerFileScoringTask.cs reports n_refined when pass 2 is + // accepted. + num_confident_peptides = n_refined; + // Overwrite the LOESS input dump with pass 2's points + // so the diagnostic reflects the calibration actually + // used. C# overwrites unconditionally on pass 2; Rust + // previously only wrote pass 1 (line 1000), leaving + // ~960 entries unreported here on Stellar Single. + crate::diagnostics::dump_loess_input(&refined_lib_rts, &refined_meas_rts); // Re-collect mass errors from refined matches mz_qc_data = MzQCData::new(config.fragment_tolerance.unit); @@ -1719,7 +1763,22 @@ fn write_scores_parquet_with_metadata( .map(|_| Float64Builder::with_capacity(n)) .collect(); - for entry in entries { + // Iterate in canonical sorted order (entry_id, charge, scan_number) so + // per-side parquets have identical physical row layout across Rust and C# + // impls. Order-sensitive consumers downstream (Stage 5 standardizer, + // SVM training) then see the same row sequence regardless of which side + // wrote the parquet, making per-side cross-tool stages bit-equal. + let mut sorted_indices: Vec = (0..entries.len()).collect(); + sorted_indices.sort_by(|&a, &b| { + entries[a] + .entry_id + .cmp(&entries[b].entry_id) + .then_with(|| entries[a].charge.cmp(&entries[b].charge)) + .then_with(|| entries[a].scan_number.cmp(&entries[b].scan_number)) + }); + + for &idx in &sorted_indices { + let entry = &entries[idx]; entry_id_b.append_value(entry.entry_id); decoy_b.append_value(entry.is_decoy); seq_b.append_value(&entry.sequence); @@ -2838,6 +2897,13 @@ pub(crate) fn rescore_per_file_loop( file_name_to_idx: &HashMap, config: &OspreyConfig, seq_interner: &mut HashSet>, + // Full set of file stems that participated in the first-join phase + // (the join that produced the reconciliation actions this loop applies). + // Used to compute the reconciliation hash written into the reconciled + // parquet's metadata so the downstream `--join-at-pass=2` merge node + // accepts it. When empty, falls back to `config.input_files` stems + // (the in-process pipeline path always has the full set there). + join_file_stems: &[String], ) -> Result { use crate::reconciliation::ReconcileAction; let total_reconciliation: usize = reconciliation_actions @@ -3023,14 +3089,41 @@ pub(crate) fn rescore_per_file_loop( ); } - let cal_params: Option = input_file.parent().and_then(|input_dir| { + // Hard-error on missing or unreadable calibration JSON. The + // calibration sidecar is required by Stage 6 for MS1/MS2 mass + // calibration of the rescore search; silently proceeding with + // uncalibrated MS scoring would produce wrong scores without + // surfacing the cause. In-process this should always exist + // (Stage 2 just wrote it); workers already gate on it in + // `run_rescore` so this catch is belt-and-suspenders for both + // paths. + let cal_params: CalibrationParams = { + let input_dir = input_file.parent().ok_or_else(|| { + OspreyError::config(format!( + "rescore_per_file_loop: cannot derive parent directory from input path `{}`. \ + Stage 6 needs to locate the Stage 1-4 calibration sidecar next to the mzML.", + input_file.display() + )) + })?; let cal_path = calibration_path_for_input(input_file, input_dir); - if cal_path.exists() { - load_calibration(&cal_path).ok() - } else { - None + if !cal_path.exists() { + return Err(OspreyError::config(format!( + "rescore_per_file_loop: required calibration JSON not found at `{}` (input \ + file: `{}`). Stage 6 needs the Stage 1-4 calibration sidecar to rescore.", + cal_path.display(), + input_file.display() + ))); } - }); + load_calibration(&cal_path).map_err(|e| { + OspreyError::config(format!( + "rescore_per_file_loop: failed to read calibration JSON `{}`: {}. The file \ + exists but could not be parsed -- check that it was written by a matching \ + Osprey version.", + cal_path.display(), + e + )) + })? + }; // --- Re-score existing entries (consensus + reconciliation) --- let mut overlay: HashMap = HashMap::new(); @@ -3039,7 +3132,7 @@ pub(crate) fn rescore_per_file_loop( &subset_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, config, file_name, @@ -3066,7 +3159,7 @@ pub(crate) fn rescore_per_file_loop( &gap_fill_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, &gap_config, file_name, @@ -3121,7 +3214,7 @@ pub(crate) fn rescore_per_file_loop( &forced_library, &spectra, &ms1_index, - cal_params.as_ref(), + Some(&cal_params), rt_cal, config, file_name, @@ -3219,25 +3312,86 @@ pub(crate) fn rescore_per_file_loop( } } // Append gap-fill entries (remaining overlay entries beyond parquet range). - // Update the corresponding fdr_entries stubs' parquet_index to point to - // the actual Parquet row they will occupy after write-back. Without this, - // gap-fill stubs keep parquet_index = u32::MAX and Phase 2 of the next - // Percolator run can't load their features. + // The vec_idx in the overlay matches the position in fdr_entries + // (gap-fill entries were appended to fdr_entries in order). + // Track the (vec_idx, pre-sort row) pair so we can remap each + // FdrEntry stub's parquet_index to the post-canonical-sort row. let gap_start_pq_row = full_entries.len(); let mut gap_entries: Vec<(usize, CoelutionScoredEntry)> = overlay.into_iter().collect(); gap_entries.sort_by_key(|(idx, _)| *idx); + let mut gap_vec_idx_for_pre_sort_row: Vec = + Vec::with_capacity(gap_entries.len()); for (gap_offset, (vec_idx, entry)) in gap_entries.into_iter().enumerate() { - let new_pq_row = (gap_start_pq_row + gap_offset) as u32; + let pre_sort_row = gap_start_pq_row + gap_offset; full_entries.push(entry); - // The vec_idx in the overlay matches the position in fdr_entries - // (gap-fill entries were appended to fdr_entries in order) - if vec_idx < fdr_entries.len() { - fdr_entries[vec_idx].parquet_index = new_pq_row; + // Track gap-fill mapping by pre-sort row position. + // Length pads with `usize::MAX` for the upstream rows + // (which already have a real parquet_index pointing to + // the pre-sort row). + while gap_vec_idx_for_pre_sort_row.len() + gap_start_pq_row < pre_sort_row { + gap_vec_idx_for_pre_sort_row.push(usize::MAX); + } + gap_vec_idx_for_pre_sort_row.push(vec_idx); + } + + // Canonical (entry_id, charge, scan_number) sort permutation. + // Must match write_scores_parquet_with_metadata's internal + // sort key exactly; the writer then iterates in this same + // order, producing a parquet whose physical row layout is + // identical regardless of caller order. + // + // CRITICAL: every FdrEntry stub whose parquet_index points + // into this rewritten parquet must be remapped to its post + // -sort row. Upstream rows held a parquet_index equal to + // their pre-sort row position; the sort moves them, so the + // stale pre-sort index would silently fetch a different + // entry's features in the next Percolator pass. + let mut perm: Vec = (0..full_entries.len()).collect(); + perm.sort_by(|&a, &b| { + full_entries[a] + .entry_id + .cmp(&full_entries[b].entry_id) + .then_with(|| full_entries[a].charge.cmp(&full_entries[b].charge)) + .then_with(|| { + full_entries[a] + .scan_number + .cmp(&full_entries[b].scan_number) + }) + }); + // Inverse permutation: pre_sort_row -> post_sort_row. + let mut pre_to_post: Vec = vec![0u32; full_entries.len()]; + for (post, &pre) in perm.iter().enumerate() { + pre_to_post[pre] = post as u32; + } + // Remap every FdrEntry stub for this file. Upstream stubs + // already had parquet_index pointing to their pre-sort row. + // Gap-fill stubs still have parquet_index = u32::MAX; for + // them we look up the pre-sort row via the vec_idx mapping + // built above (vec_idx -> pre_sort_row inverse). + let mut vec_idx_to_gap_pre_sort: HashMap = HashMap::new(); + for (offset, &vec_idx) in gap_vec_idx_for_pre_sort_row.iter().enumerate() { + if vec_idx != usize::MAX { + vec_idx_to_gap_pre_sort.insert(vec_idx, gap_start_pq_row + offset); + } + } + for (vec_idx, fdr_entry) in fdr_entries.iter_mut().enumerate() { + let pre_sort_row = if fdr_entry.parquet_index == u32::MAX { + // Gap-fill stub: look up pre-sort row by vec_idx. + match vec_idx_to_gap_pre_sort.get(&vec_idx) { + Some(&row) => row, + None => continue, // gap-fill stub without an overlay entry? + } + } else { + fdr_entry.parquet_index as usize + }; + if pre_sort_row < pre_to_post.len() { + fdr_entry.parquet_index = pre_to_post[pre_sort_row]; } } + // Write back to Parquet with reconciliation metadata - let recon_metadata = build_reconciled_metadata(config); + let recon_metadata = build_reconciled_metadata(config, Some(join_file_stems)); let codec = parquet_compression_codec(config.parquet_compression); if let Err(e) = write_scores_parquet_with_metadata( cache_path, @@ -3257,7 +3411,7 @@ pub(crate) fn rescore_per_file_loop( let scores_path = per_file_cache_paths.get(file_name.as_str()); if let Some(cache_path) = scores_path { if let Ok(full_entries) = load_scores_parquet(cache_path) { - let recon_metadata = build_reconciled_metadata(config); + let recon_metadata = build_reconciled_metadata(config, Some(join_file_stems)); let codec = parquet_compression_codec(config.parquet_compression); if let Err(e) = write_scores_parquet_with_metadata( cache_path, @@ -3285,6 +3439,11 @@ pub(crate) fn rescore_per_file_loop( } pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { + // Sliding timer for [STAGE-WALL] per-stage perf markers. Reset at + // each stage boundary (stage1to4 -> stage5 -> stage6 -> stage7 -> blib). + // Parsed by Measure-Pipeline.ps1 for cross-impl perf comparison. + let mut stage_marker = std::time::Instant::now(); + // Validate config: --join-only feeds per-file state from --input-scores // parquets and ignores --input-files; default + --no-join modes both // require mzML inputs. @@ -4110,6 +4269,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { percolator::compute_fdr_from_stubs(&mut per_file_entries, config.run_fdr, None); } } else { + log::info!( + "[STAGE-WALL] stage1to4: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + // Dispatch FDR control based on method log::info!(""); log::info!( @@ -4158,11 +4323,16 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // = 1.0` and silently break protein-rescue parity in the worker. } - // Stage 5 diagnostic dump. Gated by OSPREY_DUMP_PERCOLATOR=1; exits the - // process when OSPREY_PERCOLATOR_ONLY=1 is also set. Writes all 4 q-values - // plus the SVM score and PEP for every FdrEntry, before compaction drops - // any rows, so the cross-impl diff sees both targets and decoys. + // Stage 5 diagnostic dump. Gated by OSPREY_DUMP_PERCOLATOR=1. Writes + // all 4 q-values plus the SVM score and PEP for every FdrEntry, before + // compaction drops any rows, so the cross-impl diff sees both targets + // and decoys. crate::diagnostics::dump_stage5_percolator(&per_file_entries); + // OSPREY_PERCOLATOR_ONLY exits after Stage 5 work completes, + // independently of whether the dump ran. Lets us measure production + // stage5 wall without paying the dump cost. Matches the C# decoupling + // in OspreySharp/Tasks/FirstJoinTask.cs. + osprey_core::diagnostics::exit_if_only("OSPREY_PERCOLATOR_ONLY", "Stage 5"); // First-pass protein FDR (picked-protein, Savitski 2015). // @@ -4352,13 +4522,30 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // 1st-pass-derived q-values, swap entry.score over to 2nd-pass scores so // the second-pass q-value recompute below sees the right rank ordering. // No-op if no 2nd-pass sidecar exists for a given file (best-effort). + // + // When NO 2nd-pass sidecar is present for ANY file (the HPC + // distribution case — the per-file rescore worker writes reconciled + // parquets but does NOT write 2nd-pass sidecars, since 2nd-pass + // Percolator is a cross-file join step that the per-file worker + // can't perform), we run 2nd-pass Percolator here so the merge + // node's output matches the straight-through pipeline. Without this, + // --join-at-pass=2 would use stale first-pass scores for rescored + // entries (which had their scores reset to 0 during reconciliation + // and never re-scored), and the final blib would silently lose + // ~25% of the precursors a straight-through run produces. if config.expect_reconciled_input { let mut reloaded = 0usize; + let mut missing_sidecars = 0usize; for (file_name, entries) in per_file_entries.iter_mut() { - if let Some(p2) = pass2_sidecar_paths.get(file_name) { + let p2 = pass2_sidecar_paths.get(file_name); + if let Some(p2) = p2 { if load_fdr_scores_sidecar(p2, entries, 2) { reloaded += 1; + } else { + missing_sidecars += 1; } + } else { + missing_sidecars += 1; } } log::debug!( @@ -4366,6 +4553,61 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { reloaded, per_file_entries.len() ); + if missing_sidecars > 0 { + log::info!( + "--join-at-pass=2: {}/{} file(s) lack a 2nd-pass sidecar — running second-pass FDR \ + to compute scores from reconciled features (HPC distribution path).", + missing_sidecars, + per_file_entries.len() + ); + log::info!(""); + log::info!("Second-pass FDR"); + let empty_overlay: RescoreOverlay = HashMap::new(); + match config.fdr_method { + FdrMethod::Percolator => { + run_percolator_fdr( + &mut per_file_entries, + &per_file_cache_paths, + &config, + &empty_overlay, + Some(&first_pass_base_ids), + )?; + } + FdrMethod::Mokapot => { + run_mokapot_fdr( + &mut per_file_entries, + &per_file_cache_paths, + &mokapot, + &pin_files, + &mokapot_dir, + &config, + )?; + } + FdrMethod::Simple => { + for (_, entries) in per_file_entries.iter_mut() { + apply_simple_fdr(entries, config.run_fdr)?; + } + } + } + crate::trace::log_fdr_qvalues(&per_file_entries, "second-pass"); + + // Persist the 2nd-pass scores we just computed as a resume + // cache for future --join-at-pass=2 invocations against the + // same reconciled parquets. Mirrors the persist call inside + // the in-process reconciliation block below; gated by + // multi-file because reconciliation only fires there. + let has_reconciliation = config.reconciliation.enabled && per_file_entries.len() > 1; + if has_reconciliation { + log::debug!("Persisting 2nd-pass FDR scores as resume cache..."); + let _ = persist_fdr_scores( + &per_file_entries, + &config, + fdr_scores_path_pass2, + "2nd-pass", + 2, + ); + } + } } // Post-FDR re-scoring: multi-charge consensus + inter-replicate reconciliation. @@ -4457,6 +4699,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { } }); + log::info!( + "[STAGE-WALL] stage5: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + // 1. Multi-charge consensus: compute per-file rescore targets // Groups by (peptide, file). If at least one charge state passes FDR, // the best-scoring charge state defines the consensus peak; other @@ -4632,6 +4880,18 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // joined work above. let search_hash = config.search_parameter_hash(); let library_hash = config.library_identity_hash(); + // The full set of file stems that participated in this first-join + // phase. Embedded in every per-file reconciliation.json so the + // downstream per-file rescore worker can compute the same + // reconciliation hash the `--join-at-pass=2` merge node will + // validate against (which is computed over all parquets at the + // 2nd-join node, not the worker's single input). Without this + // the worker writes a single-stem hash and `--join-at-pass=2` + // rejects the reconciled parquet. + let join_file_stems: Vec = per_file_entries + .iter() + .map(|(fname, _)| fname.clone()) + .collect(); // Pre-group reconciliation actions by file name to avoid the // O(num_files * num_actions) walk that the previous // implementation performed inside `from_planner_output` (one @@ -4669,6 +4929,7 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { refined_calibrations.get(file_name), &search_hash, &library_hash, + &join_file_stems, ); if let Err(e) = crate::reconciliation_io::write_reconciliation_file(&recon_path, &recon_file) @@ -4733,6 +4994,11 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { &file_name_to_idx, &config, &mut seq_interner, + // In-process pipeline: config.input_files already has the + // full set, so the empty-slice override makes + // build_reconciled_metadata fall back to + // config.reconciliation_parameter_hash(). + &[], )?; // Cross-impl bisection seam: dump the per-precursor q-values @@ -4752,6 +5018,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { total_gap_cwt, total_gap_forced, ); + log::info!( + "[STAGE-WALL] stage6: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + log::info!(""); log::info!("Second-pass FDR"); match config.fdr_method { @@ -4785,22 +5057,27 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { } } - // Persist second-pass SVM scores to sidecar files (after reconciliation block closes) + // Persist second-pass SVM scores to sidecar files (after reconciliation block closes). + // Always written when FDR ran, regardless of reconciliation. In single-file + // mode there is no rescore step, so the persisted scores equal the 1st-pass + // scores; this is fine semantically (per_file_entries already carry the + // authoritative final-pass scores) and required for cross-impl parity with + // the C# port, whose Stage 7 (--join-at-pass=2) reads the 2nd-pass sidecar + // unconditionally. Without an unconditional write here, single-file per-side + // cross-impl runs hit a load-cached-vs-retrain asymmetry in Stage 7 protein + // FDR even when every upstream stage is bit-equal. if !can_skip_fdr { - let has_reconciliation = config.reconciliation.enabled && config.input_files.len() > 1; - if has_reconciliation { - log::debug!("Persisting 2nd-pass FDR scores..."); - // 2nd-pass is a resume-only optimization (skip-Percolator on - // reruns); the return value is logged via per-file warnings - // already, so no need to escalate here. - let _ = persist_fdr_scores( - &per_file_entries, - &config, - fdr_scores_path_pass2, - "2nd-pass", - 2, - ); - } + log::debug!("Persisting 2nd-pass FDR scores..."); + // 2nd-pass is a resume-only optimization (skip-Percolator on + // reruns); the return value is logged via per-file warnings + // already, so no need to escalate here. + let _ = persist_fdr_scores( + &per_file_entries, + &config, + fdr_scores_path_pass2, + "2nd-pass", + 2, + ); } // Protein parsimony (always runs) + optional second-pass picked-protein FDR. @@ -4838,6 +5115,21 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { .map(|e| e.modified_sequence.to_string()) .collect(); + // Cross-impl bisection: dump sorted detected_peptides for diff against + // the C# port. Gated by env var; zero overhead when unset. + if std::env::var("OSPREY_DUMP_DETECTED_PEPTIDES").as_deref() == Ok("1") { + let mut sorted: Vec = detected_peptides.iter().cloned().collect(); + sorted.sort(); + let _ = std::fs::write( + "rust_stage7_detected_peptides.txt", + sorted.join("\n") + "\n", + ); + log::info!( + "[DIAG] Wrote rust_stage7_detected_peptides.txt ({} entries)", + sorted.len() + ); + } + let parsimony = protein::build_protein_parsimony( &library, config.shared_peptides, @@ -5119,6 +5411,12 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // to the final destination. This avoids SQLite locking issues on network // filesystems. if !plan_entries.is_empty() { + log::info!( + "[STAGE-WALL] stage7: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); + stage_marker = std::time::Instant::now(); + log::info!("Writing blib to {}", config.output_blib.display()); let final_path = &config.output_blib; @@ -5131,6 +5429,11 @@ pub fn run_analysis(mut config: OspreyConfig) -> Result<()> { // Move to final destination (safe copy for network filesystems) osprey_core::copy_and_verify(&temp_path, final_path)?; + + log::info!( + "[STAGE-WALL] blib: {:.1}s", + stage_marker.elapsed().as_secs_f64() + ); } else { log::warn!("No peptides passed FDR threshold, skipping blib output"); } @@ -5208,6 +5511,38 @@ fn run_percolator_fdr( log::debug!("Running native Percolator FDR on coelution entries"); + // Sort each file's entries by the parquet canonical key (entry_id, + // charge, scan_number, parquet_index) so the SVM working-set + // selection sees the same order regardless of upstream operation + // history *and* matches the parquet's physical row layout. The + // 1st-pass input is already entry_id-sorted via deduplicate_pairs + // (pipeline.rs:6123), but the post-rescore pool that feeds 2nd-pass + // Percolator can have gap-fill entries appended after the sorted + // pre-existing rows. The parquet writer uses (entry_id, charge, + // scan_number) as its sort key, but on Stellar 3-file the post- + // reconciliation pool has 94 per-file groups of 2 entries sharing + // ALL THREE of (entry_id, charge, scan_number) (gap-fill rescore + // landing on the same scan as an original row with a *different* + // rt_deviation). With only those three keys, ties leave a stable + // sort taking input order on each side -- but C# `List.Sort` is + // unstable, so the same input produced a swapped order at every + // tied group, drifting the 2nd-pass standardizer mean by 1 ULP on + // rt_deviation and cascading through all downstream SVM weights. + // Including `parquet_index` as a final tie-break gives the same + // total order on both sides regardless of underlying sort + // stability, since the parquets are byte-identical cross-impl and + // each entry's parquet_index points to its row in that shared + // canonical layout. + for (_, entries) in per_file_entries.iter_mut() { + entries.sort_by(|a, b| { + a.entry_id + .cmp(&b.entry_id) + .then_with(|| a.charge.cmp(&b.charge)) + .then_with(|| a.scan_number.cmp(&b.scan_number)) + .then_with(|| a.parquet_index.cmp(&b.parquet_index)) + }); + } + let total_entries: usize = per_file_entries.iter().map(|(_, e)| e.len()).sum(); let n_files = per_file_entries.len(); let perc_config = percolator::PercolatorConfig { diff --git a/crates/osprey/src/reconciliation_io.rs b/crates/osprey/src/reconciliation_io.rs index fa66b9b..cb295e3 100644 --- a/crates/osprey/src/reconciliation_io.rs +++ b/crates/osprey/src/reconciliation_io.rs @@ -141,9 +141,27 @@ impl<'a> Formatter for RoundtripPrettyFormatter<'a> { /// Top-level JSON envelope. Field declaration order is alphabetical so /// `serde_json::to_writer_pretty` emits keys in alphabetical order and /// matches the C# emitter byte-for-byte. +/// +/// `file_stems` (added 2026-05-19, format_version 2) carries the full +/// set of file stems that participated in the first-join phase. The +/// per-file rescore worker (`--join-at-pass=1 --no-join`) consumes only +/// its own parquet, so its `config.input_files` reflects a single +/// stem. To produce a reconciled parquet that the downstream +/// `--join-at-pass=2` merge node can validate (its +/// `reconciliation_parameter_hash` is computed over all parquets it +/// receives), the worker uses this list to compute the same multi-stem +/// hash the original join used. v1 files (no `file_stems`) deserialize +/// with an empty list; the worker will then fall back to its own +/// `config.input_files` stems, preserving v1 behavior. #[derive(Debug, Clone, Serialize, Deserialize)] #[serde(deny_unknown_fields)] pub struct ReconciliationFile { + /// Sorted list of file stems that participated in the first-join phase. + /// `#[serde(default)]` deserializes v1 files (which omit this field) as + /// an empty list; the worker code falls back to `config.input_files` + /// stems in that case for backward compatibility. + #[serde(default)] + pub file_stems: Vec, pub forced_integration_actions: Vec, pub format_version: u32, pub gap_fill_targets: Vec, @@ -198,7 +216,16 @@ pub struct RefinedRtCalibrationJson { } /// Current format version. Bump on incompatible schema changes. -pub const RECONCILIATION_FORMAT_VERSION: u32 = 1; +/// +/// v1: initial format. +/// v2: added `file_stems` so per-file rescore workers can compute the +/// reconciliation parameter hash that the downstream +/// `--join-at-pass=2` merge node expects (which is computed over +/// all files participating in the join, not the worker's single +/// parquet). Old v1 files deserialize with empty `file_stems` via +/// `#[serde(default)]`; the worker falls back to its +/// `config.input_files` stems for those (preserving v1 behavior). +pub const RECONCILIATION_FORMAT_VERSION: u32 = 2; impl ReconciliationFile { /// Build the wire envelope for a single file. Filters @@ -211,6 +238,7 @@ impl ReconciliationFile { /// should prefer [`Self::from_planner_output_pre_grouped`] together /// with a single up-front grouping to keep the total cost /// O(num_actions) rather than O(num_files * num_actions). + #[allow(clippy::too_many_arguments)] pub fn from_planner_output( file_name: &str, file_entries: &[FdrEntry], @@ -219,6 +247,7 @@ impl ReconciliationFile { refined_calibration: Option<&RTCalibration>, search_hash: &str, library_hash: &str, + file_stems: &[String], ) -> Self { let file_actions: Vec<(usize, &ReconcileAction)> = reconciliation_actions .iter() @@ -237,6 +266,7 @@ impl ReconciliationFile { refined_calibration, search_hash, library_hash, + file_stems, ) } @@ -246,6 +276,7 @@ impl ReconciliationFile { /// included or not — they are filtered here either way). Use this /// when emitting envelopes for many files so the per-file cost stays /// O(actions_for_this_file) rather than O(total_actions). + #[allow(clippy::too_many_arguments)] pub fn from_planner_output_pre_grouped( file_entries: &[FdrEntry], file_actions: &[(usize, &ReconcileAction)], @@ -253,7 +284,14 @@ impl ReconciliationFile { refined_calibration: Option<&RTCalibration>, search_hash: &str, library_hash: &str, + file_stems: &[String], ) -> Self { + // Sort and deduplicate the file_stems for deterministic output + // — matches the sort+dedup used by reconciliation_parameter_hash + // on the read side, so the worker hash matches the merge node hash. + let mut stems: Vec = file_stems.to_vec(); + stems.sort(); + stems.dedup(); let mut use_cwt: Vec = Vec::new(); let mut forced: Vec = Vec::new(); for (vec_idx, action) in file_actions { @@ -317,6 +355,7 @@ impl ReconciliationFile { }); Self { + file_stems: stems, forced_integration_actions: forced, format_version: RECONCILIATION_FORMAT_VERSION, gap_fill_targets: gap, @@ -423,6 +462,7 @@ mod tests { fn sample_file() -> ReconciliationFile { ReconciliationFile { + file_stems: vec!["fileA".to_string(), "fileB".to_string()], forced_integration_actions: vec![ ForcedIntegrationEntry { entry_id: 200, diff --git a/crates/osprey/src/rescore.rs b/crates/osprey/src/rescore.rs index b58e67b..f3d22ae 100644 --- a/crates/osprey/src/rescore.rs +++ b/crates/osprey/src/rescore.rs @@ -81,6 +81,13 @@ pub struct RescoreInputs { /// gap-fill stubs reuse the same `Arc` identities that the /// hydrated stubs already hold. pub seq_interner: HashSet>, + /// Full set of file stems that participated in the first-join phase + /// that produced the reconciliation actions in this RescoreInputs. + /// Read from `reconciliation.json`'s `file_stems` field (v2+); + /// empty for v1 files. Used to compute the reconciliation hash + /// written into the reconciled parquet's metadata so the downstream + /// `--join-at-pass=2` merge node will accept it. + pub join_file_stems: Vec, } impl RescoreInputs { @@ -131,6 +138,12 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { let mut refined_calibrations: HashMap = HashMap::new(); let mut per_file_gap_fill: HashMap> = HashMap::new(); let mut reconciliation_actions: HashMap<(String, usize), ReconcileAction> = HashMap::new(); + // Captured from the first per-file reconciliation.json envelope; every + // file's envelope must declare the same join file_stems (all came + // from the same first-join run). Used downstream to compute the + // reconciliation hash the worker stamps into the reconciled parquet's + // metadata. + let mut join_file_stems: Vec = Vec::new(); for parquet_path in parquet_paths { // The parquet stem (with `.scores` stripped) is the canonical @@ -184,6 +197,32 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { )) })?; + // Capture the join file_stems from the first envelope. Every + // file's envelope should declare the same list (they all came + // from the same first-join phase). Validate consistency + // file-to-file so an accidental mix of boundary files from + // different join runs fails loudly here rather than producing + // a reconciled parquet whose hash doesn't match anything. + if join_file_stems.is_empty() { + join_file_stems = envelope.file_stems.clone(); + join_file_stems.sort(); + join_file_stems.dedup(); + } else { + let mut these = envelope.file_stems.clone(); + these.sort(); + these.dedup(); + if !envelope.file_stems.is_empty() && these != join_file_stems { + return Err(OspreyError::config(format!( + "hydrate_for_rescore: reconciliation.json file_stems mismatch — {} \ + declares {:?}, expected {:?} (boundary files from different first-join \ + runs cannot be mixed in one rescore worker invocation).", + recon_path.display(), + these, + join_file_stems, + ))); + } + } + // 3a. Build entry_id → vec_idx map from the loaded stubs so the // planner's entry_id-keyed actions can be rehomed onto the // in-memory shape (file_name, vec_idx) the rescore engine @@ -303,6 +342,7 @@ pub fn hydrate_for_rescore(config: &OspreyConfig) -> Result { refined_calibrations, per_file_gap_fill, seq_interner, + join_file_stems, }) } @@ -383,19 +423,23 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( }; let cal_path = calibration_path_for_input(input_file, input_dir); if !cal_path.exists() { - continue; + return Err(OspreyError::config(format!( + "run_rescore: required calibration JSON not found at {} (input file: {}). \ + Stage 6 needs the Stage 1-4 calibration sidecar to rescore; without it the \ + worker would silently produce no reconciliation re-scores. Run Stages 1-4 \ + first or fix the path.", + cal_path.display(), + input_file.display() + ))); } - let cal_params = match load_calibration(&cal_path) { - Ok(p) => p, - Err(e) => { - log::warn!( - "run_rescore: failed to read calibration JSON {}: {}", - cal_path.display(), - e - ); - continue; - } - }; + let cal_params = load_calibration(&cal_path).map_err(|e| { + OspreyError::config(format!( + "run_rescore: failed to read calibration JSON {}: {}. The file exists but \ + could not be parsed — check that it was written by a matching Osprey version.", + cal_path.display(), + e + )) + })?; if let Some(ref mp) = cal_params.rt_calibration.model_params { match RTCalibration::from_model_params(mp, cal_params.rt_calibration.residual_sd) { Ok(rt_cal) => { @@ -420,6 +464,7 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( refined_calibrations, per_file_gap_fill, mut seq_interner, + join_file_stems, } = hydrate_for_rescore(&config)?; // Cross-impl bisection seam (mirrors the dump call from @@ -463,6 +508,30 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( } } + // Union with entries that have reconciliation actions. The planner + // emits actions for entries that pass FDR via cross-file consensus + // rescue (`compute_consensus_rts` upgrades a peptide if it passes + // FDR in any file in the experiment, even if it fails locally). + // Without this union, the worker would drop those entries via + // local compaction and the planner's actions would be silently + // discarded — the per-file rescore would diverge from the + // straight-through in-process pipeline. The base_ids come from + // the pre-compaction `per_file_entries` keyed by (file_name, idx) + // in `reconciliation_actions_pre`. + { + let entries_by_name: HashMap<&str, &Vec> = per_file_entries + .iter() + .map(|(name, entries)| (name.as_str(), entries)) + .collect(); + for ((file_name, idx), _action) in reconciliation_actions_pre.iter() { + if let Some(entries) = entries_by_name.get(file_name.as_str()) { + if let Some(e) = entries.get(*idx) { + first_pass_base_ids.insert(e.entry_id & 0x7FFF_FFFF); + } + } + } + } + // Save (file, entry_id) → action before per_file_entries shrinks // so we can rebuild (file, new_vec_idx) → action below. Use a // file_name → entries lookup map so the per-action lookup is O(1) @@ -519,8 +588,8 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( } log::info!( - "Worker compaction: {} -> {} entries ({} survived peptide-FDR or protein-rescue), \ - {} reconciliation actions retained", + "Worker compaction: {} -> {} entries ({} survived peptide-FDR, protein-rescue, \ + or cross-file consensus rescue), {} reconciliation actions retained", entries_before, entries_after, first_pass_base_ids.len(), @@ -570,6 +639,13 @@ pub fn run_rescore(config: OspreyConfig, library: Vec) -> Result<( &file_name_to_idx, &config, &mut seq_interner, + // The per-file rescore worker's config.input_files reflects + // only this worker's single parquet. To stamp the reconciled + // parquet with the multi-file reconciliation hash the + // downstream `--join-at-pass=2` merge node will validate + // against, pass the full join file_stems list extracted from + // reconciliation.json. + &join_file_stems, )?; // Cross-impl bisection seam: dump the per-precursor q-values