mzdeisotope_map/
feature_fit.rs

1/*! Description of isotopic pattern fits */
2use std::cmp::Ordering;
3
4use chemical_elements::isotopic_pattern::TheoreticalIsotopicPattern;
5
6use mzsignal::smooth::moving_average;
7use mzdeisotope::scorer::ScoreType;
8use mzpeaks::{feature::TimeInterval, feature_map::FeatureMapLike, CoordinateRange};
9
10use crate::traits::FeatureMapType;
11
12#[derive(Debug, Default, Clone, Copy, PartialEq, PartialOrd)]
13pub struct MapCoordinate {
14    pub coord: f64,
15    pub time: f64,
16}
17
18impl MapCoordinate {
19    pub fn new(coord: f64, time: f64) -> Self {
20        Self { coord, time }
21    }
22}
23
24#[derive(Debug, Clone)]
25pub struct FeatureSetFit {
26    pub features: Vec<Option<usize>>,
27    pub start: MapCoordinate,
28    pub end: MapCoordinate,
29    pub score: ScoreType,
30    pub theoretical: TheoreticalIsotopicPattern,
31    pub mz: f64,
32    pub neutral_mass: f64,
33    pub charge: i32,
34    pub missing_features: usize,
35    pub n_points: usize,
36    pub monoisotopic_index: usize,
37    pub scores: Vec<ScoreType>,
38    pub times: Vec<f64>,
39}
40
41impl PartialOrd for FeatureSetFit {
42    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
43        Some(self.cmp(other))
44    }
45}
46
47impl Ord for FeatureSetFit {
48    fn cmp(&self, other: &Self) -> Ordering {
49        self.score.total_cmp(&other.score)
50    }
51}
52
53impl PartialEq for FeatureSetFit {
54    fn eq(&self, other: &Self) -> bool {
55        self.features == other.features
56            && self.score == other.score
57            && self.theoretical == other.theoretical
58            && self.charge == other.charge
59    }
60}
61
62impl Eq for FeatureSetFit {}
63
64impl FeatureSetFit {
65    #[allow(clippy::too_many_arguments)]
66    pub fn new(
67        features: Vec<Option<usize>>,
68        theoretical: TheoreticalIsotopicPattern,
69        start: MapCoordinate,
70        end: MapCoordinate,
71        score: ScoreType,
72        charge: i32,
73        missing_features: usize,
74        neutral_mass: f64,
75        n_points: usize,
76        scores: Vec<ScoreType>,
77        times: Vec<f64>,
78    ) -> Self {
79        let mz = theoretical.peaks.first().map(|p| p.mz).unwrap();
80        Self {
81            features,
82            theoretical,
83            start,
84            end,
85            score,
86            charge,
87            missing_features,
88            monoisotopic_index: 0,
89            mz,
90            neutral_mass,
91            n_points,
92            scores,
93            times,
94        }
95    }
96
97    pub fn count_null_features(&self) -> usize {
98        self.features.iter().map(|f| f.is_none() as usize).sum()
99    }
100
101    pub fn has_multiple_real_features(&self) -> bool {
102        (self.features.len() - self.count_null_features()) > 1
103    }
104
105    pub fn len(&self) -> usize {
106        self.features.len()
107    }
108
109    pub fn is_empty(&self) -> bool {
110        self.features.is_empty()
111    }
112
113    pub fn find_bounds<Y: Clone>(&self, feature_map: &FeatureMapType<Y>, detection_threshold: f32) -> CoordinateRange<Y> {
114        let mut start_time: f64 = f64::INFINITY;
115        let mut end_time: f64 = f64::INFINITY;
116
117        for (f, p) in self.features.iter().zip(self.theoretical.iter()) {
118            if let Some(i) = f {
119                let f = feature_map.get_item(*i);
120                if p.intensity() >= detection_threshold {
121                    start_time = f.start_time().unwrap().min(start_time);
122                    end_time = f.end_time().unwrap().min(end_time);
123                }
124
125            } else {
126                continue;
127            }
128        }
129
130        CoordinateRange::new(Some(start_time), Some(end_time))
131    }
132
133    pub fn find_separation<Y: Clone>(&self, feature_map: &FeatureMapType<Y>, detection_threshold: f32) -> (CoordinateRange<Y>, Vec<ScoreSegment>) {
134        let mut time_range = self.find_bounds(feature_map, detection_threshold);
135        if self.n_points > 0 {
136            let mut segments = Vec::new();
137            let mut last_score = ScoreType::INFINITY;
138            let mut scores = Vec::with_capacity(self.scores.len());
139            scores.resize(self.scores.len(), 0.0);
140            moving_average::<f32, 1>(&self.scores, &mut scores);
141            let mut begin_i = 0;
142            for (i, score) in scores.iter().copied().enumerate() {
143                if score > 0.0 && last_score < 0.0 {
144                    begin_i = i
145                } else if score < 0.0 && last_score > 0.0 {
146                    let end_i = i;
147                    segments.push(ScoreSegment::new(
148                        begin_i, end_i,
149                        scores[begin_i..end_i].iter().copied().sum::<f32>()
150                    ));
151                    begin_i = i;
152                }
153                last_score = score;
154            }
155            let end_i = scores.len().saturating_sub(1);
156            segments.push(ScoreSegment::new(
157                begin_i,
158                end_i,
159                scores[begin_i..end_i].iter().sum()
160            ));
161            let segment = segments.iter().max_by(|a, b| {a.score.total_cmp(&b.score)}).copied().unwrap();
162            if time_range.start.unwrap() < self.times[segment.start] && segment.start != 0 {
163                time_range.start = Some(self.times[segment.start]);
164            }
165            if time_range.end.unwrap() > self.times[segment.end] {
166                time_range.end = Some(self.times[segment.end])
167            }
168            (time_range, segments)
169        } else {
170            (time_range, Vec::new())
171        }
172    }
173}
174
175
176#[derive(Debug, Default, Clone, Copy, PartialEq)]
177pub struct ScoreSegment {
178    pub start: usize,
179    pub end: usize,
180    pub score: f32
181}
182
183impl ScoreSegment {
184    pub fn new(start: usize, end: usize, score: f32) -> Self {
185        Self { start, end, score }
186    }
187}
mzdeisotope_map/feature_fit.rs

mzdeisotope_map/
feature_fit.rs
