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
62fn centroid_record(mut rec: SpectrumRecord, min_intensity: f32) -> SpectrumRecord {
63    if rec.scan_mode == Some(ScanMode::Centroid) {
64        return rec;
65    }
66
67    let (mz, intensity, inv_mobility_per_peak) = pick_peaks(
68        &rec.mz,
69        &rec.intensity,
70        rec.inv_mobility_per_peak.as_deref(),
71        min_intensity,
72    );
73
74    rec.mz = mz;
75    rec.intensity = intensity;
76    rec.inv_mobility_per_peak = inv_mobility_per_peak;
77    rec.scan_mode = Some(ScanMode::Centroid);
78    // These were derived from the profile arrays; let the mzML writer /
79    // Arrow bridge recompute them from the new centroided arrays instead of
80    // carrying stale values forward (see `SpectrumRecord`'s field docs).
81    rec.total_ion_current = None;
82    rec.base_peak_mz = None;
83    rec.base_peak_intensity = None;
84    rec.low_mz = None;
85    rec.high_mz = None;
86    rec
87}
88
89/// Local-maxima peak picking: a point is a picked peak if it is no smaller
90/// than both neighbors and strictly larger than at least one of them (this
91/// also correctly picks a single-point spectrum, `n == 1`, since the "beats
92/// a neighbor" check is vacuously satisfied). The picked m/z (and, when
93/// present, inverse mobility) is the intensity-weighted centroid over the
94/// apex and its immediate neighbors; the picked intensity is the apex
95/// height. This is intentionally simple - see `timsrust-centroid` /
96/// `pyteomics` for more sophisticated peer approaches if this ever needs to
97/// improve on plain local-maxima picking.
98fn pick_peaks(
99    mz: &[f64],
100    intensity: &[f32],
101    inv_mobility_per_peak: Option<&[f32]>,
102    min_intensity: f32,
103) -> (Vec<f64>, Vec<f32>, Option<Vec<f32>>) {
104    let n = mz.len();
105    let mut out_mz = Vec::new();
106    let mut out_intensity = Vec::new();
107    let mut out_im = inv_mobility_per_peak.map(|_| Vec::new());
108
109    for i in 0..n {
110        let no_smaller_than_left = i == 0 || intensity[i] >= intensity[i - 1];
111        let no_smaller_than_right = i == n - 1 || intensity[i] >= intensity[i + 1];
112        let beats_a_neighbor = n == 1
113            || (i > 0 && intensity[i] > intensity[i - 1])
114            || (i < n - 1 && intensity[i] > intensity[i + 1]);
115        if !(no_smaller_than_left && no_smaller_than_right && beats_a_neighbor) {
116            continue;
117        }
118        if intensity[i] < min_intensity {
119            continue;
120        }
121
122        let lo = i.saturating_sub(1);
123        let hi = (i + 1).min(n - 1);
124        let mut weighted_mz = 0.0f64;
125        let mut weighted_im = 0.0f64;
126        let mut weight_sum = 0.0f64;
127        for j in lo..=hi {
128            let w = f64::from(intensity[j]);
129            weighted_mz += mz[j] * w;
130            weight_sum += w;
131            if let Some(im) = inv_mobility_per_peak {
132                weighted_im += f64::from(im[j]) * w;
133            }
134        }
135        let centroid_mz = if weight_sum > 0.0 {
136            weighted_mz / weight_sum
137        } else {
138            mz[i]
139        };
140        out_mz.push(centroid_mz);
141        out_intensity.push(intensity[i]);
142        if let Some(out) = out_im.as_mut() {
143            #[allow(clippy::cast_possible_truncation)]
144            let centroid_im = if weight_sum > 0.0 {
145                (weighted_im / weight_sum) as f32
146            } else {
147                inv_mobility_per_peak.expect("out_im is Some only when inv_mobility_per_peak is")[i]
148            };
149            out.push(centroid_im);
150        }
151    }
152
153    (out_mz, out_intensity, out_im)
154}
155
156#[cfg(test)]
157mod tests {
158    use super::*;
159    use crate::conformance::assert_source_invariants;
160    use crate::{CvTerm, MsPower, Polarity};
161
162    struct OneSpectrumSource {
163        meta: RunMetadata,
164        rec: Option<SpectrumRecord>,
165    }
166
167    impl OneSpectrumSource {
168        fn new(rec: SpectrumRecord) -> Self {
169            Self {
170                meta: RunMetadata {
171                    source_file_name: "toy.raw".into(),
172                    source_file_format: CvTerm::new("MS:1000563", "Thermo RAW format"),
173                    native_id_format: CvTerm::new("MS:1000768", "Thermo nativeID format"),
174                    instrument: CvTerm::new("MS:1001911", "Q Exactive"),
175                    software_name: "toy-writer".into(),
176                    software_version: "0.0.0".into(),
177                    start_timestamp: None,
178                    mobility_array_kind: None,
179                },
180                rec: Some(rec),
181            }
182        }
183    }
184
185    impl SpectrumSource for OneSpectrumSource {
186        fn run_metadata(&self) -> RunMetadata {
187            self.meta.clone()
188        }
189
190        fn iter_spectra<'a>(&'a mut self) -> Box<dyn Iterator<Item = SpectrumRecord> + 'a> {
191            Box::new(self.rec.take().into_iter())
192        }
193
194        fn spectrum_count_hint(&self) -> Option<usize> {
195            Some(usize::from(self.rec.is_some()))
196        }
197    }
198
199    fn profile_spectrum() -> SpectrumRecord {
200        // Two synthetic Gaussian-ish bumps centered at mz=100 and mz=200,
201        // known peak apexes so we can assert the picked centroids land near
202        // them. Compared only against this hand-built synthetic input, not
203        // any vendor tool or reference converter (clean-room rule).
204        let mz: Vec<f64> = vec![
205            99.0, 99.5, 100.0, 100.5, 101.0, // bump 1, apex at 100.0
206            199.0, 199.5, 200.0, 200.5, 201.0, // bump 2, apex at 200.0
207        ];
208        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];
209        SpectrumRecord {
210            index: 0,
211            scan_number: 1,
212            native_id: "scan=1".into(),
213            ms_level: MsPower::Ms1.ms_level(),
214            polarity: Some(Polarity::Positive),
215            scan_mode: Some(ScanMode::Profile),
216            analyzer: None,
217            filter: None,
218            retention_time_sec: 1.0,
219            total_ion_current: Some(intensity.iter().map(|&v| v as f64).sum()),
220            base_peak_mz: Some(100.0),
221            base_peak_intensity: Some(100.0),
222            low_mz: Some(99.0),
223            high_mz: Some(201.0),
224            ion_injection_time_ms: None,
225            inv_mobility: None,
226            precursor: None,
227            mz,
228            intensity,
229            inv_mobility_per_peak: None,
230        }
231    }
232
233    #[test]
234    fn centroids_profile_spectrum_at_known_apexes() {
235        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
236        let recs: Vec<_> = src.iter_spectra().collect();
237        assert_eq!(recs.len(), 1);
238        let rec = &recs[0];
239
240        assert_eq!(rec.scan_mode, Some(ScanMode::Centroid));
241        assert!(
242            rec.mz.len() < 10,
243            "expected fewer peaks after centroiding, got {}",
244            rec.mz.len()
245        );
246        assert_eq!(rec.mz.len(), rec.intensity.len());
247
248        let near_100 = rec.mz.iter().any(|&m| (m - 100.0).abs() < 0.6);
249        let near_200 = rec.mz.iter().any(|&m| (m - 200.0).abs() < 0.6);
250        assert!(near_100, "no picked peak near mz=100: {:?}", rec.mz);
251        assert!(near_200, "no picked peak near mz=200: {:?}", rec.mz);
252
253        // Stale profile-derived summary fields must be cleared, not carried
254        // forward, so the writer recomputes them from the new arrays.
255        assert_eq!(rec.total_ion_current, None);
256        assert_eq!(rec.base_peak_mz, None);
257    }
258
259    #[test]
260    fn already_centroided_spectrum_passes_through_unchanged() {
261        let mut rec = profile_spectrum();
262        rec.scan_mode = Some(ScanMode::Centroid);
263        let original_mz = rec.mz.clone();
264
265        let mut src = Centroided::new(OneSpectrumSource::new(rec));
266        let recs: Vec<_> = src.iter_spectra().collect();
267        assert_eq!(recs[0].mz, original_mz);
268    }
269
270    #[test]
271    fn wrapped_source_still_satisfies_conformance_invariants() {
272        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
273        let n = assert_source_invariants(&mut src).expect("conformance");
274        assert_eq!(n, 1);
275    }
276
277    #[test]
278    fn composes_with_the_streaming_mzml_writer() {
279        // Centroided is just another SpectrumSource, so it drops straight
280        // into write_mzml with no special-casing on the writer's part.
281        let mut src = Centroided::new(OneSpectrumSource::new(profile_spectrum()));
282        let mut buf = Vec::new();
283        crate::write_mzml(&mut src, &mut buf).expect("write_mzml");
284        let xml = String::from_utf8(buf).expect("utf8");
285
286        assert!(xml.contains(r#"<spectrumList count="1""#));
287        assert!(xml.contains(
288            r#"<cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum" value=""/>"#
289        ));
290        assert!(
291            !xml.contains(
292                r#"<cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum" value=""/>"#
293            ),
294            "output should not still claim profile mode after centroiding"
295        );
296    }
297}
openmassspec_core/centroid.rs

openmassspec_core/
centroid.rs
