Loading...
Searching...
No Matches
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 "pappsomspp/core/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(Enums::Isotope isotope, unsigned int total, unsigned int heavy)
{
pappso_double abundance = 1;
switch(isotope)
{
case Enums::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 Enums::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 Enums::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 Enums::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 Enums::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 Enums::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 Enums::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 Enums::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<Enums::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(Enums::Isotope::O18) +
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::O17);
int number_of_fixed_sulfur = m_peptide.get()->getNumberOfIsotope(Enums::Isotope::S33) +
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::S34) +
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::S36);
m_ratio = isotopem_ratio_cache(Enums::Isotope::C13,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::C) -
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::C13),
m_mapIsotope.at(Enums::Isotope::C13));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::N15,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::N) -
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::N15),
m_mapIsotope.at(Enums::Isotope::N15));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::O18,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::O) -
number_of_fixed_oxygen,
m_mapIsotope.at(Enums::Isotope::O18));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::O17,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::O) -
number_of_fixed_oxygen,
m_mapIsotope.at(Enums::Isotope::O17));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::S33,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::S) -
number_of_fixed_sulfur,
m_mapIsotope.at(Enums::Isotope::S33));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::S34,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::S) -
number_of_fixed_sulfur,
m_mapIsotope.at(Enums::Isotope::S34));
m_ratio *= isotopem_ratio_cache(Enums::Isotope::S36,
m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::S) -
number_of_fixed_sulfur,
m_mapIsotope.at(Enums::Isotope::S36));
// qDebug() << "Aa::getMass() begin";
m_mass = m_peptide.get()->getMass();
m_mass += (DIFFC12C13 * m_mapIsotope.at(Enums::Isotope::C13));
m_mass += (DIFFS32S33 * m_mapIsotope.at(Enums::Isotope::S33));
m_mass += (DIFFS32S34 * m_mapIsotope.at(Enums::Isotope::S34));
m_mass += (DIFFS32S36 * m_mapIsotope.at(Enums::Isotope::S36));
m_mass += (DIFFH1H2 * m_mapIsotope.at(Enums::Isotope::H2));
m_mass += (DIFFO16O18 * m_mapIsotope.at(Enums::Isotope::O18));
m_mass += (DIFFN14N15 * m_mapIsotope.at(Enums::Isotope::N15));
m_mass += (DIFFO16O17 * m_mapIsotope.at(Enums::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(
Enums::Isotope::H2,
(m_peptide.get()->getNumberOfAtom(Enums::AtomIsotopeSurvey::H) + charge) -
m_peptide.get()->getNumberOfIsotope(Enums::Isotope::H2),
m_mapIsotope.at(Enums::Isotope::H2));
}
int
PeptideNaturalIsotope::getNumberOfAtom(Enums::AtomIsotopeSurvey atom) const
{
return m_peptide.get()->getNumberOfAtom(atom);
}
int
PeptideNaturalIsotope::getNumberOfIsotope(Enums::Isotope isotope) const
{
return m_mapIsotope.at(isotope) + m_peptide.get()->getNumberOfIsotope(isotope);
}
const std::map<Enums::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(Enums::Isotope::C13) + m_mapIsotope.at(Enums::Isotope::H2) +
m_mapIsotope.at(Enums::Isotope::O17) + m_mapIsotope.at(Enums::Isotope::O18) * 2 +
m_mapIsotope.at(Enums::Isotope::N15) + m_mapIsotope.at(Enums::Isotope::S33) +
(m_mapIsotope.at(Enums::Isotope::S34) * 2) + (m_mapIsotope.at(Enums::Isotope::S36) * 4);
}
} // namespace pappso
pappso::PeptideNaturalIsotope::m_mass
pappso_double m_mass
Definition peptidenaturalisotope.h:73
pappso::PeptideNaturalIsotope::getSequence
virtual const QString getSequence() const override
amino acid sequence without modification
Definition peptidenaturalisotope.cpp:327
pappso::PeptideNaturalIsotope::getIsotopeMap
const std::map< Enums::Isotope, int > & getIsotopeMap() const
Definition peptidenaturalisotope.cpp:307
pappso::PeptideNaturalIsotope::getNumberOfIsotope
virtual int getNumberOfIsotope(Enums::Isotope isotope) const override
get the number of isotopes C13, H2, O17, O18, N15, S33, S34, S36 in the molecule
Definition peptidenaturalisotope.cpp:301
pappso::PeptideNaturalIsotope::m_peptide
const PeptideInterfaceSp m_peptide
Definition peptidenaturalisotope.h:69
pappso::PeptideNaturalIsotope::m_ratio
pappso_double m_ratio
Definition peptidenaturalisotope.h:72
pappso::PeptideNaturalIsotope::PeptideNaturalIsotope
PeptideNaturalIsotope(const PeptideInterfaceSp &peptide, const std::map< Enums::Isotope, int > &map_isotope)
Definition peptidenaturalisotope.cpp:195
pappso::PeptideNaturalIsotope::isPalindrome
virtual bool isPalindrome() const override
tells if the peptide sequence is a palindrome
Definition peptidenaturalisotope.cpp:314
pappso::PeptideNaturalIsotope::size
virtual unsigned int size() const override
Definition peptidenaturalisotope.cpp:321
pappso::PeptideNaturalIsotope::getIsotopeNumber
virtual unsigned int getIsotopeNumber() const
Definition peptidenaturalisotope.cpp:333
pappso::PeptideNaturalIsotope::~PeptideNaturalIsotope
virtual ~PeptideNaturalIsotope()
Definition peptidenaturalisotope.cpp:269
pappso::PeptideNaturalIsotope::getIntensityRatio
pappso_double getIntensityRatio(unsigned int charge) const
Definition peptidenaturalisotope.cpp:283
pappso::PeptideNaturalIsotope::getNumberOfAtom
virtual int getNumberOfAtom(Enums::AtomIsotopeSurvey atom) const override
get the number of atom C, O, N, H in the molecule
Definition peptidenaturalisotope.cpp:295
pappso::PeptideNaturalIsotope::m_mapIsotope
const std::map< Enums::Isotope, int > m_mapIsotope
Definition peptidenaturalisotope.h:71
pappso::PeptideNaturalIsotope::getMass
pappso_double getMass() const override
Definition peptidenaturalisotope.cpp:275
pappso::Enums::AtomIsotopeSurvey
AtomIsotopeSurvey
Definition types.h:110
pappso::Enums::AtomIsotopeSurvey::C
@ C
Definition types.h:111
pappso::Enums::AtomIsotopeSurvey::S
@ S
Definition types.h:115
pappso::Enums::AtomIsotopeSurvey::N
@ N
Definition types.h:114
pappso::Enums::AtomIsotopeSurvey::H
@ H
Definition types.h:112
pappso::Enums::AtomIsotopeSurvey::O
@ O
Definition types.h:113
pappso::Enums::Isotope
Isotope
Definition types.h:126
pappso::Enums::Isotope::S36
@ S36
Definition types.h:139
pappso::Enums::Isotope::C13
@ C13
Definition types.h:128
pappso::Enums::Isotope::N15
@ N15
Definition types.h:135
pappso::Enums::Isotope::S33
@ S33
Definition types.h:137
pappso::Enums::Isotope::O18
@ O18
Definition types.h:133
pappso::Enums::Isotope::O17
@ O17
Definition types.h:132
pappso::Enums::Isotope::S34
@ S34
Definition types.h:138
pappso::Enums::Isotope::H2
@ H2
Definition types.h:130
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappso::isotopem_ratio_cache
pappso_double isotopem_ratio_cache(Enums::Isotope isotope, unsigned int total, unsigned int heavy)
Definition peptidenaturalisotope.cpp:95
pappso::DIFFS32S33
const pappso_double DIFFS32S33(32.9714589101 - MASSSULFUR)
pappso::Combinations
uint64_t Combinations(unsigned int n, unsigned int k)
Definition peptidenaturalisotope.cpp:47
pappso::DIFFS32S34
const pappso_double DIFFS32S34(33.9678670300 - MASSSULFUR)
pappso::ratioO17_cache
pappso_double ratioO17_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:90
pappso::DIFFO16O17
const pappso_double DIFFO16O17(16.99913150 - MASSOXYGEN)
pappso::ABUNDANCES36
const pappso_double ABUNDANCES36(0.00010999120070394368536836893213148869108408689498901367187500)
pappso::PeptideInterfaceSp
std::shared_ptr< const PeptideInterface > PeptideInterfaceSp
Definition peptideinterface.h:40
pappso::ratioS34_cache
pappso_double ratioS34_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:88
pappso::ABUNDANCEN15
const pappso_double ABUNDANCEN15(0.00364198543205827118818262988497735932469367980957031250000000)
pappso::DIFFS32S36
const pappso_double DIFFS32S36(35.9670812000 - MASSSULFUR)
pappso::ratioC13_cache
pappso_double ratioC13_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:85
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition types.h:60
pappso::isotopem_ratio
pappso_double isotopem_ratio(pappso_double abundance, unsigned int total, unsigned int heavy)
Definition peptidenaturalisotope.cpp:79
pappso::ratioH2_cache
pappso_double ratioH2_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:92
pappso::ABUNDANCEC13
const pappso_double ABUNDANCEC13(0.01078805814953308406245469086570665240287780761718750000000000)
pappso::ABUNDANCEO17
const pappso_double ABUNDANCEO17(0.00038099847600609595965615028489992255344986915588378906250000)
pappso::ABUNDANCEH2
const pappso_double ABUNDANCEH2(0.00011570983569203332000374651045149221317842602729797363281250)
pappso::ratioO18_cache
pappso_double ratioO18_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:91
pappso::ABUNDANCES34
const pappso_double ABUNDANCES34(0.04252059835213182203972337447339668869972229003906250000000000)
pappso::ratioS36_cache
pappso_double ratioS36_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:87
pappso::ABUNDANCEO18
const pappso_double ABUNDANCEO18(0.00205139179443282221315669744399201590567827224731445312500000)
pappso::AtomIsotope
AtomIsotope
Definition peptidenaturalisotope.cpp:70
pappso::AtomIsotope::C
@ C
Definition peptidenaturalisotope.cpp:71
pappso::AtomIsotope::S
@ S
Definition peptidenaturalisotope.cpp:75
pappso::AtomIsotope::N
@ N
Definition peptidenaturalisotope.cpp:74
pappso::AtomIsotope::H
@ H
Definition peptidenaturalisotope.cpp:72
pappso::AtomIsotope::O
@ O
Definition peptidenaturalisotope.cpp:73
pappso::combinations_cache
uint64_t combinations_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:44
pappso::DIFFO16O18
const pappso_double DIFFO16O18(17.9991610 - MASSOXYGEN)
pappso::ABUNDANCES33
const pappso_double ABUNDANCES33(0.00751939844812414937003097747947322204709053039550781250000000)
pappso::DIFFN14N15
const pappso_double DIFFN14N15(15.0001088982 - MASSNITROGEN)
pappso::ratioN15_cache
pappso_double ratioN15_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:86
pappso::ratioS33_cache
pappso_double ratioS33_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition peptidenaturalisotope.cpp:89
pappso::DIFFC12C13
const pappso_double DIFFC12C13(1.0033548378)
pappso::DIFFH1H2
const pappso_double DIFFH1H2(2.0141017778 - MPROTIUM)
pappsoexception.h
CACHE_ARRAY_SIZE
#define CACHE_ARRAY_SIZE
Definition peptidenaturalisotope.cpp:42
peptidenaturalisotope.h