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        // Build intermediate structures with raw Vecs first
66        let mut splits_raw: BTreeMap<i32, (Vec<f64>, Vec<f64>, Vec<i32>)> = 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            let entry = splits_raw.entry(tmp_key).or_insert_with(|| (Vec::new(), Vec::new(), Vec::new()));
75            entry.0.push(mz);
76            entry.1.push(intensity);
77            entry.2.push(tof);
78        }
79
80        if overlapping {
81            let mut splits_offset: BTreeMap<i32, (Vec<f64>, Vec<f64>, Vec<i32>)> = BTreeMap::new();
82
83            for (i, &mmz) in self.spectrum.mz_spectrum.mz.iter().enumerate() {
84                let intensity = self.spectrum.mz_spectrum.intensity[i];
85                let tof = self.spectrum.index[i];
86
87                let tmp_key = -((mmz + window_length / 2.0) / window_length).floor() as i32;
88
89                let entry = splits_offset.entry(tmp_key).or_insert_with(|| (Vec::new(), Vec::new(), Vec::new()));
90                entry.0.push(mmz);
91                entry.1.push(intensity);
92                entry.2.push(tof);
93            }
94
95            for (key, val) in splits_offset {
96                let entry = splits_raw.entry(key).or_insert_with(|| (Vec::new(), Vec::new(), Vec::new()));
97                entry.0.extend(val.0);
98                entry.1.extend(val.1);
99                entry.2.extend(val.2);
100            }
101        }
102
103        // Convert to TimsSpectrum and filter
104        let mut splits: BTreeMap<i32, TimsSpectrum> = BTreeMap::new();
105        for (key, (mz_vec, intensity_vec, index_vec)) in splits_raw {
106            if mz_vec.len() >= min_peaks {
107                let max_intensity = intensity_vec.iter().cloned().max_by(
108                    |a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)
109                ).unwrap_or(0.0);
110                if max_intensity >= min_intensity {
111                    let spectrum = IndexedMzSpectrum::new(index_vec, mz_vec, intensity_vec);
112                    splits.insert(key, TimsSpectrum::new(
113                        self.frame_id, self.scan, self.retention_time, self.mobility,
114                        self.ms_type.clone(), spectrum
115                    ));
116                }
117            }
118        }
119
120        splits
121    }
122}
123
124// implement default (empty TimsSpectrum) constructor
125impl Default for TimsSpectrum {
126    fn default() -> Self {
127        TimsSpectrum { frame_id: 0, scan: 0, retention_time: 0.0, mobility: 0.0, ms_type: MsType::Unknown, spectrum: IndexedMzSpectrum::default() }
128    }
129}
130
131impl std::ops::Add for TimsSpectrum {
132    type Output = Self;
133
134    fn add(self, other: Self) -> TimsSpectrum {
135        assert_eq!(self.frame_id, other.frame_id);
136
137        let average_mobility = (self.mobility + other.mobility) / 2.0;
138        let average_retention_time = (self.retention_time + other.retention_time) / 2.0;
139        let average_scan_floor = ((self.scan as f64 + other.scan as f64) / 2.0) as i32;
140
141        let mut combined_map: BTreeMap<i64, (f64, i32, i32)> = BTreeMap::new();
142        let quantize = |mz: f64| -> i64 { (mz * 1_000_000.0).round() as i64 };
143
144        for ((mz, intensity), index) in self.spectrum.mz_spectrum.mz.iter().zip(self.spectrum.mz_spectrum.intensity.iter()).zip(self.spectrum.index.iter()) {
145            let key = quantize(*mz);
146            combined_map.insert(key, (*intensity, *index, 1)); // Initialize count as 1
147        }
148
149        for ((mz, intensity), index) in other.spectrum.mz_spectrum.mz.iter().zip(other.spectrum.mz_spectrum.intensity.iter()).zip(other.spectrum.index.iter()) {
150            let key = quantize(*mz);
151            combined_map.entry(key).and_modify(|e| {
152                e.0 += *intensity; // Sum intensity
153                e.1 += *index;     // Sum index
154                e.2 += 1;          // Increment count
155            }).or_insert((*intensity, *index, 1));
156        }
157
158        let mz_combined: Vec<f64> = combined_map.keys().map(|&key| key as f64 / 1_000_000.0).collect();
159        let intensity_combined: Vec<f64> = combined_map.values().map(|(intensity, _, _)| *intensity).collect();
160        let index_combined: Vec<i32> = combined_map.values().map(|(_, index, count)| index / count).collect(); // Average index
161
162        let spectrum = IndexedMzSpectrum { index: index_combined, mz_spectrum: MzSpectrum::new(mz_combined, intensity_combined) };
163        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 }
164    }
165}
166
167impl Display for TimsSpectrum {
168    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
169        write!(f, "TimsSpectrum(frame_id: {}, scan_id: {}, retention_time: {}, mobility: {}, spectrum: {})", self.frame_id, self.scan, self.retention_time, self.mobility, self.spectrum)
170    }
171}
mscore/timstof/spectrum.rs

mscore/timstof/
spectrum.rs
