Loading...
Searching...
No Matches
filterlowintensitysignalremoval.cpp
Go to the documentation of this file.
/* BEGIN software license
*
* msXpertSuite - mass spectrometry software suite
* -----------------------------------------------
* Copyright(C) 2009,...,2021 Filippo Rusconi
*
* http://www.msxpertsuite.org
*
* This file is part of the msXpertSuite project.
*
* The msXpertSuite project is the successor of the massXpert project. This
* project now includes various independent modules:
*
* - massXpert, model polymer chemistries and simulate mass spectrometric data;
* - mineXpert, a powerful TIC chromatogram/mass spectrum viewer/miner;
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* END software license
*/
#include <qmath.h>
#include <limits>
#include <iterator>
#include <QVector>
#include <QDebug>
#include "pappsomspp/core/exception/exceptionnotrecognized.h"
#include "filterlowintensitysignalremoval.h"
#include "../../trace/maptrace.h"
namespace pappso
{
FilterLowIntensitySignalRemoval::FilterLowIntensitySignalRemoval(double mean,
double std_dev,
double threshold)
{
// qDebug() << "Building with parameters: [ mean:" << mean << "- std_dev:" << std_dev
// << "- threshold" << threshold << "]";
m_noiseMean = mean;
m_noiseStdDev = std_dev;
m_threshold = threshold;
}
FilterLowIntensitySignalRemoval::FilterLowIntensitySignalRemoval(const QString &parameters)
{
// qDebug() << "Building with parameters:" << parameters;
buildFilterFromString(parameters);
}
FilterLowIntensitySignalRemoval::FilterLowIntensitySignalRemoval(
const FilterLowIntensitySignalRemoval &other)
{
m_noiseMean = other.m_noiseMean;
m_noiseStdDev = other.m_noiseStdDev;
m_threshold = other.m_threshold;
}
FilterLowIntensitySignalRemoval::~FilterLowIntensitySignalRemoval()
{
}
FilterLowIntensitySignalRemoval &
FilterLowIntensitySignalRemoval::operator=(const FilterLowIntensitySignalRemoval &other)
{
if(&other == this)
return *this;
m_noiseMean = other.m_noiseMean;
m_noiseStdDev = other.m_noiseStdDev;
m_threshold = other.m_threshold;
return *this;
}
void
FilterLowIntensitySignalRemoval::buildFilterFromString(const QString &parameters)
{
// qDebug();
// Typical string: "FilterLowIntensitySignalRemoval|2000"
if(parameters.startsWith(QString("%1|").arg(name())))
{
QStringList params = parameters.split("|").back().split(";");
if(params.size() != 3)
throw pappso::ExceptionNotRecognized(
QString("Building of FilterLowIntensitySignalRemoval from string '%1' failed. Expected "
"name and three parameters.")
.arg(parameters));
m_noiseMean = params.at(0).toDouble();
m_noiseStdDev = params.at(1).toDouble();
m_threshold = params.at(2).toDouble();
}
else
{
throw pappso::ExceptionNotRecognized(
QString("Building of FilterLowIntensitySignalRemoval from string %1 failed")
.arg(parameters));
}
}
std::size_t
FilterLowIntensitySignalRemoval::detectClusterApices(const Trace &trace)
{
// We want to iterate in the trace and check if we can get the apicess
// isolated one by one. An apex is nothing more than the top cusp of a
// triangle. We only store apices that have a value > m_threshold.
m_clusters.clear();
std::size_t trace_size = trace.size();
if(trace_size <= 2)
{
// qDebug() << "The original trace has less than 3 points. Returning it "
//"without modification.";
return m_clusters.size();
}
// else
// qDebug() << "The original trace has" << trace_size << "data points";
// Seed the system with the first point of the trace.
m_curIter = trace.cbegin();
m_prevIter = trace.cbegin();
// qDebug() << "First trace point: "
//<< "(" << m_prevIter->x << "," << m_prevIter->y << ");";
// We still have not found any apex!
m_prevApex = trace.cend();
// We will need to know if we were ascending to an apex.
bool was_ascending_to_apex = false;
// Prepare a first apex cluster.
ApicesSPtr cluster_apices = std::make_shared<ClusterApices>();
// Now that we have seeded the system with the first point of the original
// trace, go to the next one.
++m_curIter;
// And now iterate in the original trace and make sure we detect all the
// apicess.
while(m_curIter != trace.cend())
{
// qDebug() << "Current trace point: "
//<< "(" << m_curIter->x << "," << m_curIter->y << ");";
// We monitor if we are going up or down a peak.
if(m_curIter->y > m_prevIter->y)
// We are ascending a peak. We do not know if we are at the apex.
{
// qDebug().noquote() << "We are ascending to an apex.\n";
was_ascending_to_apex = true;
}
else
// We are descending a peak.
{
// qDebug().noquote() << "Descending a peak. ";
// There are two situations:
//
// 1. Either we were ascending to an apex, and m_prev is that apex,
//
// 2. Or we were not ascending to an apex and in fact all we are doing
// is going down an apex that occurred more than one trace point ago.
if(!was_ascending_to_apex)
// We were not ascending to an apex.
{
// Do nothing here.
// qDebug().noquote()
//<< "But, we were not ascending to an apex, so do nothing.\n";
}
else
// We are effectively descending a peak right after the apex was
// reached at the previous iteration.
{
m_curApex = m_prevIter;
// qDebug().noquote()
//<< "And, we were ascending to an apex, so "
//"m_curApex has become m_prevIter: ("
//<< m_curApex->x << "," << m_curApex->y << ")\n";
// We might have two situations:
// 1. We had already encountered an apex.
// 2. We had not yet encountered an apex.
if(m_prevApex != trace.cend())
// We had already encountered an apex.
{
// Was that apex far on the left of the current apex ?
if(m_curApex->x - m_prevApex->x > INTRA_CLUSTER_INTER_PEAK_DISTANCE)
// The distance is not compatible with both apices to belong
// to the same apices cluster.
{
// We are creating a new isotopic apices.
// But, since another apex had been encountered already,
// that means that an isotopic apices was cooking
// already. We must store it.
if(!cluster_apices->size())
qFatal("Cannot be that the apices has no apex.");
// Store the crafted apices cluster.
m_clusters.push_back(cluster_apices);
// qDebug().noquote()
//<< "There was a previous apex already, BUT "
//"outside of apices range. "
//"Pushing the cooking apices that has size:"
//<< cluster_apices->size();
// Now create a brand new apices for later work.
cluster_apices = std::make_shared<ClusterApices>();
// qDebug() << "Created a brand new apices.";
// We only start the new apices with the current apex if
// that apex is above the threshold.
if(m_curApex->y > m_threshold)
{
// And start it with the current apex.
cluster_apices->push_back(m_curApex);
// qDebug()
//<< "Since the current apex is above the threshold, "
//"we PUSH it to the newly created apices: ("
//<< m_curApex->x << "," << m_curApex->y << ")";
m_prevApex = m_curApex;
// qDebug() << "Set prev apex to be cur apex.";
}
else
{
// qDebug()
//<< "Since the current apex is below the threshold, "
//"we DO NOT push it to the newly created "
//"apices: ("
//<< m_curApex->x << "," << m_curApex->y << ")";
// Since previous apex went to the closed apices, we
// need to reset it.
m_prevApex = trace.cend();
// qDebug() << "Since the previous apex went to the "
//"closed apices, and cur apex has too "
//"small an intensity, we reset prev apex "
//"to trace.cend().";
}
}
else
// The distance is compatible with both apices to belong to
// the same isotopic apices.
{
// But we only push back the current apex to the apices
// if its intensity is above the threshold.
if(m_curApex->y > m_threshold)
// The current apex was above the threshold
{
cluster_apices->push_back(m_curApex);
// qDebug().noquote()
//<< "There was an apex already inside of apices "
//"range. "
//"AND, since the current apex was above the "
//"threshold, we indeed push it to apices.\n";
// qDebug().noquote()
//<< "Current apex PUSHED: " << m_curApex->x << ", "
//<< m_curApex->y;
m_prevApex = m_curApex;
// qDebug() << "We set prev apex to be cur apex.";
}
else
{
// qDebug().noquote()
//<< "There was an apex already inside of apices "
//"range. "
//"BUT, since the current apex was below the "
//"threshold, we do not push it to apices.\n";
// qDebug().noquote()
//<< "Current apex NOT pushed: " << m_curApex->x
//<< ", " << m_curApex->y;
}
}
}
else
// No apex was previously found. We are fillin-up a new isotopic
// apices.
{
if(m_curApex->y > m_threshold)
// We can actually add that apex to a new isotopic apices.
{
if(cluster_apices->size())
qCritical(
"At this point, the apices should be new and "
"empty.");
cluster_apices->push_back(m_curApex);
// qDebug().noquote()
//<< "No previous apex was found. Since current apex' "
//"intensity is above threshold, we push it back to "
//"the "
//"apices.\n";
// Store current apex as previous apex for next rounds.
m_prevApex = m_curApex;
// qDebug().noquote() << "And thus we store the current "
//"apex as previous apex.\n";
}
else
{
// qDebug().noquote()
//<< "No previous apex was found. Since current apex' "
//"intensity is below threshold, we do nothing.\n";
}
}
}
// End of
// ! if(!was_ascending_to_apex)
// That is, we were ascending to an apex.
// Tell what it is: we were not ascending to an apex.
was_ascending_to_apex = false;
}
// End of
// ! if(m_curIter->y > m_prevIter->y)
// That is, we are descending a peak.
// At this point, prepare next round.
// qDebug().noquote()
//<< "Preparing next round, with m_prevIter = m_curIter and ++index.\n";
m_prevIter = m_curIter;
++m_curIter;
}
// End of
// while(index < trace_size)
// At this point, if a apices had been cooking, add it.
if(cluster_apices->size())
m_clusters.push_back(cluster_apices);
return m_clusters.size();
}
Trace::const_iterator
FilterLowIntensitySignalRemoval::backwardFindApex(const Trace &trace,
Trace::const_iterator iter,
double distance_threshold)
{
// We receive an iterator that points to an apex. We want iterate back in the
// trace working copy and look if there are apices that are distant less than
// 1.1~Th.
Trace::const_iterator init_iter = iter;
if(iter == trace.cbegin())
return iter;
// Seed the previous iter to iter, because we'll move from there right away.
Trace::const_iterator prev_iter = init_iter;
Trace::const_iterator last_apex_iter = init_iter;
// We will need to know if we were ascending to an apex.
bool was_ascending_to_apex = false;
// Now that we have seeded the system, we can move iter one data point:
--iter;
while(iter != trace.cbegin())
{
// If we are already outside of distance_threshold, return the last apex
// that was found (or initial iter if none was encountered)
if(abs(init_iter->x - iter->x) >= distance_threshold)
return last_apex_iter;
if(iter->y > prev_iter->y)
{
// New data point has greater intensity. Just record that fact.
was_ascending_to_apex = true;
}
else
{
// New data point has smaller intensity. We are descending a peak.
if(was_ascending_to_apex)
{
// We had an apex at previous iter. We are inside the distance
// threshold. That is good. But we still could find another apex
// while moving along the trace that is in the distance threshold.
// This is why we keep going, but we store the previous iter as
// the last encountered apex.
last_apex_iter = prev_iter;
}
}
prev_iter = iter;
// Move.
--iter;
}
// qDebug() << "Init m/z: " << init_iter->x
//<< "Returning m/z: " << last_apex_iter->x
//<< "at distance:" << std::distance(last_apex_iter, init_iter);
// At this point last_apex_iter is the same as the initial iter.
return last_apex_iter;
}
Trace::const_iterator
FilterLowIntensitySignalRemoval::forwardFindApex(const Trace &trace,
Trace::const_iterator iter,
double distance_threshold)
{
// We receive an iterator that points to an apex. We want iterate back in the
// trace working copy and look if there are apices that are distant less than
// 1.1~Th.
Trace::const_iterator init_iter = iter;
if(iter == trace.cend())
return iter;
// Seed the previous iter to iter, because we'll move from there right away.
Trace::const_iterator prev_iter = init_iter;
Trace::const_iterator last_apex_iter = init_iter;
// We will need to know if we were ascending to an apex.
bool was_ascending_to_apex = false;
// Now that we have seeded the system, we can move iter one data point:
++iter;
while(iter != trace.cend())
{
// If we are already outside of distance_threshold, return the last apex
// that was found (or initial iter if none was encountered)
// FIXME: maybe we should compare prev_iter->x with iter->x so that we
// continue moving if all the succcessive apices found are each one in the
// distance_threshold from the previous one ?
if(abs(init_iter->x - iter->x) >= distance_threshold)
return last_apex_iter;
if(iter->y > prev_iter->y)
{
// New data point has greater intensity. Just record that fact.
was_ascending_to_apex = true;
}
else
{
// New data point has smaller intensity. We are descending a peak.
if(was_ascending_to_apex)
{
// We had an apex at previous iter. We are inside the distance
// threshold. That is good. But we still could find another apex
// while moving along the trace that is in the distance threshold.
// This is why we keep going, but we store the previous iter as
// the last encountered apex.
last_apex_iter = prev_iter;
}
}
prev_iter = iter;
// Move.
++iter;
}
// At this point last_apex_iter is the same as the initial iter.
return last_apex_iter;
}
Trace
FilterLowIntensitySignalRemoval::reconstructTrace(const Trace &trace)
{
// We start from the vector of apices and try to remake a high resolution
// trace out of these apices.
MapTrace map_trace;
// The general idea is that apices were detected, and only apices having
// their intensity above the threshold were stored. That means that we need to
// add points to the trace to reconstruct a meaningful trace. Indeed, imagine
// a heavy peptide isotopic cluster: the first peak's apex is below threshold
// and was not stored. The second peak is also below. But the third isotopic
// cluster peak is above and was stored.
//
// How do we reconstruct the trace to have all these points that were
// preceding the first isotopic cluster apex that was detected.
//
// The same applies to the last peaks of the cluster that are below the
// threshold whil the preceeding ones were above!
// The general idea is to iterate in the vector of apices and for each apex
// that is encountered, ask if there were apices
// m_clusters is a vector that contains vectors of Trace::const_iter.
// Each vector in m_clusters should represent all the apices of a cluster.
// qDebug() << "Reconstructing trace with " << m_clusters.size() <<
// "clusters.";
Trace::const_iterator left_begin_iter = trace.cend();
Trace::const_iterator right_end_iter = trace.cend();
for(auto &&cluster_apices : m_clusters)
{
// cluster_apices is a vector of Trace::const_iterator. If we want to
// reconstruct the Trace, we need to iterate through all the DataPoint
// objects in between cluster_apices.begin() and cluster_apices.end().
Trace::const_iterator begin_iter = *(cluster_apices->begin());
Trace::const_iterator end_iter = *(std::prev(cluster_apices->end()));
// qDebug() << "Iterating in a new cluster apices with boundaries:"
//<< begin_iter->x << "-" << end_iter->x;
left_begin_iter = backwardFindApex(trace, begin_iter, INTRA_CLUSTER_INTER_PEAK_DISTANCE);
// qDebug() << "After backwardFindApex, left_begin_iter points to:"
//<< left_begin_iter->toString() << "with distance:"
//<< std::distance(left_begin_iter, begin_iter);
right_end_iter = forwardFindApex(trace, end_iter, 1.5 * INTRA_CLUSTER_INTER_PEAK_DISTANCE);
for(Trace::const_iterator iter = left_begin_iter; iter != right_end_iter; ++iter)
{
map_trace[iter->x] = iter->y;
}
// Now reset the left and right iters to intensity 0 to avoid having
// disgraceful oblique lines connnecting trace segments.
map_trace[left_begin_iter->x] = 0;
map_trace[std::prev(right_end_iter)->x] = 0;
}
return map_trace.toTrace();
}
Trace &
FilterLowIntensitySignalRemoval::filter(Trace &trace) const
{
// qDebug();
// Horrible hack to have a non const filtering process.
return const_cast<FilterLowIntensitySignalRemoval *>(this)->nonConstFilter(trace);
}
Trace &
FilterLowIntensitySignalRemoval::nonConstFilter(Trace &trace)
{
// qDebug();
if(trace.size() <= 2)
{
// qDebug() << "The original trace has less than 3 points. Returning it "
//"without modification.";
return trace;
}
// else
// qDebug() << "The original trace had" << trace.size() << "data points";
std::size_t cluster_count = detectClusterApices(trace);
// qDebug() << "Number of detected cluster apices: " << cluster_count;
if(!cluster_count)
return trace;
// At this point we want to work on the apices and reconstruct a full
// trace.
Trace reconstructed_trace = reconstructTrace(trace);
trace = std::move(reconstructed_trace);
return trace;
}
double
FilterLowIntensitySignalRemoval::getThreshold() const
{
return m_threshold;
}
//! Return a string with the textual representation of the configuration data.
QString
FilterLowIntensitySignalRemoval::toString() const
{
return QString("%1|%2;%3;%4")
.arg(name())
.arg(QString::number(m_noiseMean, 'f', 2))
.arg(QString::number(m_noiseStdDev, 'f', 2))
.arg(QString::number(m_threshold, 'f', 2));
}
QString
FilterLowIntensitySignalRemoval::name() const
{
return "FilterLowIntensitySignalRemoval";
}
} // namespace pappso
pappso::ExceptionNotRecognized
excetion to use when an item type is not recognized
Definition exceptionnotrecognized.h:42
pappso::FilterLowIntensitySignalRemoval
Redefines the floor intensity of the Trace.
Definition filterlowintensitysignalremoval.h:73
pappso::FilterLowIntensitySignalRemoval::m_curApex
TraceCIter m_curApex
Definition filterlowintensitysignalremoval.h:133
pappso::FilterLowIntensitySignalRemoval::nonConstFilter
Trace & nonConstFilter(Trace &data_points)
Definition filterlowintensitysignalremoval.cpp:614
pappso::FilterLowIntensitySignalRemoval::ApicesSPtr
std::shared_ptr< ClusterApices > ApicesSPtr
Definition filterlowintensitysignalremoval.h:138
pappso::FilterLowIntensitySignalRemoval::backwardFindApex
Trace::const_iterator backwardFindApex(const Trace &trace, Trace::const_iterator iter, double distance_threshold)
Definition filterlowintensitysignalremoval.cpp:397
pappso::FilterLowIntensitySignalRemoval::m_curIter
TraceCIter m_curIter
Definition filterlowintensitysignalremoval.h:134
pappso::FilterLowIntensitySignalRemoval::name
QString name() const override
Definition filterlowintensitysignalremoval.cpp:665
pappso::FilterLowIntensitySignalRemoval::FilterLowIntensitySignalRemoval
FilterLowIntensitySignalRemoval(double mean, double std_dev, double threshold)
Definition filterlowintensitysignalremoval.cpp:52
pappso::FilterLowIntensitySignalRemoval::getThreshold
double getThreshold() const
Definition filterlowintensitysignalremoval.cpp:646
pappso::FilterLowIntensitySignalRemoval::m_noiseMean
double m_noiseMean
Definition filterlowintensitysignalremoval.h:113
pappso::FilterLowIntensitySignalRemoval::reconstructTrace
Trace reconstructTrace(const Trace &trace)
Definition filterlowintensitysignalremoval.cpp:535
pappso::FilterLowIntensitySignalRemoval::filter
Trace & filter(Trace &data_points) const override
Definition filterlowintensitysignalremoval.cpp:604
pappso::FilterLowIntensitySignalRemoval::m_prevApex
TraceCIter m_prevApex
Definition filterlowintensitysignalremoval.h:132
pappso::FilterLowIntensitySignalRemoval::forwardFindApex
Trace::const_iterator forwardFindApex(const Trace &trace, Trace::const_iterator iter, double distance_threshold)
Definition filterlowintensitysignalremoval.cpp:466
pappso::FilterLowIntensitySignalRemoval::m_threshold
double m_threshold
Definition filterlowintensitysignalremoval.h:112
pappso::FilterLowIntensitySignalRemoval::detectClusterApices
std::size_t detectClusterApices(const Trace &trace)
Definition filterlowintensitysignalremoval.cpp:131
pappso::FilterLowIntensitySignalRemoval::operator=
FilterLowIntensitySignalRemoval & operator=(const FilterLowIntensitySignalRemoval &other)
Definition filterlowintensitysignalremoval.cpp:88
pappso::FilterLowIntensitySignalRemoval::INTRA_CLUSTER_INTER_PEAK_DISTANCE
static constexpr double INTRA_CLUSTER_INTER_PEAK_DISTANCE
Definition filterlowintensitysignalremoval.h:116
pappso::FilterLowIntensitySignalRemoval::~FilterLowIntensitySignalRemoval
virtual ~FilterLowIntensitySignalRemoval()
Definition filterlowintensitysignalremoval.cpp:82
pappso::FilterLowIntensitySignalRemoval::m_noiseStdDev
double m_noiseStdDev
Definition filterlowintensitysignalremoval.h:114
pappso::FilterLowIntensitySignalRemoval::m_prevIter
TraceCIter m_prevIter
Definition filterlowintensitysignalremoval.h:135
pappso::FilterLowIntensitySignalRemoval::buildFilterFromString
void buildFilterFromString(const QString &strBuildParams) override
build this filter using a string
Definition filterlowintensitysignalremoval.cpp:102
pappso::FilterLowIntensitySignalRemoval::toString
QString toString() const override
Return a string with the textual representation of the configuration data.
Definition filterlowintensitysignalremoval.cpp:654
pappso::FilterLowIntensitySignalRemoval::m_clusters
std::vector< ApicesSPtr > m_clusters
Definition filterlowintensitysignalremoval.h:141
pappso::MapTrace
Definition maptrace.h:36
pappso::MapTrace::toTrace
Trace toTrace() const
Definition maptrace.cpp:214
pappso::Trace
A simple container of DataPoint instances.
Definition trace.h:152
exceptionnotrecognized.h
filterlowintensitysignalremoval.h
maptrace.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39