peptidenaturalisotope.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/peptide/peptidenaturalisotope.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
 ******************************************************************************/
 
#include "peptidenaturalisotope.h"
#include "../pappsoexception.h"
 
#include <cmath>
#include <QDebug>
 
using namespace std;
 
namespace pappso
{
 
#define CACHE_ARRAY_SIZE 500
 
uint64_t combinations_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
 
uint64_t
Combinations(unsigned int n, unsigned int k)
{
  if(k > n)
    return 0;
 
  uint64_t r = 1;
  if((n < CACHE_ARRAY_SIZE) && (combinations_cache[n][k] != 0))
    {
      return combinations_cache[n][k];
    }
  for(unsigned int d = 1; d <= k; ++d)
    {
      r *= n--;
      r /= d;
    }
  if(n < CACHE_ARRAY_SIZE)
    {
      combinations_cache[n][k] = r;
    }
  return r;
}
 
enum class AtomIsotope
{
  C,
  H,
  O,
  N,
  S
};
 
pappso_double
isotopem_ratio(pappso_double abundance, unsigned int total, unsigned int heavy)
{
  return (pow(abundance, heavy) * pow((double)1 - abundance, (total - heavy)) *
          (double)Combinations(total, heavy));
}
 
pappso_double ratioC13_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioN15_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioS36_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioS34_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioS33_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioO17_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioO18_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE] = {0};
pappso_double ratioH2_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]  = {0};
 
pappso_double
isotopem_ratio_cache(Isotope isotope, unsigned int total, unsigned int heavy)
{
  pappso_double abundance = 1;
  switch(isotope)
    {
      case Isotope::H2:
        abundance = ABUNDANCEH2;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioH2_cache[total][heavy] == 0)
              {
                ratioH2_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioH2_cache[total][heavy];
          }
        break;
      case Isotope::C13:
        abundance = ABUNDANCEC13;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioC13_cache[total][heavy] == 0)
              {
                ratioC13_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioC13_cache[total][heavy];
          }
        break;
      case Isotope::N15:
        abundance = ABUNDANCEN15;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioN15_cache[total][heavy] == 0)
              {
                ratioN15_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioN15_cache[total][heavy];
          }
        break;
      case Isotope::O18:
        abundance = ABUNDANCEO18;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioO18_cache[total][heavy] == 0)
              {
                ratioO18_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioO18_cache[total][heavy];
          }
        break;
      case Isotope::O17:
        abundance = ABUNDANCEO17;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioO17_cache[total][heavy] == 0)
              {
                ratioO17_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioO17_cache[total][heavy];
          }
        break;
      case Isotope::S33:
        abundance = ABUNDANCES33;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioS33_cache[total][heavy] == 0)
              {
                ratioS33_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioS33_cache[total][heavy];
          }
        break;
      case Isotope::S34:
        abundance = ABUNDANCES34;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioS34_cache[total][heavy] == 0)
              {
                ratioS34_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioS34_cache[total][heavy];
          }
        break;
      case Isotope::S36:
        abundance = ABUNDANCES36;
        if(total < CACHE_ARRAY_SIZE)
          {
            if(ratioS36_cache[total][heavy] == 0)
              {
                ratioS36_cache[total][heavy] = isotopem_ratio(abundance, total, heavy);
              }
            return ratioS36_cache[total][heavy];
          }
        break;
      default:
        break;
    }
  return isotopem_ratio(abundance, total, heavy);
}
 
 
PeptideNaturalIsotope::PeptideNaturalIsotope(const PeptideInterfaceSp &peptide,
                                             const std::map<Isotope, int> &map_isotope)
  : m_peptide(peptide), m_mapIsotope(map_isotope)
{
  //_abundance = ((_number_of_carbon - number_of_C13) * ABUNDANCEC12) +
  //(number_of_C13 * ABUNDANCEC13); p = pow(0.01, i)*pow(0.99, (c-i))*comb(c,i)
  // qDebug()<< "pow" << pow(ABUNDANCEC13, number_of_C13)*pow(1-ABUNDANCEC13,
  // (_number_of_carbon-number_of_C13));
  // qDebug() <<"conb" << Combinations(_number_of_carbon,number_of_C13);
 
  // CHNO
  //_probC13 = pow(ABUNDANCEC13, number_of_C13)*pow((double)1-ABUNDANCEC13,
  //(_number_of_carbon-number_of_C13))* (double)
  // Combinations(_number_of_carbon,number_of_C13);
  // qDebug() <<"_probC13" <<_probC13;
 
  // number of fixed Oxygen atoms (already labelled, not natural) :
  int number_of_fixed_oxygen = m_peptide.get()->getNumberOfIsotope(Isotope::O18) +
                               m_peptide.get()->getNumberOfIsotope(Isotope::O17);
  int number_of_fixed_sulfur = m_peptide.get()->getNumberOfIsotope(Isotope::S33) +
                               m_peptide.get()->getNumberOfIsotope(Isotope::S34) +
                               m_peptide.get()->getNumberOfIsotope(Isotope::S36);
 
  m_ratio = isotopem_ratio_cache(Isotope::C13,
                                 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::C) -
                                   m_peptide.get()->getNumberOfIsotope(Isotope::C13),
                                 m_mapIsotope.at(Isotope::C13));
  m_ratio *= isotopem_ratio_cache(Isotope::N15,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::N) -
                                    m_peptide.get()->getNumberOfIsotope(Isotope::N15),
                                  m_mapIsotope.at(Isotope::N15));
  m_ratio *= isotopem_ratio_cache(Isotope::O18,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
                                    number_of_fixed_oxygen,
                                  m_mapIsotope.at(Isotope::O18));
  m_ratio *= isotopem_ratio_cache(Isotope::O17,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
                                    number_of_fixed_oxygen,
                                  m_mapIsotope.at(Isotope::O17));
  m_ratio *= isotopem_ratio_cache(Isotope::S33,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
                                    number_of_fixed_sulfur,
                                  m_mapIsotope.at(Isotope::S33));
  m_ratio *= isotopem_ratio_cache(Isotope::S34,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
                                    number_of_fixed_sulfur,
                                  m_mapIsotope.at(Isotope::S34));
  m_ratio *= isotopem_ratio_cache(Isotope::S36,
                                  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
                                    number_of_fixed_sulfur,
                                  m_mapIsotope.at(Isotope::S36));
 
 
  // qDebug() << "Aa::getMass() begin";
  m_mass = m_peptide.get()->getMass();
  m_mass += (DIFFC12C13 * m_mapIsotope.at(Isotope::C13));
  m_mass += (DIFFS32S33 * m_mapIsotope.at(Isotope::S33));
  m_mass += (DIFFS32S34 * m_mapIsotope.at(Isotope::S34));
  m_mass += (DIFFS32S36 * m_mapIsotope.at(Isotope::S36));
  m_mass += (DIFFH1H2 * m_mapIsotope.at(Isotope::H2));
  m_mass += (DIFFO16O18 * m_mapIsotope.at(Isotope::O18));
  m_mass += (DIFFN14N15 * m_mapIsotope.at(Isotope::N15));
  m_mass += (DIFFO16O17 * m_mapIsotope.at(Isotope::O17));
 
 
  // qDebug() << "Aa::getMass() end " << mass;
}
 
PeptideNaturalIsotope::PeptideNaturalIsotope(const PeptideNaturalIsotope &other)
  : m_peptide(other.m_peptide), m_mapIsotope(other.m_mapIsotope)
{
  m_ratio = other.m_ratio;
}
 
PeptideNaturalIsotope::~PeptideNaturalIsotope()
{
}
 
 
pappso_double
PeptideNaturalIsotope::getMass() const
{
 
  return m_mass;
}
 
 
pappso_double
PeptideNaturalIsotope::getIntensityRatio(unsigned int charge) const
{
 
  return m_ratio *
         isotopem_ratio_cache(Isotope::H2,
                              (m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::H) + charge) -
                                m_peptide.get()->getNumberOfIsotope(Isotope::H2),
                              m_mapIsotope.at(Isotope::H2));
}
 
 
int
PeptideNaturalIsotope::getNumberOfAtom(AtomIsotopeSurvey atom) const
{
  return m_peptide.get()->getNumberOfAtom(atom);
}
 
int
PeptideNaturalIsotope::getNumberOfIsotope(Isotope isotope) const
{
  return m_mapIsotope.at(isotope) + m_peptide.get()->getNumberOfIsotope(isotope);
}
 
const std::map<Isotope, int> &
PeptideNaturalIsotope::getIsotopeMap() const
{
  return m_mapIsotope;
}
 
 
bool
PeptideNaturalIsotope::isPalindrome() const
{
  return m_peptide.get()->isPalindrome();
}
 
 
unsigned int
PeptideNaturalIsotope::size() const
{
  return m_peptide.get()->size();
}
 
const QString
PeptideNaturalIsotope::getSequence() const
{
  return m_peptide.get()->getSequence();
}
 
unsigned int
PeptideNaturalIsotope::getIsotopeNumber() const
{
  // only count variable (natural) isotope
  return m_mapIsotope.at(Isotope::C13) + m_mapIsotope.at(Isotope::H2) +
         m_mapIsotope.at(Isotope::O17) + m_mapIsotope.at(Isotope::O18) * 2 +
         m_mapIsotope.at(Isotope::N15) + m_mapIsotope.at(Isotope::S33) +
         (m_mapIsotope.at(Isotope::S34) * 2) + (m_mapIsotope.at(Isotope::S36) * 4);
}
} // namespace pappso