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