filterchargedeconvolution.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/filters/filtertdfcorrectpeak.cpp
 * \date 30/09/2020
 * \author Thomas Renne
 * \brief Sum peaks and transform mz to fit charge = 1
 */
 
/*******************************************************************************
 * Copyright (c) 2020 Olivier Langella <Olivier.Langella@u-psud.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/>.
 *
 ******************************************************************************/
#include "filterchargedeconvolution.h"
#include <QDebug>
#include "pappsomspp/core/exception/exceptionnotrecognized.h"
#include "pappsomspp/core/mzrange.h"
 
using namespace pappso;
 
pappso::FilterChargeDeconvolution::FilterChargeDeconvolution(const QString &strBuildParams)
{
  buildFilterFromString(strBuildParams);
  // qInfo() << "ChargeDeconvolution created";
}
 
 
FilterChargeDeconvolution::FilterChargeDeconvolution(PrecisionPtr precision_ptr)
  : m_precisionPtrZ1(precision_ptr)
{
  m_diffC12C13_z1  = DIFFC12C13;
  m_diffC12C13_z2  = DIFFC12C13 / 2;
  m_precisionPtrZ1 = precision_ptr;
  m_precisionPtrZ2 =
    pappso::PrecisionFactory::getPrecisionPtrFractionInstance(m_precisionPtrZ1, 0.5);
 
  // pappso::PrecisionFactory::getPrecisionDividedBy(m_precisionPtrZ1, 2):
 
  // qInfo() << m_precisionPtrZ2->getNominal();
}
 
FilterChargeDeconvolution::FilterChargeDeconvolution(const FilterChargeDeconvolution &other)
  : m_precisionPtrZ1(other.m_precisionPtrZ1), m_precisionPtrZ2(other.m_precisionPtrZ2)
{
  m_diffC12C13_z1 = DIFFC12C13;
  m_diffC12C13_z2 = DIFFC12C13 / 2;
}
 
FilterChargeDeconvolution::~FilterChargeDeconvolution()
{
  qDebug() << "ChargeDeconvolution destroyed";
}
 
 
void
pappso::FilterChargeDeconvolution::buildFilterFromString(const QString &strBuildParams)
{
  //"chargeDeconvolution|0.02dalton"
  qDebug();
  if(strBuildParams.startsWith("chargeDeconvolution|"))
    {
      QStringList params = strBuildParams.split("|").back().split(";", Qt::SkipEmptyParts);
 
      QString precision = params.at(0);
      m_precisionPtrZ1  = PrecisionFactory::fromString(
        precision.replace("dalton", " dalton").replace("ppm", " ppm").replace("res", " res"));
      qDebug();
      m_precisionPtrZ2 =
        pappso::PrecisionFactory::getPrecisionPtrFractionInstance(m_precisionPtrZ1, 0.5);
 
 
      m_diffC12C13_z1 = DIFFC12C13;
      m_diffC12C13_z2 = DIFFC12C13 / 2;
    }
  else
    {
      throw pappso::ExceptionNotRecognized(
        QString("building chargeDeconvolution from string %1 is not possible").arg(strBuildParams));
    }
  qDebug();
}
 
 
QString
pappso::FilterChargeDeconvolution::name() const
{
  return "chargeDeconvolution";
}
 
 
QString
pappso::FilterChargeDeconvolution::toString() const
{
  QString strCode = QString("%1|%2").arg(name()).arg(m_precisionPtrZ1->toString());
 
  strCode.replace(" ", "");
 
  return strCode;
}
 
Trace &
FilterChargeDeconvolution::filter(Trace &data_points) const
{
  qDebug();
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> data_points_info;
  data_points.sortY(Enums::SortOrder::descending);
  qDebug() << data_points.size();
  Trace new_trace;
 
  for(auto &data_point : data_points)
    {
      addDataPointToList(data_points_info, data_point);
    }
  computeBestChargeOfDataPoint(data_points_info);
 
  //   qDebug() << data_points_info.size();
  computeIsotopeDeconvolution(data_points_info);
  //   qDebug() << data_points_info.size();
  transformToMonoChargedForAllDataPoint(data_points_info);
  for(DataPointInfoSp &dpi : data_points_info)
    {
      //       qDebug() << dpi->new_mono_charge_data_point.x << dpi->z_charge;
      new_trace.push_back(dpi->new_mono_charge_data_point);
    }
 
  new_trace.sortX();
  data_points = std::move(new_trace);
  qDebug() << data_points.size();
  qDebug();
  return data_points;
}
 
void
pappso::FilterChargeDeconvolution::addDataPointToList(
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> &points,
  pappso::DataPoint &data_point) const
{
  DataPointInfoSp new_dpi(std::make_shared<DataPointInfo>());
 
  new_dpi->data_point = data_point;
  MzRange range1(data_point.x + m_diffC12C13_z1, m_precisionPtrZ1);
  new_dpi->z1_range = std::pair<double, double>(range1.lower(), range1.upper());
  MzRange range2(data_point.x + m_diffC12C13_z2, m_precisionPtrZ2);
  new_dpi->z2_range = std::pair<double, double>(range2.lower(), range2.upper());
  addDataPointRefByExclusion(points, new_dpi);
  points.push_back(new_dpi);
}
 
void
pappso::FilterChargeDeconvolution::addDataPointRefByExclusion(
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> &points,
  FilterChargeDeconvolution::DataPointInfoSp &new_dpi) const
{
  // add datapoint which match the mz_range = 1 to z1_list
  auto i_z1 = points.begin(), end = points.end();
  while(i_z1 != end)
    {
      // get the datapoint which match the range
      i_z1 = std::find_if(i_z1, end, [&new_dpi](DataPointInfoSp dpi) {
        return (new_dpi->data_point.x >= dpi->z1_range.first &&
                new_dpi->data_point.x <= dpi->z1_range.second);
      });
      if(i_z1 != end)
        {
          // add the datapoint to the list and add the parent pointer
          i_z1->get()->z1_vect.push_back(new_dpi);
          new_dpi->parent           = *i_z1;
          DataPointInfoSp parent_z1 = i_z1->get()->parent.lock();
          while(parent_z1 != nullptr)
            {
              parent_z1.get()->z1_vect.push_back(new_dpi);
              parent_z1 = parent_z1->parent.lock();
            }
          i_z1++;
        }
    }
 
  // add datapoint which match the mz_range = 2 to z2_list
  auto i_z2 = points.begin();
  while(i_z2 != end)
    {
      // get the datapoint which match the range
      i_z2 = std::find_if(i_z2, end, [&new_dpi](DataPointInfoSp dpi) {
        return (new_dpi->data_point.x >= dpi->z2_range.first &&
                new_dpi->data_point.x <= dpi->z2_range.second);
      });
      if(i_z2 != end)
        {
          // add the datapoint to the list and add the parent pointer
          i_z2->get()->z2_vect.push_back(new_dpi);
          new_dpi->parent           = *i_z2;
          DataPointInfoSp parent_z2 = i_z2->get()->parent.lock();
          while(parent_z2 != nullptr)
            {
              parent_z2.get()->z2_vect.push_back(new_dpi);
              parent_z2 = parent_z2->parent.lock();
            }
          i_z2++;
        }
    }
}
 
void
pappso::FilterChargeDeconvolution::computeBestChargeOfDataPoint(
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> &data_points_info) const
{
  for(DataPointInfoSp &data_point_info : data_points_info)
    {
      if(data_point_info.get()->z1_vect.size() >= 1 && data_point_info.get()->z2_vect.size() == 0)
        {
          for(std::weak_ptr<DataPointInfo> other : data_point_info.get()->z1_vect)
            {
              other.lock()->z_charge = 1;
            }
          data_point_info.get()->z_charge = 1;
        }
      else if(data_point_info.get()->z1_vect.size() == 0 &&
              data_point_info.get()->z2_vect.size() >= 1)
        {
          for(std::weak_ptr<DataPointInfo> other : data_point_info.get()->z2_vect)
            {
              other.lock()->z_charge = 2;
            }
          data_point_info.get()->z_charge = 2;
        }
      else if(data_point_info.get()->z1_vect.size() >= 1 &&
              data_point_info.get()->z2_vect.size() >= 1)
        {
          for(std::weak_ptr<DataPointInfo> other : data_point_info.get()->z2_vect)
            {
              other.lock()->z_charge = 2;
            }
          data_point_info.get()->z_charge = 2;
        }
      else
        {
          if(data_point_info.get()->z_charge == -1)
            {
              data_point_info.get()->z_charge = 0;
            }
        }
    }
}
 
void
pappso::FilterChargeDeconvolution::computeIsotopeDeconvolution(
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> &data_points_info) const
{
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> deconvoluted_points_info;
 
  for(DataPointInfoSp &data_point_info : data_points_info)
    {
      if(data_point_info->parent.lock() == nullptr)
        {
          DataPointInfoSp deconvoluted_point(std::make_shared<DataPointInfo>());
 
          deconvoluted_point->z_charge                   = data_point_info->z_charge;
          deconvoluted_point->new_mono_charge_data_point = data_point_info->data_point;
 
          if(data_point_info->z_charge == 1)
            {
 
              for(std::weak_ptr<DataPointInfo> data : data_point_info->z1_vect)
                {
                  deconvoluted_point->new_mono_charge_data_point.y += data.lock()->data_point.y;
                }
            }
          else if(data_point_info->z_charge == 2)
            {
              for(std::weak_ptr<DataPointInfo> data : data_point_info->z2_vect)
                {
                  deconvoluted_point->new_mono_charge_data_point.y += data.lock()->data_point.y;
                }
            }
          else // if z.charge == 0
            {
              deconvoluted_point->new_mono_charge_data_point = data_point_info->data_point;
            }
          deconvoluted_points_info.push_back(deconvoluted_point);
        }
    }
  data_points_info = deconvoluted_points_info;
}
 
void
pappso::FilterChargeDeconvolution::transformToMonoChargedForAllDataPoint(
  std::vector<FilterChargeDeconvolution::DataPointInfoSp> &data_points_info) const
{
  for(DataPointInfoSp &dpi : data_points_info)
    {
      if(dpi->z_charge == 2)
        {
          dpi->new_mono_charge_data_point.x += dpi->new_mono_charge_data_point.x - MHPLUS;
        }
    }
}