thermo_parity/
thermo_parity.rs

1//! Compare the native Thermo `.raw` reader against an mzML conversion of the
2//! same run (msconvert with vendor peak picking).
3//!
4//! ```text
5//! cargo run --release --example thermo_parity -- ms1  RUN.raw RUN.mzML[.gz]
6//! cargo run --release --example thermo_parity -- ms2  RUN.raw RUN.mzML[.gz]
7//! cargo run --release --example thermo_parity -- load RUN.raw|RUN.mzML
8//! ```
9//!
10//! `ms1` pairs MS1 spectra in RT order and reports scan counts, RT, per-scan
11//! peak counts and m/z / intensity agreement. `ms2` compares the DIA isolation
12//! windows each reader recovers. `load` reads MS1 once and reports the time.
13
14use std::collections::BTreeMap;
15use std::path::Path;
16use std::time::Instant;
17
18use koth_ms::config::FileConfig;
19use koth_ms::io;
20use koth_ms::models::Spectrum;
21
22fn main() -> anyhow::Result<()> {
23    let args: Vec<String> = std::env::args().collect();
24    match args.get(1).map(String::as_str) {
25        Some("ms1") => ms1(Path::new(&args[2]), Path::new(&args[3])),
26        Some("ms2") => ms2(Path::new(&args[2]), Path::new(&args[3])),
27        Some("load") => load(Path::new(&args[2])),
28        _ => anyhow::bail!("usage: thermo_parity ms1|ms2 RAW MZML | load PATH"),
29    }
30}
31
32fn load(path: &Path) -> anyhow::Result<()> {
33    let t = Instant::now();
34    let mut n = 0usize;
35    let mut peaks = 0usize;
36    for s in io::stream_spectra(path, &FileConfig::default())? {
37        let s = s?;
38        n += 1;
39        peaks += s.peaks.len();
40    }
41    println!(
42        "load {}: {n} MS1 spectra, {peaks} peaks, {:.2} s",
43        path.display(),
44        t.elapsed().as_secs_f64()
45    );
46    Ok(())
47}
48
49fn read_timed(path: &Path) -> anyhow::Result<Vec<Spectrum>> {
50    let t = Instant::now();
51    let spectra = io::read_spectra(path, &FileConfig::default())?;
52    println!(
53        "read {:>6} MS1 spectra in {:6.2} s from {}",
54        spectra.len(),
55        t.elapsed().as_secs_f64(),
56        path.display()
57    );
58    Ok(spectra)
59}
60
61fn ms1(raw: &Path, mzml: &Path) -> anyhow::Result<()> {
62    let a = read_timed(raw)?;
63    let b = read_timed(mzml)?;
64    println!("MS1 spectra: raw {} mzML {}", a.len(), b.len());
65    let n = a.len().min(b.len());
66    let mut max_rt = 0.0f64;
67    let mut same_count = 0usize;
68    let (mut peaks_a, mut peaks_b) = (0usize, 0usize);
69    let mut max_ppm = 0.0f64;
70    let mut max_rel_int = 0.0f64;
71    let mut compared = 0usize;
72    let mut faims_mismatch = 0usize;
73    for (sa, sb) in a.iter().zip(b.iter()).take(n) {
74        max_rt = max_rt.max((sa.retention_time - sb.retention_time).abs());
75        peaks_a += sa.peaks.len();
76        peaks_b += sb.peaks.len();
77        if sa.faims_cv != sb.faims_cv {
78            faims_mismatch += 1;
79        }
80        if sa.peaks.len() != sb.peaks.len() {
81            continue;
82        }
83        same_count += 1;
84        for (pa, pb) in sa.peaks.iter().zip(&sb.peaks) {
85            let ppm = ((pa.mz - pb.mz) as f64).abs() / pb.mz as f64 * 1e6;
86            max_ppm = max_ppm.max(ppm);
87            let rel =
88                ((pa.intensity - pb.intensity) as f64).abs() / (pb.intensity as f64).max(1e-12);
89            max_rel_int = max_rel_int.max(rel);
90            compared += 1;
91        }
92    }
93    println!("paired {n}; max |dRT| = {:.3e} min", max_rt);
94    println!(
95        "identical peak count in {same_count}/{n} scans; total peaks raw {peaks_a} mzML {peaks_b}"
96    );
97    println!(
98        "over {compared} paired peaks: max |dm/z| = {max_ppm:.4} ppm, max rel |dI| = {max_rel_int:.3e}"
99    );
100    println!("FAIMS CV mismatches: {faims_mismatch}");
101    Ok(())
102}
103
104fn ms2(raw: &Path, mzml: &Path) -> anyhow::Result<()> {
105    let t = Instant::now();
106    let a = io::thermo::read_thermo_ms2(raw)?;
107    println!(
108        "raw: {} MS2 spectra in {:.2} s",
109        a.len(),
110        t.elapsed().as_secs_f64()
111    );
112    let b: Vec<Spectrum> = io::mzml::stream_mzml_ms2(mzml)?.collect();
113    println!("mzML: {} MS2 spectra", b.len());
114    let windows = |v: &[Spectrum]| {
115        let mut m: BTreeMap<(i64, i64, i64), usize> = BTreeMap::new();
116        for s in v {
117            *m.entry(s.isolation_window.unwrap().key()).or_default() += 1;
118        }
119        m
120    };
121    let (wa, wb) = (windows(&a), windows(&b));
122    println!("distinct windows: raw {} mzML {}", wa.len(), wb.len());
123    let fmt = |k: &(i64, i64, i64)| {
124        format!(
125            "{:.3} [{:.3}, {:.3}]",
126            k.0 as f64 / 1e4,
127            k.1 as f64 / 1e4,
128            k.2 as f64 / 1e4
129        )
130    };
131    for (k, c) in wa.iter().take(5) {
132        println!("  raw  {} x{c}", fmt(k));
133    }
134    for (k, c) in wb.iter().take(5) {
135        println!("  mzML {} x{c}", fmt(k));
136    }
137    // Pair scans in RT order and compare windows and peaks.
138    let n = a.len().min(b.len());
139    let (mut same_w, mut same_n, mut max_rt, mut max_bound) = (0usize, 0usize, 0.0f64, 0.0f64);
140    for (sa, sb) in a.iter().zip(&b).take(n) {
141        let (wa, wb) = (sa.isolation_window.unwrap(), sb.isolation_window.unwrap());
142        max_rt = max_rt.max((sa.retention_time - sb.retention_time).abs());
143        let d = (wa.target - wb.target)
144            .abs()
145            .max((wa.lower - wb.lower).abs())
146            .max((wa.upper - wb.upper).abs());
147        max_bound = max_bound.max(d);
148        if wa.key() == wb.key() {
149            same_w += 1;
150        }
151        if sa.peaks.len() == sb.peaks.len() {
152            same_n += 1;
153        }
154    }
155    println!(
156        "paired {n}: identical window key {same_w}, identical peak count {same_n}, \
157         max |dRT| {max_rt:.3e} min, max window bound diff {max_bound:.4} Th"
158    );
159    Ok(())
160}
thermo_parity/thermo_parity.rs

thermo_parity/
thermo_parity.rs
