peptideisotopespectrummatch.cpp

Go to the documentation of this file.

/*
 * *******************************************************************************
 * * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.fr>.
 * *
 * * This file is part of MassChroqPRM.
 * *
 * *     MassChroqPRM 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.
 * *
 * *     MassChroqPRM 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 MassChroqPRM.  If not, see <http://www.gnu.org/licenses/>.
 * *
 * * Contributors:
 * *     Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
 * implementation
 * ******************************************************************************/
 
 
#include "peptideisotopespectrummatch.h"
#include "../peptide/peptidenaturalisotopelist.h"
 
namespace pappso
{
PeptideIsotopeSpectrumMatch::PeptideIsotopeSpectrumMatch(
  const MassSpectrum &spectrum,
  const PeptideSp &peptideSp,
  unsigned int parent_charge,
  PrecisionPtr precision,
  const std::list<PeptideIon> &ion_type_list_in,
  unsigned int max_isotope_number,
  [[maybe_unused]] unsigned int max_isotope_rank)
  : _precision(precision)
{
 
  try
    {
      _peak_ion_match_list.clear();
      qDebug() << "peptideSp.get()->getSequence()=" << peptideSp.get()->getSequence()
               << " max_isotope_number=" << max_isotope_number
               << " spectrum.size=" << spectrum.size() << " parent_charge=" << parent_charge
               << " ion_type_list.size=" << ion_type_list_in.size()
               << " max_isotope_rank=" << max_isotope_rank;
      PeptideFragmentIonListBase fragmentIonList(peptideSp, ion_type_list_in);
      std::list<PeptideIon> ion_type_list(ion_type_list_in);
 
      ion_type_list.sort([](const PeptideIon &iona, const PeptideIon &ionb) {
        if(iona == PeptideIon::b)
          return true;
        if(ionb == PeptideIon::b)
          return false;
        if(iona == PeptideIon::y)
          return true;
        if(ionb == PeptideIon::y)
          return false;
        return iona > ionb;
      });
 
      qDebug() << " fragmentIonList.size()=" << fragmentIonList.size();
      std::vector<DataPoint> peak_list(spectrum.begin(), spectrum.end());
      for(auto ion_type : ion_type_list)
        {
          const std::list<PeptideFragmentIonSp> peptide_fragment_ion_list =
            fragmentIonList.getPeptideFragmentIonSp(ion_type);
          qDebug() << " peptide_fragment_ion_list.size()=" << peptide_fragment_ion_list.size();
 
          for(unsigned int charge = 1; charge <= parent_charge; charge++)
            {
              qDebug();
              for(auto &&peptide_fragment_ion : peptide_fragment_ion_list)
                {
                  qDebug() << "peptide_fragment_ion.get()->getMz(charge)="
                           << peptide_fragment_ion.get()->getMz(charge);
                  PeptideNaturalIsotopeList isotope_list(peptide_fragment_ion);
                  for(unsigned int isotope_number = 0; isotope_number <= max_isotope_number;
                      isotope_number++)
                    {
                      qDebug();
                      PeptideNaturalIsotopeAverage *p_isotopeIon;
                      if(max_isotope_rank == 0)
                        {
                          // fast constructor, no need to compute isotope ratio
                          p_isotopeIon = new PeptideNaturalIsotopeAverage(
                            peptide_fragment_ion, isotope_number, charge, precision);
                        }
                      else
                        {
                          p_isotopeIon = new PeptideNaturalIsotopeAverage(
                            isotope_list, max_isotope_rank, isotope_number, charge, precision);
                        }
 
                      qDebug() << "max_isotope_rank=" << max_isotope_rank
                               << " isotope_number=" << isotope_number << " charge=" << charge
                               << " precision=" << precision
                               << " p_isotopeIon->getMz()=" << p_isotopeIon->getMz() << " "
                               << isotope_number << " " << p_isotopeIon->toString();
 
                      std::vector<DataPoint>::iterator it_peak =
                        getBestPeakIterator(peak_list, *p_isotopeIon);
                      if(it_peak != peak_list.end())
                        {
                          _peak_ion_match_list.push_back(
                            PeakIonIsotopeMatch(*it_peak,
                                                p_isotopeIon->makePeptideNaturalIsotopeAverageSp(),
                                                peptide_fragment_ion));
                          peak_list.erase(it_peak);
 
                          qDebug()
                            << isotope_number << " " << _peak_ion_match_list.back().toString();
                        }
                      delete p_isotopeIon;
                    }
                }
            }
        }
    }
  catch(PappsoException &exception_pappso)
    {
      QString errorStr =
        QObject::tr("ERROR building PeptideIsotopeSpectrumMatch, PAPPSO exception:\n%1")
          .arg(exception_pappso.qwhat());
      qDebug() << "PeptideIsotopeSpectrumMatch::PeptideIsotopeSpectrumMatch "
                  "PappsoException :\n"
               << errorStr;
      throw PappsoException(errorStr);
    }
  catch(std::exception &exception_std)
    {
      QString errorStr =
        QObject::tr("ERROR building PeptideIsotopeSpectrumMatch, std exception:\n%1")
          .arg(exception_std.what());
      qDebug() << "PeptideIsotopeSpectrumMatch::PeptideIsotopeSpectrumMatch "
                  "std::exception :\n"
               << errorStr;
      throw PappsoException(errorStr);
    }
 
  qDebug() << "_peak_ion_match_list.size()=" << _peak_ion_match_list.size();
}
 
PeptideIsotopeSpectrumMatch::PeptideIsotopeSpectrumMatch(
  const MassSpectrum &spectrum,
  std::vector<PeptideNaturalIsotopeAverageSp> v_peptideIsotopeList,
  std::vector<PeptideFragmentIonSp> v_peptideIonList,
  PrecisionPtr precision)
  : _precision(precision)
{
  qDebug() << " begin";
  if(v_peptideIsotopeList.size() != v_peptideIonList.size())
    {
      throw PappsoException(
        QObject::tr("v_peptideIsotopeList.size() %1 != v_peptideIonList.size() %2")
          .arg(v_peptideIsotopeList.size())
          .arg(v_peptideIonList.size()));
    }
 
  auto isotopeIt = v_peptideIsotopeList.begin();
  auto ionIt     = v_peptideIonList.begin();
  std::vector<DataPoint> peak_list(spectrum.begin(), spectrum.end());
 
  while(isotopeIt != v_peptideIsotopeList.end())
    {
      std::vector<DataPoint>::iterator it_peak =
        getBestPeakIterator(peak_list, *(isotopeIt->get()));
      if(it_peak != peak_list.end())
        {
          _peak_ion_match_list.push_back(PeakIonIsotopeMatch(*it_peak, *isotopeIt, *ionIt));
          peak_list.erase(it_peak);
        }
      isotopeIt++;
      ionIt++;
    }
  qDebug() << "_peak_ion_match_list.size()=" << _peak_ion_match_list.size() << " end";
}
 
 
PeptideIsotopeSpectrumMatch::PeptideIsotopeSpectrumMatch(const PeptideIsotopeSpectrumMatch &other)
  : _precision(other._precision), _peak_ion_match_list(other._peak_ion_match_list)
{
  qDebug();
}
 
PeptideIsotopeSpectrumMatch::~PeptideIsotopeSpectrumMatch()
{
}
 
 
std::vector<DataPoint>::iterator
PeptideIsotopeSpectrumMatch::getBestPeakIterator(std::vector<DataPoint> &peak_list,
                                                 const PeptideNaturalIsotopeAverage &ion) const
{
  // qDebug();
  std::vector<DataPoint>::iterator itpeak   = peak_list.begin();
  std::vector<DataPoint>::iterator itend    = peak_list.end();
  std::vector<DataPoint>::iterator itselect = peak_list.end();
 
  pappso_double best_intensity = 0;
 
  while(itpeak != itend)
    {
      // qDebug() << itpeak->x;
      if(ion.matchPeak(itpeak->x))
        {
          if(itpeak->y > best_intensity)
            {
              best_intensity = itpeak->y;
              itselect       = itpeak;
            }
        }
      itpeak++;
    }
  // qDebug();
  return (itselect);
}
 
const std::list<PeakIonIsotopeMatch> &
PeptideIsotopeSpectrumMatch::getPeakIonIsotopeMatchList() const
{
  qDebug() << "_peak_ion_match_list.size()=" << _peak_ion_match_list.size();
  return _peak_ion_match_list;
}
 
std::size_t
PeptideIsotopeSpectrumMatch::size() const
{
  return _peak_ion_match_list.size();
}
PeptideIsotopeSpectrumMatch::const_iterator
PeptideIsotopeSpectrumMatch::begin() const
{
  return _peak_ion_match_list.begin();
}
PeptideIsotopeSpectrumMatch::const_iterator
PeptideIsotopeSpectrumMatch::end() const
{
  return _peak_ion_match_list.end();
}
 
void
PeptideIsotopeSpectrumMatch::dropPeaksLackingMonoisotope()
{
  qDebug();
  _peak_ion_match_list.sort([](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
    if(a.getPeptideIonType() < b.getPeptideIonType())
      return true;
    if(a.getPeptideFragmentIonSp().get()->size() < b.getPeptideFragmentIonSp().get()->size())
      return true;
    if(a.getCharge() < b.getCharge())
      return true;
    if(a.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() <
       b.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber())
      return true;
    return false;
  });
  PeptideIon ion_type      = PeptideIon::b;
  std::size_t nserie       = 0;
  std::size_t isotopeserie = 0;
  unsigned int charge      = 0;
  for(std::list<PeakIonIsotopeMatch>::iterator it = _peak_ion_match_list.begin();
      it != _peak_ion_match_list.end();
      it++)
    {
      if((nserie != it->getPeptideFragmentIonSp().get()->size()) ||
         (ion_type != it->getPeptideIonType()) || (charge != it->getCharge()))
        {
          ion_type     = it->getPeptideIonType();
          isotopeserie = 0;
          nserie       = it->getPeptideFragmentIonSp().get()->size();
          charge       = it->getCharge();
        }
      if(isotopeserie <= it->getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber())
        {
          isotopeserie = it->getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber();
        }
      else
        {
          it = _peak_ion_match_list.erase(it);
        }
    }
  qDebug();
}
} // namespace pappso