libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
pappso::HighResPeakPicker Class Reference

#include <highrespeakpicker.h>

Classes

struct  Parameters
struct  PeakRange

Public Member Functions

 HighResPeakPicker (Parameters &parameters)
virtual ~HighResPeakPicker ()
void pick (const Trace &trace, Trace &picked_peaks, QList< PeakRange > &peak_ranges)

Protected Attributes

Parameters m_parameters

Detailed Description

Definition at line 24 of file highrespeakpicker.h.

Constructor & Destructor Documentation

◆ HighResPeakPicker()

pappso::HighResPeakPicker::HighResPeakPicker ( Parameters & parameters)

Definition at line 24 of file highrespeakpicker.cpp.

24 : m_parameters(parameters)
25{
26}

References m_parameters.

◆ ~HighResPeakPicker()

pappso::HighResPeakPicker::~HighResPeakPicker ( )
virtual

Definition at line 28 of file highrespeakpicker.cpp.

29{
30}

Member Function Documentation

◆ pick()

void pappso::HighResPeakPicker::pick ( const Trace & trace,
Trace & picked_peaks,
QList< PeakRange > & peak_ranges )

Definition at line 33 of file highrespeakpicker.cpp.

36{
37 picked_peaks.clear();
38
39 qsizetype trace_size = trace.size();
40
41 if(trace_size < 5)
42 return;
43
44 // signal-to-noise estimation
45 LocalSignalToNoiseEstimator::Parameters signal_to_noise_estimator_parameters;
46 LocalSignalToNoiseEstimator local_signal_estimator(signal_to_noise_estimator_parameters);
47
48 if(m_parameters.signalToNoise > 0.0)
49 {
50 local_signal_estimator.initialize(trace);
51 }
52
53 // Iterate in the trace and search for local maxima, that is, roughly,
54 // DataPoint instances that are higher than previous and next point.
55 for(qsizetype iter = 2; iter < trace_size - 2; ++iter)
56 {
57 double center_data_point_mz = trace.at(iter).x;
58 double center_data_point_int = trace.at(iter).y;
59 double left_data_point_mz = trace.at(iter - 1).x;
60 double left_data_point_int = trace.at(iter - 1).y;
61 double right_data_point_mz = trace.at(iter + 1).x;
62 double right_data_point_int = trace.at(iter + 1).y;
63
64 // Do not interpolate when the left or right DataPoint is zero-intensity.
65 if(std::fabs(left_data_point_int) <
66 std::numeric_limits<double>::epsilon())
67 {
68 continue;
69 }
70 if(std::fabs(right_data_point_int) <
71 std::numeric_limits<double>::epsilon())
72 {
73 continue;
74 }
75
76 double current_snt_estimate = 0.0;
77 double current_snt_estimate_left_1 = 0.0;
78 double current_snt_estimate_right_1 = 0.0;
79
80 if(m_parameters.signalToNoise > 0.0)
81 {
82 current_snt_estimate =
83 local_signal_estimator.getSignalToNoiseRatio(iter);
84
85 // qDebug() << "For current center point at iter:" << iter
86 // << "with mz:" << center_data_point_mz << "of intensity"
87 // << center_data_point_int
88 // << "the snt estimate is:" << current_snt_estimate;
89
90 current_snt_estimate_left_1 =
91 local_signal_estimator.getSignalToNoiseRatio(iter - 1);
92 current_snt_estimate_right_1 =
93 local_signal_estimator.getSignalToNoiseRatio(iter + 1);
94 }
95
96 // look for peak cores meeting MZ and intensity/SNT criteria
97 if((center_data_point_int > left_data_point_int) &&
98 (center_data_point_int > right_data_point_int) &&
99 (current_snt_estimate >= m_parameters.signalToNoise) &&
100 (current_snt_estimate_left_1 >= m_parameters.signalToNoise) &&
101 (current_snt_estimate_right_1 >= m_parameters.signalToNoise))
102 {
103 // First check if the currently iterated center center point
104 // is surrounded by more intense peaks on its left-left or on its
105 // right-right, because that might indicate signal instability instead
106 // of an actual maximum. In this case, we do not consider the data
107 // point.
108
109 double current_snt_estimate_left_2 = 0.0;
110 double current_snt_estimate_right_2 = 0.0;
111
112 if(m_parameters.signalToNoise > 0.0)
113 {
114 current_snt_estimate_left_2 =
115 local_signal_estimator.getSignalToNoiseRatio(iter - 2);
116 current_snt_estimate_right_2 =
117 local_signal_estimator.getSignalToNoiseRatio(iter + 2);
118 }
119
120 // Checking signal-to-noise
121 if((iter > 1) && (iter + 2 < static_cast<qsizetype>(trace.size())) &&
122 (left_data_point_int < trace.at(iter - 2).y) &&
123 (right_data_point_int < trace.at(iter + 2).y) &&
124 (current_snt_estimate_left_2 >= m_parameters.signalToNoise) &&
125 (current_snt_estimate_right_2 >= m_parameters.signalToNoise))
126 {
127 // If the point left of the left data point (with respect to
128 // current center data point that was found to be a local maximum)
129 // and if the point right of the right data point (again, with
130 // respect to current center data point that was found to be a
131 // local maximum) and both of these left_left and right_right
132 // datapoints are above the signal to noise ratio, then we are
133 // experiencing signal instability.
134 ++iter;
135 continue;
136 }
137
138 QMap<double, double> peak_raw_data;
139
140 peak_raw_data[center_data_point_mz] = center_data_point_int;
141 peak_raw_data[left_data_point_mz] = left_data_point_int;
142 peak_raw_data[right_data_point_mz] = right_data_point_int;
143
144 // peak core found, now extend it
145 // to the left
146 qsizetype k = 2;
147
148 // No need to extend peak if previous intensity was zero
149 bool previous_zero_left = false;
150 qsizetype missing_left_count = 0;
151 // Index of the left boundary for the spline interpolation
152 qsizetype left_boundary = iter - 1;
153
154 while((k <= iter) && // prevent underflow
155 (iter - k + 1 > 0) && !previous_zero_left &&
156 (trace.at(iter - k).y <= peak_raw_data.begin().value()))
157 {
158 double current_snt_estimate_left_k = 0.0;
159
160 if(m_parameters.signalToNoise > 0.0)
161 {
162 current_snt_estimate_left_k =
163 local_signal_estimator.getSignalToNoiseRatio(iter - k);
164 }
165
166 if(current_snt_estimate_left_k >= m_parameters.signalToNoise)
167 {
168 peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
169 }
170 {
171 ++missing_left_count;
172 if(missing_left_count <= m_parameters.missingPeakCount)
173 {
174 peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
175 }
176 }
177
178 previous_zero_left = (trace.at(iter - k).y == 0);
179 left_boundary = iter - k;
180 ++k;
181 }
182
183 // to the right
184 k = 2;
185
186 // No need to extend peak if previous intensity was zero
187 bool previous_zero_right = false;
188 qsizetype missing_right_count = 0;
189 // Index of the right boundary for the spline interpolation
190 qsizetype right_boundary(iter + 1);
191
192 while((iter + k < static_cast<qsizetype>(trace.size())) &&
193 !previous_zero_right &&
194 (trace.at(iter + k).y <= peak_raw_data.last()))
195 {
196 double current_snt_estimate_right_k = 0.0;
197
198 if(m_parameters.signalToNoise > 0.0)
199 {
200 current_snt_estimate_right_k =
201 local_signal_estimator.getSignalToNoiseRatio(iter + k);
202 }
203
204 if(current_snt_estimate_right_k >= m_parameters.signalToNoise)
205 {
206 peak_raw_data[trace.at(iter + k).x] = trace.at(iter + k).y;
207 }
208 else
209 {
210 ++missing_right_count;
211 if(missing_right_count <= m_parameters.missingPeakCount)
212 {
213 peak_raw_data[trace.at(iter + k).x] =
214 trace.at(iter + k).y;
215 }
216 }
217
218 previous_zero_right = (trace.at(iter + k).y == 0);
219 right_boundary = iter + k;
220 ++k;
221 }
222
223 // skip if the minimal number of 3 points for fitting is not reached
224 if(peak_raw_data.size() < 3)
225 {
226 continue;
227 }
228
229 // qDebug() << "peak_raw_data size:" << peak_raw_data.size();
230
231 CubicSplineModel cubic_spline_model(peak_raw_data);
232 // qDebug() << "Stack-allocated:" << &cubic_spline_model;
233
234
235 // qDebug() << "The cubic spline model has "
236 // << cubic_spline_model.getKnots().size() << "knots.";
237
238 // calculate maximum by evaluating the spline's 1st derivative
239 // (bisection method)
240 double max_peak_mz = center_data_point_mz;
241 double max_peak_int = center_data_point_int;
242 double threshold = 1e-6;
243 spline_bisection(cubic_spline_model,
244 left_data_point_mz,
245 right_data_point_mz,
246 max_peak_mz,
247 max_peak_int,
248 threshold);
249
250 // save picked peak into output spectrum
251 DataPoint data_point;
252
253 PeakRange peak_range;
254
255 data_point.x = max_peak_mz;
256 data_point.y = max_peak_int;
257 peak_range.mz_start = trace.at(left_boundary).x;
258 peak_range.mz_stop = trace.at(right_boundary).x;
259
260 // qDebug() << "Now appending new data point:" <<
261 // data_point.toString();
262
263 picked_peaks.append(data_point);
264
265 peak_ranges.push_back(peak_range);
266
267 // jump over profile data points that have been considered already
268 iter += k - 1;
269 }
270 }
271
272 // qDebug() << "Now returning" << picked_peaks.size()
273 // << "picked centroids:" << picked_peaks.toString();
274}
pappso::spline_bisection
void spline_bisection(const CubicSplineModel &spline_model, double const mz_at_left, double const mz_at_right, double &center_peak_mz, double &center_peak_intensity, double const threshold)
Definition cubicsplinemodel.cpp:241

References pappso::Trace::append(), pappso::LocalSignalToNoiseEstimator::getSignalToNoiseRatio(), pappso::LocalSignalToNoiseEstimator::initialize(), m_parameters, pappso::HighResPeakPicker::PeakRange::mz_start, pappso::HighResPeakPicker::PeakRange::mz_stop, pappso::spline_bisection(), pappso::DataPoint::x, pappso::DataPoint::y, and pappso::y.

Member Data Documentation

◆ m_parameters

Parameters pappso::HighResPeakPicker::m_parameters
protected

Definition at line 46 of file highrespeakpicker.h.

Referenced by HighResPeakPicker(), and pick().

The documentation for this class was generated from the following files: