Loading...
Searching...
No Matches
highrespeakpicker.cpp
Go to the documentation of this file.
// Copyright 2026 Filippo Rusconi
// Inspired by work by Lars Nilse in OpenMS
/////////////////////// stdlib includes
#include <limits>
/////////////////////// Qt includes
#include <QList>
#include <QDebug>
/////////////////////// pappsomspp includes
/////////////////////// Local includes
#include "pappsomspp/core/processing/detection/highrespeakpicker.h"
#include "pappsomspp/core/processing/detection/cubicsplinemodel.h"
#include "pappsomspp/core/processing/detection/localsignaltonoiseestimator.h"
namespace pappso
{
HighResPeakPicker::HighResPeakPicker(Parameters &parameters) : m_parameters(parameters)
{
}
HighResPeakPicker::~HighResPeakPicker()
{
}
void
HighResPeakPicker::pick(const Trace &trace, Trace &picked_peaks, QList<PeakRange> &peak_ranges)
{
picked_peaks.clear();
qsizetype trace_size = trace.size();
if(trace_size < 5)
return;
// signal-to-noise estimation
LocalSignalToNoiseEstimator::Parameters signal_to_noise_estimator_parameters;
LocalSignalToNoiseEstimator local_signal_estimator(signal_to_noise_estimator_parameters);
if(m_parameters.signalToNoise > 0.0)
{
// Computes, for each DataPoint in the trace, an estimage of local
// signal-to-ratio value that is consulted below.
local_signal_estimator.initialize(trace);
}
// Iterate in the trace and search for local maxima, that is, roughly,
// DataPoint instances that are higher than previous and next point.
for(qsizetype iter = 2; iter < trace_size - 2; ++iter)
{
// Skip any data point that is of an intensity below min_intensity
if(trace.at(iter).y < m_parameters.minimumIntensity)
{
// qDebug() << "Now skipping low intensity peak at m/z:" << trace.at(iter).x
// << "with intensity:" << trace.at(iter).y;
continue;
}
double center_data_point_mz = trace.at(iter).x;
double center_data_point_int = trace.at(iter).y;
double left_data_point_mz = trace.at(iter - 1).x;
double left_data_point_int = trace.at(iter - 1).y;
double right_data_point_mz = trace.at(iter + 1).x;
double right_data_point_int = trace.at(iter + 1).y;
// Do not interpolate when the left or right DataPoint is zero-intensity.
if(std::fabs(left_data_point_int) <
std::numeric_limits<double>::epsilon())
{
continue;
}
if(std::fabs(right_data_point_int) <
std::numeric_limits<double>::epsilon())
{
continue;
}
double current_snt_estimate = 0.0;
double current_snt_estimate_left_1 = 0.0;
double current_snt_estimate_right_1 = 0.0;
if(m_parameters.signalToNoise > 0.0)
{
current_snt_estimate =
local_signal_estimator.getSignalToNoiseRatio(iter);
// qDebug() << "For current center point at iter:" << iter
// << "with mz:" << center_data_point_mz << "of intensity"
// << center_data_point_int
// << "the snt estimate is:" << current_snt_estimate;
current_snt_estimate_left_1 =
local_signal_estimator.getSignalToNoiseRatio(iter - 1);
current_snt_estimate_right_1 =
local_signal_estimator.getSignalToNoiseRatio(iter + 1);
}
// Look for "peaking" data points meeting MZ and intensity/SNT criteria
// The currently iterated data point should at least have its
// intensity greater than both the left and right data points intensities.
// If the signal-to-noise checks are OK, than the data point becomes a peak "core",
// that is the top point of the peak shape.
if((center_data_point_int > left_data_point_int) &&
(center_data_point_int > right_data_point_int) &&
(current_snt_estimate >= m_parameters.signalToNoise) &&
(current_snt_estimate_left_1 >= m_parameters.signalToNoise) &&
(current_snt_estimate_right_1 >= m_parameters.signalToNoise))
{
// First check if the currently iterated center center point
// is surrounded by more intense peaks on its left-left or on its
// right-right, because that might indicate signal instability instead
// of an actual maximum. In this case, we do not consider the data
// point.
double current_snt_estimate_left_2 = 0.0;
double current_snt_estimate_right_2 = 0.0;
if(m_parameters.signalToNoise > 0.0)
{
current_snt_estimate_left_2 =
local_signal_estimator.getSignalToNoiseRatio(iter - 2);
current_snt_estimate_right_2 =
local_signal_estimator.getSignalToNoiseRatio(iter + 2);
}
// Checking signal-to-noise
if((iter > 1) && (iter + 2 < static_cast<qsizetype>(trace.size())) &&
(left_data_point_int < trace.at(iter - 2).y) &&
(right_data_point_int < trace.at(iter + 2).y) &&
(current_snt_estimate_left_2 >= m_parameters.signalToNoise) &&
(current_snt_estimate_right_2 >= m_parameters.signalToNoise))
{
// If the point left of the left data point (with respect to
// current center data point that was found to be a local maximum)
// and if the point right of the right data point (again, with
// respect to current center data point that was found to be a
// local maximum) and both of these left_left and right_right
// datapoints are above the signal to noise ratio, then we are
// experiencing signal instability.
++iter;
continue;
}
QMap<double, double> peak_raw_data;
peak_raw_data[center_data_point_mz] = center_data_point_int;
peak_raw_data[left_data_point_mz] = left_data_point_int;
peak_raw_data[right_data_point_mz] = right_data_point_int;
// Peak core found, now extend it to the left
qsizetype k = 2;
// No need to extend peak if previous intensity was zero
bool previous_zero_left = false;
qsizetype missing_left_count = 0;
// Index of the left boundary for the spline interpolation
qsizetype left_boundary = iter - 1;
while((k <= iter) && // Prevent underflow
(iter - k + 1 > 0) && !previous_zero_left &&
(trace.at(iter - k).y <= peak_raw_data.begin().value()))
{
double current_snt_estimate_left_k = 0.0;
if(m_parameters.signalToNoise > 0.0)
{
current_snt_estimate_left_k =
local_signal_estimator.getSignalToNoiseRatio(iter - k);
}
if(current_snt_estimate_left_k >= m_parameters.signalToNoise)
{
peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
}
{
++missing_left_count;
if(missing_left_count <= m_parameters.missingPeakCount)
{
peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
}
}
previous_zero_left = (trace.at(iter - k).y == 0);
left_boundary = iter - k;
++k;
}
// Now do the same on the right
k = 2;
// No need to extend peak if previous intensity was zero
bool previous_zero_right = false;
qsizetype missing_right_count = 0;
// Index of the right boundary for the spline interpolation
qsizetype right_boundary(iter + 1);
while((iter + k < static_cast<qsizetype>(trace.size())) &&
!previous_zero_right &&
(trace.at(iter + k).y <= peak_raw_data.last()))
{
double current_snt_estimate_right_k = 0.0;
if(m_parameters.signalToNoise > 0.0)
{
current_snt_estimate_right_k =
local_signal_estimator.getSignalToNoiseRatio(iter + k);
}
if(current_snt_estimate_right_k >= m_parameters.signalToNoise)
{
peak_raw_data[trace.at(iter + k).x] = trace.at(iter + k).y;
}
else
{
++missing_right_count;
if(missing_right_count <= m_parameters.missingPeakCount)
{
peak_raw_data[trace.at(iter + k).x] =
trace.at(iter + k).y;
}
}
previous_zero_right = (trace.at(iter + k).y == 0);
right_boundary = iter + k;
++k;
}
// Skip peak altogether if the minimal number of 3 points for fitting is not reached
if(peak_raw_data.size() < 3)
{
continue;
}
// qDebug() << "peak_raw_data size:" << peak_raw_data.size();
CubicSplineModel cubic_spline_model(peak_raw_data);
// qDebug() << "Stack-allocated:" << &cubic_spline_model;
// qDebug() << "The cubic spline model has "
// << cubic_spline_model.getKnots().size() << "knots.";
// calculate maximum by evaluating the spline's 1st derivative
// (bisection method)
double max_peak_mz = center_data_point_mz;
double max_peak_int = center_data_point_int;
double threshold = 1e-6;
spline_bisection(cubic_spline_model,
left_data_point_mz,
right_data_point_mz,
max_peak_mz,
max_peak_int,
threshold);
// save picked peak into output spectrum
DataPoint data_point;
PeakRange peak_range;
data_point.x = max_peak_mz;
data_point.y = max_peak_int;
peak_range.mz_start = trace.at(left_boundary).x;
peak_range.mz_stop = trace.at(right_boundary).x;
// qDebug() << "Now appending new data point:" <<
// data_point.toString();
picked_peaks.append(data_point);
peak_ranges.push_back(peak_range);
// jump over profile data points that have been considered already
iter += k - 1;
}
}
// qDebug() << "Now returning" << picked_peaks.size()
// << "picked centroids:" << picked_peaks.toString();
}
} // namespace pappso
pappso::CubicSplineModel
Definition cubicsplinemodel.h:24
pappso::HighResPeakPicker::pick
void pick(const Trace &trace, Trace &picked_peaks, QList< PeakRange > &peak_ranges)
Definition highrespeakpicker.cpp:33
pappso::HighResPeakPicker::~HighResPeakPicker
virtual ~HighResPeakPicker()
Definition highrespeakpicker.cpp:28
pappso::HighResPeakPicker::m_parameters
Parameters m_parameters
Definition highrespeakpicker.h:47
pappso::HighResPeakPicker::HighResPeakPicker
HighResPeakPicker(Parameters &parameters)
Definition highrespeakpicker.cpp:24
pappso::LocalSignalToNoiseEstimator
Definition localsignaltonoiseestimator.h:26
pappso::LocalSignalToNoiseEstimator::initialize
void initialize(const Trace &trace)
Definition localsignaltonoiseestimator.cpp:65
pappso::LocalSignalToNoiseEstimator::getSignalToNoiseRatio
double getSignalToNoiseRatio(qsizetype index) const
Definition localsignaltonoiseestimator.cpp:351
pappso::Trace
A simple container of DataPoint instances.
Definition trace.h:152
pappso::Trace::append
size_t append(const DataPoint &data_point)
appends a datapoint and return new size
Definition trace.cpp:657
cubicsplinemodel.h
highrespeakpicker.h
localsignaltonoiseestimator.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
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
pappso::PeptideIonCter::y
@ y
Definition peptide.h:87
pappso::DataPoint
Definition datapoint.h:22
pappso::DataPoint::x
pappso_double x
Definition datapoint.h:24
pappso::DataPoint::y
pappso_double y
Definition datapoint.h:25
pappso::HighResPeakPicker::Parameters
Definition highrespeakpicker.h:29
pappso::HighResPeakPicker::PeakRange
Definition highrespeakpicker.h:39
pappso::HighResPeakPicker::PeakRange::mz_start
double mz_start
Definition highrespeakpicker.h:40
pappso::HighResPeakPicker::PeakRange::mz_stop
double mz_stop
Definition highrespeakpicker.h:41
pappso::LocalSignalToNoiseEstimator::Parameters
Definition localsignaltonoiseestimator.h:38