HPC chain Stage 7 correctness: 2nd-pass Percolator, file_stems-aware reconciliation hash, Stage 6 hard-error gates

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)

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<f64> = pairs.iter().map(|(x, _)| *x).collect();
         let y: Vec<f64> = pairs.iter().map(|(_, y)| *y).collect();

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<String> = 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<String> = 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())
     }

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<Vec<usize>> = 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<u32, usize> = HashMap::new();
+    let mut best_decoy: HashMap<u32, usize> = 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<usize> = 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<Vec<usize>> =
+        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<bool> = dedup_indices.iter().map(|&i| labels[i]).collect();
+            let dedup_entry_ids: Vec<u32> = dedup_indices.iter().map(|&i| entry_ids[i]).collect();
+            let dedup_peptides: Vec<String> =
+                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<f64> = 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<u32, f64> = 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, <feature_name_0>..<feature_name_N>`
+/// 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<usize> = (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

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 }