openmassspec_core/
centroid.rs

1//! Vendor-neutral profile-to-centroid peak picking.
2//!
3//! Centroiding is a transform over the arrays a [`SpectrumSource`] already
4//! yields (`mz`/`intensity`, and optionally `inv_mobility_per_peak`) - once a
5//! [`SpectrumRecord`] exists, the operation is identical regardless of which
6//! vendor parser produced it. [`Centroided`] wraps any `SpectrumSource` and
7//! applies the transform lazily, one spectrum at a time, so it composes with
8//! the streaming mzML writer and the Arrow bridge without buffering a whole
9//! run in memory.
10
11use crate::enums::ScanMode;
12use crate::source::SpectrumSource;
13use crate::types::{ChromatogramRecord, RunMetadata, SpectrumRecord};
14
15/// Wraps a [`SpectrumSource`], centroiding every profile-mode spectrum it
16/// yields. Spectra already tagged [`ScanMode::Centroid`] pass through
17/// unchanged (idempotent).
18pub struct Centroided<S: SpectrumSource> {
19    inner: S,
20    min_intensity: f32,
21}
22
23impl<S: SpectrumSource> Centroided<S> {
24    /// Wrap `inner`, picking every local-maximum peak regardless of height.
25    pub fn new(inner: S) -> Self {
26        Self {
27            inner,
28            min_intensity: 0.0,
29        }
30    }
31
32    /// Discard picked peaks below `min_intensity` (a simple noise floor).
33    pub fn with_min_intensity(mut self, min_intensity: f32) -> Self {
34        self.min_intensity = min_intensity;
35        self
36    }
37}
38
39impl<S: SpectrumSource> SpectrumSource for Centroided<S> {
40    fn run_metadata(&self) -> RunMetadata {
41        self.inner.run_metadata()
42    }
43
44    fn iter_spectra<'a>(&'a mut self) -> Box<dyn Iterator<Item = SpectrumRecord> + 'a> {
45        let min_intensity = self.min_intensity;
46        Box::new(
47            self.inner
48                .iter_spectra()
49                .map(move |rec| centroid_record(rec, min_intensity)),
50        )
51    }
52
53    fn iter_chromatograms<'a>(&'a mut self) -> Box<dyn Iterator<Item = ChromatogramRecord> + 'a> {
54        self.inner.iter_chromatograms()
55    }
56
57    fn spectrum_count_hint(&self) -> Option<usize> {
58        self.inner.spectrum_count_hint()
59    }
60
61    fn additional_processing_steps(&self) -> Vec<(&'static str, &'static str)> {
62        let mut steps = self.inner.additional_processing_steps();
63        steps.push(("MS:1000035", "peak picking"));
64        steps
65    }
66}
67
68fn centroid_record(mut rec: SpectrumRecord, min_intensity: f32) -> SpectrumRecord {
69    if rec.scan_mode == Some(ScanMode::Centroid) {
70        return rec;
71    }
72
73    let (mz, intensity, inv_mobility_per_peak) = pick_peaks(
74        &rec.mz,
75        &rec.intensity,
76        rec.inv_mobility_per_peak.as_deref(),
77        min_intensity,
78    );
79
80    rec.mz = mz;
81    rec.intensity = intensity;
82    rec.inv_mobility_per_peak = inv_mobility_per_peak;
83    rec.scan_mode = Some(ScanMode::Centroid);
84    // These were derived from the profile arrays; let the mzML writer /
85    // Arrow bridge recompute them from the new centroided arrays instead of
86    // carrying stale values forward (see `SpectrumRecord`'s field docs).
87    rec.total_ion_current = None;
88    rec.base_peak_mz = None;
89    rec.base_peak_intensity = None;
90    rec.low_mz = None;
91    rec.high_mz = None;
92    rec
93}
94
95/// Local-maxima peak picking: a point is a picked peak if it is no smaller
96/// than both neighbors and strictly larger than at least one of them (this
97/// also correctly picks a single-point spectrum, `n == 1`, since the "beats
98/// a neighbor" check is vacuously satisfied). The picked m/z (and, when
99/// present, inverse mobility) is the intensity-weighted centroid over the
100/// apex and its immediate neighbors; the picked intensity is the apex
101/// height. This is intentionally simple - see `timsrust-centroid` /
102/// `pyteomics` for more sophisticated peer approaches if this ever needs to
103/// improve on plain local-maxima picking.
104fn pick_peaks(
105    mz: &[f64],
106    intensity: &[f32],
107    inv_mobility_per_peak: Option<&[f32]>,
108    min_intensity: f32,
109) -> (Vec<f64>, Vec<f32>, Option<Vec<f32>>) {
110    let n = mz.len();
111    let mut out_mz = Vec::new();
112    let mut out_intensity = Vec::new();
113    let mut out_im = inv_mobility_per_peak.map(|_| Vec::new());
114
115    for i in 0..n {
116        let no_smaller_than_left = i == 0 || intensity[i] >= intensity[i - 1];
117        let no_smaller_than_right = i == n - 1 || intensity[i] >= intensity[i + 1];
118        let beats_a_neighbor = n == 1
119            || (i > 0 && intensity[i] > intensity[i - 1])
120            || (i < n - 1 && intensity[i] > intensity[i + 1]);
121        if !(no_smaller_than_left && no_smaller_than_right && beats_a_neighbor) {
122            continue;
123        }
124        if intensity[i] < min_intensity {
125            continue;
126        }
127
128        let lo = i.saturating_sub(1);
129        let hi = (i + 1).min(n - 1);
130        let mut weighted_mz = 0.0f64;
131        let mut weighted_im = 0.0f64;
132        let mut weight_sum = 0.0f64;
133        for j in lo..=hi {
134            let w = f64::from(intensity[j]);
135            weighted_mz += mz[j] * w;
136            weight_sum += w;
137            if let Some(im) = inv_mobility_per_peak {
138                weighted_im += f64::from(im[j]) * w;
139            }
140        }
141        let centroid_mz = if weight_sum > 0.0 {
142            weighted_mz / weight_sum
143        } else {
144            mz[i]
145        };
146        out_mz.push(centroid_mz);
147        out_intensity.push(intensity[i]);
148        if let Some(out) = out_im.as_mut() {
149            #[allow(clippy::cast_possible_truncation)]
150            let centroid_im = if weight_sum > 0.0 {
151                (weighted_im / weight_sum) as f32
152            } else {
153                inv_mobility_per_peak.expect("out_im is Some only when inv_mobility_per_peak is")[i]
154            };
155            out.push(centroid_im);
156        }
157    }
158
159    (out_mz, out_intensity, out_im)
160}
161
162#[cfg(test)]
163mod tests {
164    use super::*;
165    use crate::conformance::assert_source_invariants;
166    use crate::{CvTerm, MsPower, Polarity};
167
168    struct OneSpectrumSource {
169        meta: RunMetadata,
170        rec: Option<SpectrumRecord>,
171    }
172
173    impl OneSpectrumSource {
174        fn new(rec: SpectrumRecord) -> Self {
175            Self {
176                meta: RunMetadata {
177                    source_file_name: "toy.raw".into(),
178                    source_file_format: CvTerm::new("MS:1000563", "Thermo RAW format"),
179                    native_id_format: CvTerm::new("MS:1000768", "Thermo nativeID format"),
180                    instrument: CvTerm::new("MS:1001911", "Q Exactive"),
181                    instrument_serial_number: None,
182                    software_name: "toy-writer".into(),
183                    software_version: "0.0.0".into(),
184                    acquisition_software_name: None,
185                    acquisition_software_version: None,
186                    start_timestamp: None,
187                    mobility_array_kind: None,
188                    analyzers: Vec::new(),
189                    extra: Default::default(),
190                },
191                rec: Some(rec),
192            }
193        }
194    }
195
196    impl SpectrumSource for OneSpectrumSource {
197        fn run_metadata(&self) -> RunMetadata {
198            self.meta.clone()
199        }
200
201        fn iter_spectra<'a>(&'a mut self) -> Box<dyn Iterator<Item = SpectrumRecord> + 'a> {
202            Box::new(self.rec.take().into_iter())
203        }
204
205        fn spectrum_count_hint(&self) -> Option<usize> {
206            Some(usize::from(self.rec.is_some()))
207        }
208    }
209
210    fn profile_spectrum() -> SpectrumRecord {
211        // Two synthetic Gaussian-ish bumps centered at mz=100 and mz=200,
212        // known peak apexes so we can assert the picked centroids land near
213        // them. Compared only against this hand-built synthetic input, not
214        // any vendor tool or reference converter (clean-room rule).
215        let mz: Vec<f64> = vec![
216            99.0, 99.5, 100.0, 100.5, 101.0, // bump 1, apex at 100.0
217            199.0, 199.5, 200.0, 200.5, 201.0, // bump 2, apex at 200.0
218        ];
219        let intensity: Vec<f32> = vec![10.0, 50.0, 100.0, 50.0, 10.0, 5.0, 40.0, 80.0, 40.0, 5.0];
220        SpectrumRecord {
221            index: 0,
222            scan_number: 1,
223            native_id: "scan=1".into(),
224            ms_level: MsPower::Ms1.ms_level(),
225            polarity: Some(Polarity::Positive),
226            scan_mode: Some(ScanMode::Profile),
227            analyzer: None,
228            acquisition_event_id: None,
229            filter: None,
230            retention_time_sec: 1.0,
231            total_ion_current: Some(intensity.iter().map(|&v| v as f64).sum()),
232            base_peak_mz: Some(100.0),
233            base_peak_intensity: Some(100.0),
234            low_mz: Some(99.0),
235            high_mz: Some(201.0),
236            ion_injection_time_ms: None,
237            inv_mobility: None,
238            faims_cv: None,
239            precursor: None,
240            mz,
241            intensity,
242            inv_mobility_per_peak: None,
243            extra: Default::default(),
244        }
245    }
246
247    #[test]
248    fn centroids_profile_spectrum_at_known_apexes() {
249        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
250        let recs: Vec<_> = src.iter_spectra().collect();
251        assert_eq!(recs.len(), 1);
252        let rec = &recs[0];
253
254        assert_eq!(rec.scan_mode, Some(ScanMode::Centroid));
255        assert!(
256            rec.mz.len() < 10,
257            "expected fewer peaks after centroiding, got {}",
258            rec.mz.len()
259        );
260        assert_eq!(rec.mz.len(), rec.intensity.len());
261
262        let near_100 = rec.mz.iter().any(|&m| (m - 100.0).abs() < 0.6);
263        let near_200 = rec.mz.iter().any(|&m| (m - 200.0).abs() < 0.6);
264        assert!(near_100, "no picked peak near mz=100: {:?}", rec.mz);
265        assert!(near_200, "no picked peak near mz=200: {:?}", rec.mz);
266
267        // Stale profile-derived summary fields must be cleared, not carried
268        // forward, so the writer recomputes them from the new arrays.
269        assert_eq!(rec.total_ion_current, None);
270        assert_eq!(rec.base_peak_mz, None);
271    }
272
273    #[test]
274    fn already_centroided_spectrum_passes_through_unchanged() {
275        let mut rec = profile_spectrum();
276        rec.scan_mode = Some(ScanMode::Centroid);
277        let original_mz = rec.mz.clone();
278
279        let mut src = Centroided::new(OneSpectrumSource::new(rec));
280        let recs: Vec<_> = src.iter_spectra().collect();
281        assert_eq!(recs[0].mz, original_mz);
282    }
283
284    #[test]
285    fn additional_processing_steps_reports_peak_picking() {
286        let src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
287        assert_eq!(
288            src.additional_processing_steps(),
289            vec![("MS:1000035", "peak picking")]
290        );
291    }
292
293    #[test]
294    fn additional_processing_steps_accumulate_through_nested_adapters() {
295        // Wrapping a Centroided source in another Centroided (a no-op in
296        // practice, since it is idempotent) must not drop the inner
297        // adapter's step - steps accumulate outer-to-inner delegation.
298        let inner = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
299        let outer = Centroided::new(inner);
300        assert_eq!(
301            outer.additional_processing_steps(),
302            vec![
303                ("MS:1000035", "peak picking"),
304                ("MS:1000035", "peak picking")
305            ]
306        );
307    }
308
309    #[test]
310    fn wrapped_source_still_satisfies_conformance_invariants() {
311        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
312        let n = assert_source_invariants(&mut src).expect("conformance");
313        assert_eq!(n, 1);
314    }
315
316    #[test]
317    fn composes_with_the_streaming_mzml_writer() {
318        // Centroided is just another SpectrumSource, so it drops straight
319        // into write_mzml with no special-casing on the writer's part.
320        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
321        let mut buf = Vec::new();
322        crate::write_mzml(&mut src, &mut buf).expect("write_mzml");
323        let xml = String::from_utf8(buf).expect("utf8");
324
325        assert!(xml.contains(r#"<spectrumList count="1""#));
326        assert!(xml.contains(
327            r#"<cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum" value=""/>"#
328        ));
329        assert!(
330            !xml.contains(
331                r#"<cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum" value=""/>"#
332            ),
333            "output should not still claim profile mode after centroiding"
334        );
335        // The dataProcessingList must record that peak picking happened,
336        // not just the blanket "Conversion to mzML" step every writer run
337        // gets (see issue: centroiding left no provenance trail).
338        assert!(xml.contains(
339            r#"<cvParam cvRef="MS" accession="MS:1000035" name="peak picking" value=""/>"#
340        ));
341    }
342}
openmassspec_core/centroid.rs

openmassspec_core/
centroid.rs
