peptidenaturalisotopelist.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/peptide/peptidenaturalisotopelist.cpp
 * \date 8/3/2015
 * \author Olivier Langella
 * \brief peptide natural isotope model
 */
 
/*******************************************************************************
 * 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
 ******************************************************************************/
// make test ARGS="-V -I 7,7"
#include <QDebug>
#include "peptidenaturalisotopelist.h"
#include "pappsomspp/core/pappsoexception.h"
 
namespace pappso
{
 
PeptideNaturalIsotopeList::PeptideNaturalIsotopeList(const PeptideInterfaceSp &peptide,
                                                     pappso_double minimum_ratio_to_compute)
  : msp_peptide(peptide)
{
  // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__ << " "
  //          << minimum_ratio_to_compute;
  int number_of_fixed_oxygen = msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::O18) +
                               msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::O17);
  int number_of_fixed_sulfur = msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::S33) +
                               msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::S34) +
                               msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::S36);
  int number_of_fixed_nitrogen = msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::N15);
  int number_of_fixed_hydrogen = msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::H2);
 
  int total_carbon(msp_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::C) -
                   msp_peptide.get()->getNumberOfIsotope(Enums::Isotope::C13));
 
  // qDebug() << "total_carbon " << total_carbon;
  // qDebug() << "total_sulfur " << total_sulfur;
  std::map<Enums::Isotope, int> map_isotope;
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::C13, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::H2, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::N15, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::O17, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::O18, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::S33, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::S34, 0));
  map_isotope.insert(std::pair<Enums::Isotope, int>(Enums::Isotope::S36, 0));
 
  for(int nbc13 = 0; nbc13 <= total_carbon; nbc13++)
    {
      map_isotope[Enums::Isotope::C13] = nbc13;
      PeptideNaturalIsotopeSp pepIsotope =
        std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
      this->msp_peptide_natural_isotope_list.push_back(pepIsotope);
      if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
        {
          break;
        }
    }
  std::list<PeptideNaturalIsotopeSp> temp_list;
 
  // ******************************************************************
  // Sulfur isotope list
  int total_sulfur(msp_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::S) -
                   number_of_fixed_sulfur);
  std::list<PeptideNaturalIsotopeSp>::iterator it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      map_isotope = it->get()->getIsotopeMap();
      for(int nbS34 = 1; nbS34 <= total_sulfur; nbS34++)
        {
          map_isotope[Enums::Isotope::S34] = nbS34;
          PeptideNaturalIsotopeSp pepIsotope =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotope);
          if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // qDebug() << "peptide " << pepIsotope.get()->getFormula(1) << "
              // " << pepIsotope.get()->getIntensityRatio(1); it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
 
  // compute S33 abundance
  temp_list.resize(0);
  it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      // qDebug() << "peptide S33 " << it->get()->getFormula(1) << " "
      // <<it->get()->getIntensityRatio(1);
      map_isotope = it->get()->getIsotopeMap();
      for(int nbS33 = 1; nbS33 <= (total_sulfur - map_isotope[Enums::Isotope::S34]); nbS33++)
        {
          map_isotope[Enums::Isotope::S33] = nbS33;
          PeptideNaturalIsotopeSp pepIsotopeS33 =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotopeS33);
          if(pepIsotopeS33.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
 
  // compute S36 abundance
  temp_list.resize(0);
  it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      map_isotope = it->get()->getIsotopeMap();
      for(int nbS36 = 1; nbS36 <= (total_sulfur - map_isotope[Enums::Isotope::S34] -
                                   map_isotope[Enums::Isotope::S33]);
          nbS36++)
        {
          map_isotope[Enums::Isotope::S36] = nbS36;
          PeptideNaturalIsotopeSp pepIsotopeS36 =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotopeS36);
          if(pepIsotopeS36.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
 
  // ******************************************************************
 
 
  // ******************************************************************
  // Hydrogen isotope list
  temp_list.resize(0);
  // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
  // total_hydrogen";
  int total_hydrogen(msp_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::H) -
                     number_of_fixed_hydrogen);
  it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
      // getIsotopeMap " << it->getFormula(1) << " " <<
      // msp_peptide_natural_isotope_list.size();
      map_isotope = it->get()->getIsotopeMap();
      for(int nbH2 = 1; nbH2 <= total_hydrogen; nbH2++)
        {
          map_isotope[Enums::Isotope::H2] = nbH2;
          PeptideNaturalIsotopeSp pepIsotope =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotope);
          if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
  // ******************************************************************
 
 
  // ******************************************************************
  // Oxygen isotope list
  temp_list.resize(0);
  // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
  // total_oxygen";
  unsigned int total_oxygen(msp_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::O) -
                            number_of_fixed_oxygen);
  it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
      // getIsotopeMap " << it->getFormula(1) << " " <<
      // msp_peptide_natural_isotope_list.size();
      map_isotope = it->get()->getIsotopeMap();
      for(unsigned int nbO18 = 1; nbO18 <= total_oxygen; nbO18++)
        {
          // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
          // nbO18 " << nbO18;
          map_isotope[Enums::Isotope::O18] = nbO18;
          PeptideNaturalIsotopeSp pepIsotope =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotope);
          if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
  // ******************************************************************
 
 
  // ******************************************************************
  // Nitrogen isotope list
  temp_list.resize(0);
  // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
  // total_nitrogen";
  unsigned int total_nitrogen(msp_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::N) -
                              number_of_fixed_nitrogen);
  it = msp_peptide_natural_isotope_list.begin();
  while(it != msp_peptide_natural_isotope_list.end())
    {
      // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
      // getIsotopeMap " << it->getFormula(1) << " " <<
      // msp_peptide_natural_isotope_list.size();
      map_isotope = it->get()->getIsotopeMap();
      for(unsigned int nbN15 = 1; nbN15 <= total_nitrogen; nbN15++)
        {
          // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
          // nbN15 " << nbN15;
          map_isotope[Enums::Isotope::N15] = nbN15;
          PeptideNaturalIsotopeSp pepIsotope =
            std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
          temp_list.push_back(pepIsotope);
          if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
            {
              // it++;
              break;
            }
        }
      it++;
    }
  msp_peptide_natural_isotope_list.insert(it, temp_list.begin(), temp_list.end());
  // ******************************************************************
 
  // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList end
  // size="<<msp_peptide_natural_isotope_list.size();
}
 
PeptideNaturalIsotopeListSp
PeptideNaturalIsotopeList::makePeptideNaturalIsotopeListSp() const
{
  return std::make_shared<PeptideNaturalIsotopeList>(*this);
}
 
PeptideNaturalIsotopeList::PeptideNaturalIsotopeList(const PeptideNaturalIsotopeList &other)
  : msp_peptide(other.msp_peptide),
    msp_peptide_natural_isotope_list(other.msp_peptide_natural_isotope_list)
{
}
 
PeptideNaturalIsotopeList::~PeptideNaturalIsotopeList()
{
}
 
const std::map<unsigned int, pappso_double>
PeptideNaturalIsotopeList::getIntensityRatioPerIsotopeNumber() const
{
  std::list<PeptideNaturalIsotopeSp>::const_iterator it = msp_peptide_natural_isotope_list.begin();
  std::map<unsigned int, pappso_double> map_isotope_number;
 
  while(it != msp_peptide_natural_isotope_list.end())
    {
      unsigned int number = it->get()->getIsotopeNumber();
      std::pair<std::map<unsigned int, pappso_double>::iterator, bool> mapnew =
        map_isotope_number.insert(std::pair<unsigned int, pappso_double>(number, 0));
      if(mapnew.second == false)
        {
          // mapit = map_isotope_number.insert(std::pair<unsigned
          // int,pappso_double>(number, 0));
        }
      mapnew.first->second += it->get()->getIntensityRatio(1);
      it++;
    }
  return map_isotope_number;
}
 
 
/** /brief get a sorted (by expected intensity) vector of isotopes of the same
 * level
 *
 * */
 
std::vector<PeptideNaturalIsotopeSp>
PeptideNaturalIsotopeList::getByIsotopeNumber(unsigned int isotope_number,
                                              unsigned int charge) const
{
  std::vector<PeptideNaturalIsotopeSp> v_isotope_list;
 
  for(auto &&isotopeSp : msp_peptide_natural_isotope_list)
    {
      if(isotopeSp.get()->getIsotopeNumber() == isotope_number)
        {
          v_isotope_list.push_back(isotopeSp);
        }
    }
  std::sort(v_isotope_list.begin(),
            v_isotope_list.end(),
            [charge](const PeptideNaturalIsotopeSp &m, const PeptideNaturalIsotopeSp &n) {
              return (m.get()->getIntensityRatio(charge) > n.get()->getIntensityRatio(charge));
            });
  return v_isotope_list;
}
 
 
/** /brief get a sorted (by expected intensity) vector of natural isotope
 * average by isotope number
 *
 * */
std::vector<PeptideNaturalIsotopeAverageSp>
getByIntensityRatioByIsotopeNumber(const PeptideInterfaceSp &peptide,
                                   unsigned int charge,
                                   PrecisionPtr precision,
                                   unsigned int isotopeNumber,
                                   pappso_double minimumIntensity)
{
 
  // qDebug() << "getByIntensityRatioByIsotopeNumber begin";
  unsigned int askedIsotopeRank;
  unsigned int maxAskedIsotopeRank = 10;
  pappso_double cumulativeRatio    = 0;
  std::vector<PeptideNaturalIsotopeAverageSp> v_isotopeAverageList;
  std::vector<PeptideNaturalIsotopeAverageSp> v_isotopeAverageListResult;
 
  std::vector<unsigned int> previousIsotopeRank;
  bool isEmpty = false;
  for(askedIsotopeRank = 1; (askedIsotopeRank < maxAskedIsotopeRank) && (!isEmpty);
      askedIsotopeRank++)
    {
      PeptideNaturalIsotopeAverage isotopeAverage(
        peptide, askedIsotopeRank, isotopeNumber, charge, precision);
      isEmpty = isotopeAverage.isEmpty();
      if(isEmpty)
        {
        }
      else
        {
          if(std::find(previousIsotopeRank.begin(),
                       previousIsotopeRank.end(),
                       isotopeAverage.getIsotopeRank()) == previousIsotopeRank.end())
            { // not Found
              previousIsotopeRank.push_back(isotopeAverage.getIsotopeRank());
              v_isotopeAverageList.push_back(isotopeAverage.makePeptideNaturalIsotopeAverageSp());
            }
        }
    }
  if(v_isotopeAverageList.size() == 0)
    return v_isotopeAverageListResult;
 
  // qDebug() << "getByIntensityRatioByIsotopeNumber comp";
  std::sort(v_isotopeAverageList.begin(),
            v_isotopeAverageList.end(),
            [](const PeptideNaturalIsotopeAverageSp &m, const PeptideNaturalIsotopeAverageSp &n) {
              return (m.get()->getIntensityRatio() > n.get()->getIntensityRatio());
            });
 
  cumulativeRatio = 0;
  auto it         = v_isotopeAverageList.begin();
  v_isotopeAverageListResult.clear();
  // qDebug() << "getByIntensityRatioByIsotopeNumber cumul";
  while((it != v_isotopeAverageList.end()) && (cumulativeRatio < minimumIntensity))
    {
      cumulativeRatio += it->get()->getIntensityRatio();
      v_isotopeAverageListResult.push_back(*it);
      it++;
    }
 
  // qDebug() << "getByIntensityRatioByIsotopeNumber end";
 
  return v_isotopeAverageListResult;
}
 
 
/** /brief get a sorted (by expected intensity) vector of natural isotope
 * average
 *
 * */
std::vector<PeptideNaturalIsotopeAverageSp>
PeptideNaturalIsotopeList::getByIntensityRatio(unsigned int charge,
                                               PrecisionPtr precision,
                                               pappso_double minimumIntensityRatio) const
{
 
  // qDebug() << "PeptideNaturalIsotopeList::getByIntensityRatio begin";
  std::vector<PeptideNaturalIsotopeAverageSp> peptide_natural_isotope_average_list;
 
  std::map<unsigned int, pappso::pappso_double> map_isotope_number =
    getIntensityRatioPerIsotopeNumber();
  std::vector<std::pair<unsigned int, pappso::pappso_double>> sorted_number_ratio;
 
  for(unsigned int i = 0; i < map_isotope_number.size(); i++)
    {
      sorted_number_ratio.push_back(
        std::pair<unsigned int, pappso::pappso_double>(i, map_isotope_number[i]));
      unsigned int asked_rank = 0;
      unsigned int given_rank = 0;
      bool more_rank          = true;
      while(more_rank)
        {
          asked_rank++;
          pappso::PeptideNaturalIsotopeAverage isotopeAverageMono(
            *this, asked_rank, i, charge, precision);
          given_rank = isotopeAverageMono.getIsotopeRank();
          if(given_rank < asked_rank)
            {
              more_rank = false;
            }
          else if(isotopeAverageMono.getIntensityRatio() == 0)
            {
              more_rank = false;
            }
          else
            {
              // isotopeAverageMono.makePeptideNaturalIsotopeAverageSp();
              peptide_natural_isotope_average_list.push_back(
                isotopeAverageMono.makePeptideNaturalIsotopeAverageSp());
            }
        }
    }
 
 
  // sort by intensity ratio
  std::sort(
    sorted_number_ratio.begin(),
    sorted_number_ratio.end(),
    [](const std::pair<unsigned int, pappso::pappso_double> &m,
       const std::pair<unsigned int, pappso::pappso_double> &n) { return (m.second > n.second); });
 
 
  double cumulativeRatio = 0;
  std::vector<unsigned int> selected_isotope_number_list;
  for(auto &pair_isotope_number : sorted_number_ratio)
    {
      if(cumulativeRatio <= minimumIntensityRatio)
        {
          selected_isotope_number_list.push_back(pair_isotope_number.first);
        }
      else
        {
          break;
        }
      cumulativeRatio += pair_isotope_number.second;
    }
 
  auto it_remove =
    std::remove_if(peptide_natural_isotope_average_list.begin(),
                   peptide_natural_isotope_average_list.end(),
                   [selected_isotope_number_list](const PeptideNaturalIsotopeAverageSp &average) {
                     auto it = std::find(selected_isotope_number_list.begin(),
                                         selected_isotope_number_list.end(),
                                         average.get()->getIsotopeNumber());
                     return (it == selected_isotope_number_list.end());
                   });
  peptide_natural_isotope_average_list.erase(it_remove, peptide_natural_isotope_average_list.end());
  return peptide_natural_isotope_average_list;
}
 
 
PeptideNaturalIsotopeList::const_iterator
PeptideNaturalIsotopeList::begin() const
{
  return msp_peptide_natural_isotope_list.begin();
}
 
PeptideNaturalIsotopeList::const_iterator
PeptideNaturalIsotopeList::end() const
{
  return msp_peptide_natural_isotope_list.end();
}
 
unsigned int
PeptideNaturalIsotopeList::size() const
{
  return msp_peptide_natural_isotope_list.size();
}
 
const PeptideInterfaceSp &
PeptideNaturalIsotopeList::getPeptideInterfaceSp() const
{
  return msp_peptide;
}
} // namespace pappso