rustms/ms/
spectrum.rs

1use std::collections::BTreeMap;
2use serde::{Deserialize, Serialize};
3
4/// Represents a mass spectrum with associated m/z values and intensities.
5#[derive(Clone, Debug, Serialize, Deserialize)]
6pub struct MzSpectrum {
7    pub mz: Vec<f64>,
8    pub intensity: Vec<f64>,
9}
10
11impl MzSpectrum {
12    /// Constructs a new `MzSpectrum`.
13    ///
14    /// # Arguments
15    ///
16    /// * `mz` - A vector of m/z values.
17    /// * `intensity` - A vector of intensity values corresponding to the m/z values.
18    ///
19    /// # Panics
20    ///
21    /// Panics if the lengths of `mz` and `intensity` are not the same. (actually, it doesn't at the moment, planning on adding this later)
22    ///
23    /// # Example
24    ///
25    /// ```rust
26    /// # use rustms::ms::spectrum::MzSpectrum;
27    /// let spectrum = MzSpectrum::new(vec![200.0, 100.0], vec![20.0, 10.0]);
28    /// assert_eq!(spectrum.mz, vec![100.0, 200.0]);
29    /// assert_eq!(spectrum.intensity, vec![10.0, 20.0]);
30    /// ```
31    pub fn new(mz: Vec<f64>, intensity: Vec<f64>) -> Self {
32        assert_eq!(mz.len(), intensity.len(), "mz and intensity vectors must have the same length");
33        // make sure mz and intensity are sorted by mz
34        let mut mz_intensity: Vec<(f64, f64)> = mz.iter().zip(intensity.iter()).map(|(m, i)| (*m, *i)).collect();
35        mz_intensity.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
36        MzSpectrum { mz: mz_intensity.iter().map(|(m, _)| *m).collect(), intensity: mz_intensity.iter().map(|(_, i)| *i).collect() }
37    }
38
39    /// Filters the m/z values and intensities based on a range of m/z values and intensities.
40    ///
41    /// # Arguments
42    ///
43    /// * `mz_min` - The minimum m/z value.
44    /// * `mz_max` - The maximum m/z value.
45    /// * `intensity_min` - The minimum intensity value.
46    /// * `intensity_max` - The maximum intensity value.
47    ///
48    /// # Returns
49    ///
50    /// * `MzSpectrum` - A new `MzSpectrum` with m/z values and intensities within the specified ranges.
51    ///
52    /// # Example
53    ///
54    /// ```rust
55    /// # use rustms::ms::spectrum::MzSpectrum;
56    /// let spectrum = MzSpectrum::new(vec![100.0, 200.0, 300.0], vec![10.0, 20.0, 30.0]);
57    /// let filtered_spectrum = spectrum.filter_ranged(150.0, 250.0, 15.0, 25.0);
58    /// assert_eq!(filtered_spectrum.mz, vec![200.0]);
59    /// assert_eq!(filtered_spectrum.intensity, vec![20.0]);
60    /// ```
61    pub fn filter_ranged(&self, mz_min: f64, mz_max: f64, intensity_min: f64, intensity_max: f64) -> Self {
62        let mut mz_vec: Vec<f64> = Vec::new();
63        let mut intensity_vec: Vec<f64> = Vec::new();
64
65        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
66            if mz_min <= *mz && *mz <= mz_max && *intensity >= intensity_min && *intensity <= intensity_max {
67                mz_vec.push(*mz);
68                intensity_vec.push(*intensity);
69            }
70        }
71        MzSpectrum { mz: mz_vec, intensity: intensity_vec }
72    }
73
74    pub fn from_collection(collection: Vec<MzSpectrum>) -> MzSpectrum {
75
76        let quantize = |mz: f64| -> i64 {
77            (mz * 1_000_000.0).round() as i64
78        };
79
80        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
81
82        for spectrum in collection {
83            for (mz, intensity) in spectrum.mz.iter().zip(spectrum.intensity.iter()) {
84                let key = quantize(*mz);
85                let entry = combined_map.entry(key).or_insert(0.0);
86                *entry += *intensity;
87            }
88        }
89
90        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
91        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
92
93        MzSpectrum { mz: mz_combined, intensity: intensity_combined }
94    }
95}
96
97impl std::ops::Add for MzSpectrum {
98    type Output = Self;
99    /// Combines two `MzSpectrum` instances by summing up the intensities of matching m/z values.
100    ///
101    /// # Description
102    /// Each m/z value is quantized to retain at least 6 decimals. If two spectra have m/z values
103    /// that quantize to the same integer value, their intensities are summed.
104    ///
105    /// # Example
106    /// ```
107    /// # use rustms::ms::spectrum::MzSpectrum;
108    /// let spectrum1 = MzSpectrum { mz: vec![100.523, 101.923], intensity: vec![10.0, 20.0] };
109    /// let spectrum2 = MzSpectrum { mz: vec![101.235, 105.112], intensity: vec![15.0, 30.0] };
110    ///
111    /// let combined = spectrum1 + spectrum2;
112    ///
113    /// assert_eq!(combined.mz, vec![100.523, 101.235, 101.923, 105.112]);
114    /// assert_eq!(combined.intensity, vec![10.0, 15.0, 20.0, 30.0]);
115    /// ```
116    fn add(self, other: Self) -> MzSpectrum {
117        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
118
119        // Helper to quantize mz to an integer key
120        let quantize = |mz: f64| -> i64 {
121            (mz * 1_000_000.0).round() as i64
122        };
123
124        // Add the m/z and intensities from the first spectrum to the map
125        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
126            let key = quantize(*mz);
127            combined_map.insert(key, *intensity);
128        }
129
130        // Combine the second spectrum into the map
131        for (mz, intensity) in other.mz.iter().zip(other.intensity.iter()) {
132            let key = quantize(*mz);
133            let entry = combined_map.entry(key).or_insert(0.0);
134            *entry += *intensity;
135        }
136
137        // Convert the combined map back into two Vec<f64>
138        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
139        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
140
141        MzSpectrum { mz: mz_combined, intensity: intensity_combined }
142    }
143}
144
145impl std::ops::Mul<f64> for MzSpectrum {
146    type Output = Self;
147    fn mul(self, scale: f64) -> Self::Output{
148        let mut scaled_intensities: Vec<f64> = vec![0.0; self.intensity.len()];
149        for (idx,intensity) in self.intensity.iter().enumerate(){
150            scaled_intensities[idx] = scale*intensity;
151        }
152        Self{ mz: self.mz.clone(), intensity: scaled_intensities}
153
154    }
155}
156
157impl std::ops::Sub for MzSpectrum {
158    type Output = Self;
159    fn sub(self, other: Self) -> Self::Output {
160        let mut combined_map: BTreeMap<i64, f64> = BTreeMap::new();
161
162        // Helper to quantize mz to an integer key
163        let quantize = |mz: f64| -> i64 {
164            (mz * 1_000_000.0).round() as i64
165        };
166
167        // Add the m/z and intensities from the first spectrum to the map
168        for (mz, intensity) in self.mz.iter().zip(self.intensity.iter()) {
169            let key = quantize(*mz);
170            combined_map.insert(key, *intensity);
171        }
172
173        // Combine the second spectrum into the map
174        for (mz, intensity) in other.mz.iter().zip(other.intensity.iter()) {
175            let key = quantize(*mz);
176            let entry = combined_map.entry(key).or_insert(0.0);
177            *entry -= *intensity;
178        }
179
180        // Convert the combined map back into two Vec<f64>
181        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
182        let intensity_combined: Vec<f64> = combined_map.values().cloned().collect();
183
184        MzSpectrum { mz: mz_combined, intensity: intensity_combined }
185    }
186}