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                },
190                rec: Some(rec),
191            }
192        }
193    }
194
195    impl SpectrumSource for OneSpectrumSource {
196        fn run_metadata(&self) -> RunMetadata {
197            self.meta.clone()
198        }
199
200        fn iter_spectra<'a>(&'a mut self) -> Box<dyn Iterator<Item = SpectrumRecord> + 'a> {
201            Box::new(self.rec.take().into_iter())
202        }
203
204        fn spectrum_count_hint(&self) -> Option<usize> {
205            Some(usize::from(self.rec.is_some()))
206        }
207    }
208
209    fn profile_spectrum() -> SpectrumRecord {
210        // Two synthetic Gaussian-ish bumps centered at mz=100 and mz=200,
211        // known peak apexes so we can assert the picked centroids land near
212        // them. Compared only against this hand-built synthetic input, not
213        // any vendor tool or reference converter (clean-room rule).
214        let mz: Vec<f64> = vec![
215            99.0, 99.5, 100.0, 100.5, 101.0, // bump 1, apex at 100.0
216            199.0, 199.5, 200.0, 200.5, 201.0, // bump 2, apex at 200.0
217        ];
218        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];
219        SpectrumRecord {
220            index: 0,
221            scan_number: 1,
222            native_id: "scan=1".into(),
223            ms_level: MsPower::Ms1.ms_level(),
224            polarity: Some(Polarity::Positive),
225            scan_mode: Some(ScanMode::Profile),
226            analyzer: None,
227            filter: None,
228            retention_time_sec: 1.0,
229            total_ion_current: Some(intensity.iter().map(|&v| v as f64).sum()),
230            base_peak_mz: Some(100.0),
231            base_peak_intensity: Some(100.0),
232            low_mz: Some(99.0),
233            high_mz: Some(201.0),
234            ion_injection_time_ms: None,
235            inv_mobility: None,
236            faims_cv: None,
237            precursor: None,
238            mz,
239            intensity,
240            inv_mobility_per_peak: None,
241        }
242    }
243
244    #[test]
245    fn centroids_profile_spectrum_at_known_apexes() {
246        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
247        let recs: Vec<_> = src.iter_spectra().collect();
248        assert_eq!(recs.len(), 1);
249        let rec = &recs[0];
250
251        assert_eq!(rec.scan_mode, Some(ScanMode::Centroid));
252        assert!(
253            rec.mz.len() < 10,
254            "expected fewer peaks after centroiding, got {}",
255            rec.mz.len()
256        );
257        assert_eq!(rec.mz.len(), rec.intensity.len());
258
259        let near_100 = rec.mz.iter().any(|&m| (m - 100.0).abs() < 0.6);
260        let near_200 = rec.mz.iter().any(|&m| (m - 200.0).abs() < 0.6);
261        assert!(near_100, "no picked peak near mz=100: {:?}", rec.mz);
262        assert!(near_200, "no picked peak near mz=200: {:?}", rec.mz);
263
264        // Stale profile-derived summary fields must be cleared, not carried
265        // forward, so the writer recomputes them from the new arrays.
266        assert_eq!(rec.total_ion_current, None);
267        assert_eq!(rec.base_peak_mz, None);
268    }
269
270    #[test]
271    fn already_centroided_spectrum_passes_through_unchanged() {
272        let mut rec = profile_spectrum();
273        rec.scan_mode = Some(ScanMode::Centroid);
274        let original_mz = rec.mz.clone();
275
276        let mut src = Centroided::new(OneSpectrumSource::new(rec));
277        let recs: Vec<_> = src.iter_spectra().collect();
278        assert_eq!(recs[0].mz, original_mz);
279    }
280
281    #[test]
282    fn additional_processing_steps_reports_peak_picking() {
283        let src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
284        assert_eq!(
285            src.additional_processing_steps(),
286            vec![("MS:1000035", "peak picking")]
287        );
288    }
289
290    #[test]
291    fn additional_processing_steps_accumulate_through_nested_adapters() {
292        // Wrapping a Centroided source in another Centroided (a no-op in
293        // practice, since it is idempotent) must not drop the inner
294        // adapter's step - steps accumulate outer-to-inner delegation.
295        let inner = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
296        let outer = Centroided::new(inner);
297        assert_eq!(
298            outer.additional_processing_steps(),
299            vec![
300                ("MS:1000035", "peak picking"),
301                ("MS:1000035", "peak picking")
302            ]
303        );
304    }
305
306    #[test]
307    fn wrapped_source_still_satisfies_conformance_invariants() {
308        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
309        let n = assert_source_invariants(&mut src).expect("conformance");
310        assert_eq!(n, 1);
311    }
312
313    #[test]
314    fn composes_with_the_streaming_mzml_writer() {
315        // Centroided is just another SpectrumSource, so it drops straight
316        // into write_mzml with no special-casing on the writer's part.
317        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
318        let mut buf = Vec::new();
319        crate::write_mzml(&mut src, &mut buf).expect("write_mzml");
320        let xml = String::from_utf8(buf).expect("utf8");
321
322        assert!(xml.contains(r#"<spectrumList count="1""#));
323        assert!(xml.contains(
324            r#"<cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum" value=""/>"#
325        ));
326        assert!(
327            !xml.contains(
328                r#"<cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum" value=""/>"#
329            ),
330            "output should not still claim profile mode after centroiding"
331        );
332        // The dataProcessingList must record that peak picking happened,
333        // not just the blanket "Conversion to mzML" step every writer run
334        // gets (see issue: centroiding left no provenance trail).
335        assert!(xml.contains(
336            r#"<cvParam cvRef="MS" accession="MS:1000035" name="peak picking" value=""/>"#
337        ));
338    }
339}
openmassspec_core/centroid.rs

openmassspec_core/
centroid.rs
