massspectrum.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/massspectrum/massspectrum.cpp
 * \date 15/3/2015
 * \author Olivier Langella
 * \brief basic mass spectrum
 */
 
/*******************************************************************************
 * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     PAPPSOms++ 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.
 *
 *     PAPPSOms++ 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 PAPPSOms++.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Contributors:
 *     Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
 *implementation
 ******************************************************************************/
 
#include <cmath>
#include <numeric>
#include <limits>
#include <vector>
#include <map>
 
#include <QDebug>
 
#include "pappsomspp/core/trace/datapoint.h"
#include "pappsomspp/core/trace/trace.h"
#include "pappsomspp/core/massspectrum/massspectrum.h"
#include "pappsomspp/core/processing/combiners/massspectrumcombiner.h"
#include "pappsomspp/core/mzrange.h"
#include "pappsomspp/core/pappsoexception.h"
 
#include "pappsomspp/core/peptide/peptidefragmentionlistbase.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "pappsomspp/core/processing/filters/filterresample.h"
 
 
int massSpectrumMetaTypeId =
  qRegisterMetaType<pappso::MassSpectrum>("pappso::MassSpectrum");
int massSpectrumPtrMetaTypeId =
  qRegisterMetaType<pappso::MassSpectrum *>("pappso::MassSpectrum *");
 
 
namespace pappso
{
 
 
MassSpectrum::MassSpectrum()
{
}
 
 
MassSpectrum::MassSpectrum(
  std::vector<std::pair<pappso_double, pappso_double>> &data_point_vector)
  : Trace::Trace(data_point_vector)
{
}
 
MassSpectrum::MassSpectrum(std::vector<DataPoint> &data_point_vector)
  : Trace::Trace(data_point_vector)
{
}
 
MassSpectrum::MassSpectrum(const Trace &other) : Trace(other)
{
}
 
MassSpectrum::MassSpectrum(const MapTrace &other) : Trace(other)
{
}
 
 
MassSpectrum::MassSpectrum(Trace &&other) : Trace(std::move(other))
{
}
 
 
MassSpectrum::MassSpectrum(const MassSpectrum &other) : Trace(other)
{
  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()";
}
 
 
MassSpectrum::MassSpectrum(MassSpectrum &&other) : Trace(std::move(other))
{
  // Specify std::move so that && reference is passed to the Trace constructor
  // that takes std::vector<DataPoint> && as rvalue reference.
 
  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()"
  //<< "Moving MassSpectrum::MassSpectrum(MassSpectrum &&)";
}
 
 
MassSpectrum::~MassSpectrum()
{
}
 
 
MassSpectrum &
MassSpectrum::operator=(const MassSpectrum &other)
{
  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()";
 
  assign(other.begin(), other.end());
  return *this;
}
 
 
MassSpectrum &
MassSpectrum::operator=(MassSpectrum &&other)
{
  vector<DataPoint>::operator=(std::move(other));
  return *this;
}
 
 
MassSpectrumSPtr
MassSpectrum::makeMassSpectrumSPtr() const
{
  return std::make_shared<MassSpectrum>(*this);
}
 
 
MassSpectrumCstSPtr
MassSpectrum::makeMassSpectrumCstSPtr() const
{
  return std::make_shared<const MassSpectrum>(*this);
}
 
 
//! Compute the total ion current of this mass spectrum
/*!
 * The sum of all the separate ion currents carried by the ions of different
 * m/z contributing to a complete mass massSpectrum or in a specified m/z
 * range of a mass massSpectrum. MS:1000285
 *
 * \return <pappso_double> The total ion current.
 */
pappso_double
MassSpectrum::totalIonCurrent() const
{
  return Trace::sumY();
}
 
 
//! Compute the total ion current of this mass spectrum
/*!
 * Convenience function that returns totalIonCurrent();
 */
pappso_double
MassSpectrum::tic() const
{
  return totalIonCurrent();
}
 
 
pappso_double
MassSpectrum::tic(double mzStart, double mzEnd)
{
  return Trace::sumY(mzStart, mzEnd);
}
 
 
//! Find the DataPoint instance having the greatest intensity (y) value.
/*!
 * \return <const DataPoint &> The data point having the maximum intensity (y)
 * value of the whole mass spectrum.
 */
const DataPoint &
MassSpectrum::maxIntensityDataPoint() const
{
  return Trace::maxYDataPoint();
}
 
 
//! Find the DataPoint instance having the smallest intensity (y) value.
/*!
 * \return <const DataPoint &> The data point having the minimum intensity (y)
 * value of the whole mass spectrum.
 */
const DataPoint &
MassSpectrum::minIntensityDataPoint() const
{
  return Trace::minYDataPoint();
}
 
 
//! Sort the DataPoint instances of this spectrum.
/*!
 * The DataPoint instances are sorted according to the x value (the m/z
 * value) and in increasing order.
 */
void
MassSpectrum::sortMz()
{
  Trace::sortX();
}
 
 
//! Tells if \c this MassSpectrum is equal to \p massSpectrum.
/*!
 * To compare \c this to \p massSpectrum, a tolerance is applied to both the x
 * and y values, that is defined using \p precision.
 *
 * \param massSpectrum Mass spectrum to compare to \c this.
 *
 * \param precision Precision to be used to perform the comparison of the x
 * and y values of the data points in \c this and \massSpectrum mass spectra.
 */
bool
MassSpectrum::equals(const MassSpectrum &other, PrecisionPtr precision) const
{
  if(size() != other.size())
    {
      qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()"
               << "The other mass spectrum size is not equal to *this size"
               << "*this size:" << size() << "trace size:" << other.size();
 
      return false;
    }
 
  PrecisionPtr precint = PrecisionFactory::getPpmInstance(10);
 
  auto trace_it = other.begin();
 
  for(auto &&data_point : *this)
    {
      qDebug() << "first:" << data_point.x << "," << data_point.y
               << " second:" << trace_it->x << "," << trace_it->y;
 
      if(!MzRange(data_point.x, precision).contains(trace_it->x))
        {
          qDebug() << "x:" << data_point.x << " != " << trace_it->x;
          return false;
        }
 
      if(!MzRange(data_point.y, precint).contains(trace_it->y))
        {
          qDebug() << "y:" << data_point.y << " != " << trace_it->y;
          return false;
        }
 
      trace_it++;
    }
 
  return true;
}
 
 
MassSpectrum
MassSpectrum::filterSum(const MzRange &range) const
{
  MassSpectrum massSpectrum;
 
  std::vector<DataPoint>::const_iterator it    = begin();
  std::vector<DataPoint>::const_iterator itEnd = end();
 
  std::vector<DataPoint>::const_reverse_iterator itRev    = rbegin();
  std::vector<DataPoint>::const_reverse_iterator itRevEnd = rend();
 
  pappso_double lower = range.lower();
  pappso_double upper = range.upper();
 
  while((it != itEnd) && (it->x <= itRev->x) && (itRev != itRevEnd))
    {
      pappso_double sumX = it->x + itRev->x;
 
      if(sumX < lower)
        {
          it++;
        }
      else if(sumX > upper)
        {
          itRev++;
        }
      else
        {
          massSpectrum.push_back(*it);
          massSpectrum.push_back(*itRev);
 
          std::vector<DataPoint>::const_reverse_iterator itRevIn = itRev;
          itRevIn++;
 
          // FIXME Attention buggy code FR 20180626.
          sumX = it->x + itRevIn->x;
          while((sumX > lower) && (it->x <= itRevIn->x) &&
                (itRevIn != itRevEnd))
            {
              sumX = it->x + itRevIn->x;
              // trace.push_back(*it);
              massSpectrum.push_back(*itRevIn);
              itRevIn++;
            }
          it++;
        }
    }
 
  // Sort all the data points in increasing order by x
  std::sort(massSpectrum.begin(),
            massSpectrum.end(),
            [](const DataPoint &a, const DataPoint &b) { return (a.x < b.x); });
 
  // Remove all the but the first element of a series of elements that are
  // considered equal. Sort of deduplication.
  std::vector<DataPoint>::iterator itEndFix =
    std::unique(massSpectrum.begin(),
                massSpectrum.end(),
                [](const DataPoint &a, const DataPoint &b) {
                  // Return true if both elements should be considered equal.
                  return (a.x == b.x) && (a.y == b.y);
                });
 
  massSpectrum.resize(std::distance(massSpectrum.begin(), itEndFix));
 
  return massSpectrum;
}
 
 
void
MassSpectrum::debugPrintValues() const
{
 
  qDebug() << size();
  for(std::size_t i = 0; i < size(); i++)
    {
      qDebug() << "(" << this->operator[](i).x << "," << this->operator[](i).y
               << ")";
    }
}
 
 
QDataStream &
operator<<(QDataStream &outstream, const MassSpectrum &massSpectrum)
{
  quint32 vector_size = massSpectrum.size();
  outstream << vector_size;
  for(auto &&peak : massSpectrum)
    {
      outstream << peak;
    }
 
  return outstream;
}
 
 
QDataStream &
operator>>(QDataStream &instream, MassSpectrum &massSpectrum)
{
 
  quint32 vector_size;
  DataPoint peak;
 
  if(!instream.atEnd())
    {
      instream >> vector_size;
      qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__
               << " vector_size=" << vector_size;
 
      for(quint32 i = 0; i < vector_size; i++)
        {
 
          if(instream.status() != QDataStream::Ok)
            {
              throw PappsoException(
                QString("error in QDataStream unserialize operator>> of "
                        "massSpectrum :\nread datastream failed status=%1 "
                        "massSpectrum "
                        "i=%2 on size=%3")
                  .arg(instream.status())
                  .arg(i)
                  .arg(vector_size));
            }
          instream >> peak;
          massSpectrum.push_back(peak);
        }
      if(instream.status() != QDataStream::Ok)
        {
          throw PappsoException(
            QString(
              "error in QDataStream unserialize operator>> of massSpectrum "
              ":\nread datastream failed status=%1")
              .arg(instream.status()));
        }
    }
 
  return instream;
}
 
 
MassSpectrum &
MassSpectrum::massSpectrumFilter(const MassSpectrumFilterInterface &filter)
{
  return filter.filter(*this);
}
 
} // namespace pappso