mscore/timstof/
spectrum.rs

1use std::collections::BTreeMap;
2use std::fmt;
3use std::fmt::{Display, Formatter};
4use crate::data::spectrum::{IndexedMzSpectrum, IndexedMzSpectrumVectorized, MsType, MzSpectrum};
5
6#[derive(Clone)]
7pub struct TimsSpectrumVectorized {
8    pub frame_id: i32,
9    pub scan: i32,
10    pub retention_time: f64,
11    pub mobility: f64,
12    pub ms_type: MsType,
13    pub vector: IndexedMzSpectrumVectorized,
14}
15
16#[derive(Clone, Debug)]
17pub struct TimsSpectrum {
18    pub frame_id: i32,
19    pub scan: i32,
20    pub retention_time: f64,
21    pub mobility: f64,
22    pub ms_type: MsType,
23    pub spectrum: IndexedMzSpectrum,
24}
25
26impl TimsSpectrum {
27    /// Creates a new `TimsSpectrum` instance.
28    ///
29    /// # Arguments
30    ///
31    /// * `frame_id` - index of frame in TDF raw file.
32    /// * `scan_id` - index of scan in TDF raw file.
33    /// * `retention_time` - The retention time in seconds.
34    /// * `mobility` - The inverse ion mobility.
35    /// * `spectrum` - A `TOFMzSpectrum` instance.
36    ///
37    /// # Examples
38    ///
39    /// ```
40    /// use mscore::data::spectrum::{IndexedMzSpectrum, MsType};
41    /// use mscore::timstof::spectrum::{TimsSpectrum};
42    ///
43    /// let spectrum = TimsSpectrum::new(1, 1, 100.0, 0.1, MsType::FragmentDda, IndexedMzSpectrum::new(vec![1000, 2000], vec![100.5, 200.5], vec![50.0, 60.0]));
44    /// ```
45    pub fn new(frame_id: i32, scan_id: i32, retention_time: f64, mobility: f64, ms_type: MsType, spectrum: IndexedMzSpectrum) -> Self {
46        TimsSpectrum { frame_id, scan: scan_id, retention_time, mobility: mobility, ms_type, spectrum }
47    }
48
49    pub fn filter_ranged(&self, mz_min: f64, mz_max: f64, intensity_min: f64, intensity_max: f64) -> Self {
50        let filtered = self.spectrum.filter_ranged(mz_min, mz_max, intensity_min, intensity_max);
51        TimsSpectrum { frame_id: self.frame_id, scan: self.scan, retention_time: self.retention_time, mobility: self.mobility, ms_type: self.ms_type.clone(), spectrum: filtered }
52    }
53
54    pub fn to_resolution(&self, resolution: i32) -> TimsSpectrum {
55        let spectrum = self.spectrum.to_resolution(resolution);
56        TimsSpectrum { frame_id: self.frame_id, scan: self.scan, retention_time: self.retention_time, mobility: self.mobility, ms_type: self.ms_type.clone(), spectrum }
57    }
58
59    pub fn vectorized(&self, resolution: i32) -> TimsSpectrumVectorized {
60        let vector = self.spectrum.vectorized(resolution);
61        TimsSpectrumVectorized { frame_id: self.frame_id, scan: self.scan, retention_time: self.retention_time, mobility: self.mobility, ms_type: self.ms_type.clone(), vector }
62    }
63
64    pub fn to_windows(&self, window_length: f64, overlapping: bool, min_peaks: usize, min_intensity: f64) -> BTreeMap<i32, TimsSpectrum> {
65
66        let mut splits: BTreeMap<i32, TimsSpectrum> = BTreeMap::new();
67
68        for (i, &mz) in self.spectrum.mz_spectrum.mz.iter().enumerate() {
69            let intensity = self.spectrum.mz_spectrum.intensity[i];
70            let tof = self.spectrum.index[i];
71
72            let tmp_key = (mz / window_length).floor() as i32;
73
74            splits.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
75                Vec::new(), Vec::new(), Vec::new()))
76            ).spectrum.mz_spectrum.mz.push(mz);
77
78            splits.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
79                Vec::new(), Vec::new(), Vec::new()))
80            ).spectrum.mz_spectrum.intensity.push(intensity);
81
82            splits.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
83                Vec::new(), Vec::new(), Vec::new()))
84            ).spectrum.index.push(tof);
85        }
86
87        if overlapping {
88            let mut splits_offset = BTreeMap::new();
89
90            for (i, &mmz) in self.spectrum.mz_spectrum.mz.iter().enumerate() {
91                let intensity = self.spectrum.mz_spectrum.intensity[i];
92                let tof = self.spectrum.index[i];
93
94                let tmp_key = -((mmz + window_length / 2.0) / window_length).floor() as i32;
95
96                splits_offset.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
97                    Vec::new(), Vec::new(), Vec::new()))
98                ).spectrum.mz_spectrum.mz.push(mmz);
99
100                splits_offset.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
101                    Vec::new(), Vec::new(), Vec::new()))
102                ).spectrum.mz_spectrum.intensity.push(intensity);
103
104                splits_offset.entry(tmp_key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
105                    Vec::new(), Vec::new(), Vec::new()))
106                ).spectrum.index.push(tof);
107            }
108
109            for (key, val) in splits_offset {
110                splits.entry(key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
111                    Vec::new(), Vec::new(), Vec::new()))
112                ).spectrum.mz_spectrum.mz.extend(val.spectrum.mz_spectrum.mz);
113
114                splits.entry(key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
115                    Vec::new(), Vec::new(), Vec::new()))
116                ).spectrum.mz_spectrum.intensity.extend(val.spectrum.mz_spectrum.intensity);
117
118                splits.entry(key).or_insert_with(|| TimsSpectrum::new(self.frame_id, self.scan, self.retention_time, self.mobility, self.ms_type.clone(), IndexedMzSpectrum::new(
119                    Vec::new(), Vec::new(), Vec::new()))
120                ).spectrum.index.extend(val.spectrum.index);
121            }
122        }
123
124        splits.retain(|_, spectrum| {
125            spectrum.spectrum.mz_spectrum.mz.len() >= min_peaks && spectrum.spectrum.mz_spectrum.intensity.iter().cloned().max_by(
126                |a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)).unwrap_or(0.0) >= min_intensity
127        });
128
129        splits
130    }
131}
132
133// implement default (empty TimsSpectrum) constructor
134impl Default for TimsSpectrum {
135    fn default() -> Self {
136        TimsSpectrum { frame_id: 0, scan: 0, retention_time: 0.0, mobility: 0.0, ms_type: MsType::Unknown, spectrum: IndexedMzSpectrum::default() }
137    }
138}
139
140impl std::ops::Add for TimsSpectrum {
141    type Output = Self;
142
143    fn add(self, other: Self) -> TimsSpectrum {
144        assert_eq!(self.frame_id, other.frame_id);
145
146        let average_mobility = (self.mobility + other.mobility) / 2.0;
147        let average_retention_time = (self.retention_time + other.retention_time) / 2.0;
148        let average_scan_floor = ((self.scan as f64 + other.scan as f64) / 2.0) as i32;
149
150        let mut combined_map: BTreeMap<i64, (f64, i32, i32)> = BTreeMap::new();
151        let quantize = |mz: f64| -> i64 { (mz * 1_000_000.0).round() as i64 };
152
153        for ((mz, intensity), index) in self.spectrum.mz_spectrum.mz.iter().zip(self.spectrum.mz_spectrum.intensity.iter()).zip(self.spectrum.index.iter()) {
154            let key = quantize(*mz);
155            combined_map.insert(key, (*intensity, *index, 1)); // Initialize count as 1
156        }
157
158        for ((mz, intensity), index) in other.spectrum.mz_spectrum.mz.iter().zip(other.spectrum.mz_spectrum.intensity.iter()).zip(other.spectrum.index.iter()) {
159            let key = quantize(*mz);
160            combined_map.entry(key).and_modify(|e| {
161                e.0 += *intensity; // Sum intensity
162                e.1 += *index;     // Sum index
163                e.2 += 1;          // Increment count
164            }).or_insert((*intensity, *index, 1));
165        }
166
167        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
168        let intensity_combined: Vec<f64> = combined_map.values().map(|(intensity, _, _)| *intensity).collect();
169        let index_combined: Vec<i32> = combined_map.values().map(|(_, index, count)| index / count).collect(); // Average index
170
171        let spectrum = IndexedMzSpectrum { index: index_combined, mz_spectrum: MzSpectrum { mz: mz_combined, intensity: intensity_combined } };
172        TimsSpectrum { frame_id: self.frame_id, scan: average_scan_floor, retention_time: average_retention_time, mobility: average_mobility, ms_type: self.ms_type.clone(), spectrum }
173    }
174}
175
176impl Display for TimsSpectrum {
177    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
178        write!(f, "TimsSpectrum(frame_id: {}, scan_id: {}, retention_time: {}, mobility: {}, spectrum: {})", self.frame_id, self.scan, self.retention_time, self.mobility, self.spectrum)
179    }
180}
mscore/timstof/spectrum.rs

mscore/timstof/
spectrum.rs
