libpappsomspp
Library for mass spectrometry
peptidenaturalisotope.cpp
Go to the documentation of this file.
1 /**
2  * \file pappsomspp/peptide/peptidenaturalisotope.cpp
3  * \date 8/3/2015
4  * \author Olivier Langella
5  * \brief peptide natural isotope model
6  */
7 
8 /*******************************************************************************
9  * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.fr>.
10  *
11  * This file is part of the PAPPSOms++ library.
12  *
13  * PAPPSOms++ is free software: you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation, either version 3 of the License, or
16  * (at your option) any later version.
17  *
18  * PAPPSOms++ is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
25  *
26  * Contributors:
27  * Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
28  *implementation
29  ******************************************************************************/
30 
31 #include "peptidenaturalisotope.h"
32 #include "../pappsoexception.h"
33 
34 #include <cmath>
35 #include <QDebug>
36 
37 using namespace std;
38 
39 namespace pappso
40 {
41 
42 #define CACHE_ARRAY_SIZE 500
43 
45 
46 uint64_t
47 Combinations(unsigned int n, unsigned int k)
48 {
49  if(k > n)
50  return 0;
51 
52  uint64_t r = 1;
53  if((n < CACHE_ARRAY_SIZE) && (combinations_cache[n][k] != 0))
54  {
55  return combinations_cache[n][k];
56  }
57  for(unsigned int d = 1; d <= k; ++d)
58  {
59  r *= n--;
60  r /= d;
61  }
62  if(n < CACHE_ARRAY_SIZE)
63  {
64  combinations_cache[n][k] = r;
65  }
66  return r;
67 }
68 
69 enum class AtomIsotope
70 {
71  C,
72  H,
73  O,
74  N,
75  S
76 };
77 
79 isotopem_ratio(pappso_double abundance, unsigned int total, unsigned int heavy)
80 {
81  return (pow(abundance, heavy) * pow((double)1 - abundance, (total - heavy)) *
82  (double)Combinations(total, heavy));
83 }
84 
93 
95 isotopem_ratio_cache(Isotope isotope, unsigned int total, unsigned int heavy)
96 {
97  pappso_double abundance = 1;
98  switch(isotope)
99  {
100  case Isotope::H2:
101  abundance = ABUNDANCEH2;
102  if(total < CACHE_ARRAY_SIZE)
103  {
104  if(ratioH2_cache[total][heavy] == 0)
105  {
106  ratioH2_cache[total][heavy] =
107  isotopem_ratio(abundance, total, heavy);
108  }
109  return ratioH2_cache[total][heavy];
110  }
111  break;
112  case Isotope::C13:
113  abundance = ABUNDANCEC13;
114  if(total < CACHE_ARRAY_SIZE)
115  {
116  if(ratioC13_cache[total][heavy] == 0)
117  {
118  ratioC13_cache[total][heavy] =
119  isotopem_ratio(abundance, total, heavy);
120  }
121  return ratioC13_cache[total][heavy];
122  }
123  break;
124  case Isotope::N15:
125  abundance = ABUNDANCEN15;
126  if(total < CACHE_ARRAY_SIZE)
127  {
128  if(ratioN15_cache[total][heavy] == 0)
129  {
130  ratioN15_cache[total][heavy] =
131  isotopem_ratio(abundance, total, heavy);
132  }
133  return ratioN15_cache[total][heavy];
134  }
135  break;
136  case Isotope::O18:
137  abundance = ABUNDANCEO18;
138  if(total < CACHE_ARRAY_SIZE)
139  {
140  if(ratioO18_cache[total][heavy] == 0)
141  {
142  ratioO18_cache[total][heavy] =
143  isotopem_ratio(abundance, total, heavy);
144  }
145  return ratioO18_cache[total][heavy];
146  }
147  break;
148  case Isotope::O17:
149  abundance = ABUNDANCEO17;
150  if(total < CACHE_ARRAY_SIZE)
151  {
152  if(ratioO17_cache[total][heavy] == 0)
153  {
154  ratioO17_cache[total][heavy] =
155  isotopem_ratio(abundance, total, heavy);
156  }
157  return ratioO17_cache[total][heavy];
158  }
159  break;
160  case Isotope::S33:
161  abundance = ABUNDANCES33;
162  if(total < CACHE_ARRAY_SIZE)
163  {
164  if(ratioS33_cache[total][heavy] == 0)
165  {
166  ratioS33_cache[total][heavy] =
167  isotopem_ratio(abundance, total, heavy);
168  }
169  return ratioS33_cache[total][heavy];
170  }
171  break;
172  case Isotope::S34:
173  abundance = ABUNDANCES34;
174  if(total < CACHE_ARRAY_SIZE)
175  {
176  if(ratioS34_cache[total][heavy] == 0)
177  {
178  ratioS34_cache[total][heavy] =
179  isotopem_ratio(abundance, total, heavy);
180  }
181  return ratioS34_cache[total][heavy];
182  }
183  break;
184  case Isotope::S36:
185  abundance = ABUNDANCES36;
186  if(total < CACHE_ARRAY_SIZE)
187  {
188  if(ratioS36_cache[total][heavy] == 0)
189  {
190  ratioS36_cache[total][heavy] =
191  isotopem_ratio(abundance, total, heavy);
192  }
193  return ratioS36_cache[total][heavy];
194  }
195  break;
196  }
197  return isotopem_ratio(abundance, total, heavy);
198 }
199 
200 
201 PeptideNaturalIsotope::PeptideNaturalIsotope(
202  const PeptideInterfaceSp &peptide, const std::map<Isotope, int> &map_isotope)
203  : m_peptide(peptide), m_mapIsotope(map_isotope)
204 {
205  //_abundance = ((_number_of_carbon - number_of_C13) * ABUNDANCEC12) +
206  //(number_of_C13 * ABUNDANCEC13); p = pow(0.01, i)*pow(0.99, (c-i))*comb(c,i)
207  // qDebug()<< "pow" << pow(ABUNDANCEC13, number_of_C13)*pow(1-ABUNDANCEC13,
208  // (_number_of_carbon-number_of_C13));
209  // qDebug() <<"conb" << Combinations(_number_of_carbon,number_of_C13);
210 
211  // CHNO
212  //_probC13 = pow(ABUNDANCEC13, number_of_C13)*pow((double)1-ABUNDANCEC13,
213  //(_number_of_carbon-number_of_C13))* (double)
214  // Combinations(_number_of_carbon,number_of_C13);
215  // qDebug() <<"_probC13" <<_probC13;
216 
217  // number of fixed Oxygen atoms (already labelled, not natural) :
218  int number_of_fixed_oxygen =
219  m_peptide.get()->getNumberOfIsotope(Isotope::O18) +
220  m_peptide.get()->getNumberOfIsotope(Isotope::O17);
221  int number_of_fixed_sulfur =
222  m_peptide.get()->getNumberOfIsotope(Isotope::S33) +
223  m_peptide.get()->getNumberOfIsotope(Isotope::S34) +
224  m_peptide.get()->getNumberOfIsotope(Isotope::S36);
225 
227  Isotope::C13,
228  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::C) -
229  m_peptide.get()->getNumberOfIsotope(Isotope::C13),
232  Isotope::N15,
233  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::N) -
234  m_peptide.get()->getNumberOfIsotope(Isotope::N15),
237  Isotope::O18,
238  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
239  number_of_fixed_oxygen,
242  Isotope::O17,
243  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
244  number_of_fixed_oxygen,
247  Isotope::S33,
248  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
249  number_of_fixed_sulfur,
252  Isotope::S34,
253  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
254  number_of_fixed_sulfur,
257  Isotope::S36,
258  m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
259  number_of_fixed_sulfur,
261 
262 
263  // qDebug() << "Aa::getMass() begin";
264  m_mass = m_peptide.get()->getMass();
273 
274 
275  // qDebug() << "Aa::getMass() end " << mass;
276 }
277 
279  : m_peptide(other.m_peptide), m_mapIsotope(other.m_mapIsotope)
280 {
281  m_ratio = other.m_ratio;
282 }
283 
285 {
286 }
287 
288 
291 {
292 
293  return m_mass;
294 }
295 
296 
299 {
300 
301  return m_ratio *
303  Isotope::H2,
304  (m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::H) + charge) -
305  m_peptide.get()->getNumberOfIsotope(Isotope::H2),
307 }
308 
309 
310 int
312 {
313  return m_peptide.get()->getNumberOfAtom(atom);
314 }
315 
316 int
318 {
319  return m_mapIsotope.at(isotope) +
320  m_peptide.get()->getNumberOfIsotope(isotope);
321 }
322 
323 const std::map<Isotope, int> &
325 {
326  return m_mapIsotope;
327 }
328 
329 
330 bool
332 {
333  return m_peptide.get()->isPalindrome();
334 }
335 
336 
337 unsigned int
339 {
340  return m_peptide.get()->size();
341 }
342 
343 const QString
345 {
346  return m_peptide.get()->getSequence();
347 }
348 
349 unsigned int
351 {
352  // only count variable (natural) isotope
356  (m_mapIsotope.at(Isotope::S34) * 2) +
357  (m_mapIsotope.at(Isotope::S36) * 4);
358 }
359 } // namespace pappso
pappso::ratioH2_cache
pappso_double ratioH2_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:92
pappso::AtomIsotopeSurvey::S
@ S
pappso::PeptideNaturalIsotope::m_mass
pappso_double m_mass
Definition: peptidenaturalisotope.h:94
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition: types.h:69
pappso::isotopem_ratio_cache
pappso_double isotopem_ratio_cache(Isotope isotope, unsigned int total, unsigned int heavy)
Definition: peptidenaturalisotope.cpp:95
pappso::ratioS34_cache
pappso_double ratioS34_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:88
pappso::Isotope::H2
@ H2
pappso::ratioO17_cache
pappso_double ratioO17_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:90
pappso::ABUNDANCEO17
const pappso_double ABUNDANCEO17(0.000372)
pappso::ABUNDANCES36
const pappso_double ABUNDANCES36(0.00020)
pappso::Isotope::S36
@ S36
pappso::PeptideNaturalIsotope::getNumberOfIsotope
virtual int getNumberOfIsotope(Isotope isotope) const override
get the number of isotopes C13, H2, O17, O18, N15, S33, S34, S36 in the molecule
Definition: peptidenaturalisotope.cpp:317
pappso::ABUNDANCES33
const pappso_double ABUNDANCES33(0.00750)
pappso::PeptideNaturalIsotope::m_peptide
const PeptideInterfaceSp m_peptide
Definition: peptidenaturalisotope.h:90
pappso::ABUNDANCEH2
const pappso_double ABUNDANCEH2(0.00015574)
pappso::AtomIsotope
AtomIsotope
Definition: peptidenaturalisotope.cpp:70
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks
Definition: aa.cpp:39
pappso::PeptideNaturalIsotope::m_ratio
pappso_double m_ratio
Definition: peptidenaturalisotope.h:93
pappso::DIFFN14N15
const pappso_double DIFFN14N15(15.0001088982 - MASSNITROGEN)
pappso::Isotope
Isotope
Definition: types.h:112
pappso::DIFFH1H2
const pappso_double DIFFH1H2(2.0141017778 - MPROTIUM)
pappso::ratioC13_cache
pappso_double ratioC13_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:85
pappso::ABUNDANCEN15
const pappso_double ABUNDANCEN15(0.003663)
pappso::PeptideNaturalIsotope::getIsotopeMap
const std::map< Isotope, int > & getIsotopeMap() const
Definition: peptidenaturalisotope.cpp:324
pappso::Isotope::N15
@ N15
pappso::PeptideNaturalIsotope
Definition: peptidenaturalisotope.h:67
pappso::Combinations
uint64_t Combinations(unsigned int n, unsigned int k)
Definition: peptidenaturalisotope.cpp:47
pappso::DIFFS32S33
const pappso_double DIFFS32S33(32.97145876 - MASSSULFUR)
pappso::PeptideInterfaceSp
std::shared_ptr< const PeptideInterface > PeptideInterfaceSp
Definition: peptideinterface.h:59
pappso::PeptideNaturalIsotope::getSequence
virtual const QString getSequence() const override
amino acid sequence without modification
Definition: peptidenaturalisotope.cpp:344
pappso::Isotope::O17
@ O17
pappso::PeptideNaturalIsotope::size
virtual unsigned int size() const override
Definition: peptidenaturalisotope.cpp:338
pappso::AtomIsotopeSurvey::H
@ H
peptidenaturalisotope.h
peptide natural isotope model
pappso::isotopem_ratio
pappso_double isotopem_ratio(pappso_double abundance, unsigned int total, unsigned int heavy)
Definition: peptidenaturalisotope.cpp:79
pappso::DIFFS32S34
const pappso_double DIFFS32S34(33.96786690 - MASSSULFUR)
pappso::AtomIsotope::C
@ C
pappso::AtomIsotope::N
@ N
pappso::PeptideNaturalIsotope::getNumberOfAtom
virtual int getNumberOfAtom(AtomIsotopeSurvey atom) const override
get the number of atom C, O, N, H in the molecule
Definition: peptidenaturalisotope.cpp:311
pappso::Isotope::S34
@ S34
pappso::DIFFC12C13
const pappso_double DIFFC12C13(1.0033548378)
pappso::AtomIsotope::H
@ H
pappso::combinations_cache
uint64_t combinations_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:44
pappso::PeptideNaturalIsotope::isPalindrome
virtual bool isPalindrome() const override
tells if the peptide sequence is a palindrome
Definition: peptidenaturalisotope.cpp:331
pappso::PeptideNaturalIsotope::PeptideNaturalIsotope
PeptideNaturalIsotope(const PeptideInterfaceSp &peptide, const std::map< Isotope, int > &map_isotope)
Definition: peptidenaturalisotope.cpp:201
pappso::DIFFS32S36
const pappso_double DIFFS32S36(35.96708076 - MASSSULFUR)
CACHE_ARRAY_SIZE
#define CACHE_ARRAY_SIZE
Definition: peptidenaturalisotope.cpp:42
pappso::AtomIsotope::O
@ O
pappso::ABUNDANCES34
const pappso_double ABUNDANCES34(0.04215)
pappso::ABUNDANCEC13
const pappso_double ABUNDANCEC13(0.011078)
pappso::ratioS36_cache
pappso_double ratioS36_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:87
pappso::PeptideNaturalIsotope::~PeptideNaturalIsotope
~PeptideNaturalIsotope()
Definition: peptidenaturalisotope.cpp:284
pappso::PeptideNaturalIsotope::getMass
pappso_double getMass() const override
Definition: peptidenaturalisotope.cpp:290
pappso::AtomIsotopeSurvey::C
@ C
pappso::Isotope::C13
@ C13
pappso::AtomIsotopeSurvey::N
@ N
pappso::PeptideNaturalIsotope::getIsotopeNumber
virtual unsigned int getIsotopeNumber() const
Definition: peptidenaturalisotope.cpp:350
pappso::ratioN15_cache
pappso_double ratioN15_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:86
pappso::PeptideNaturalIsotope::getIntensityRatio
pappso_double getIntensityRatio(unsigned int charge) const
Definition: peptidenaturalisotope.cpp:298
pappso::AtomIsotope::S
@ S
pappso::PeptideNaturalIsotope::m_mapIsotope
const std::map< Isotope, int > m_mapIsotope
Definition: peptidenaturalisotope.h:92
pappso::Isotope::O18
@ O18
pappso::DIFFO16O17
const pappso_double DIFFO16O17(16.99913150 - MASSOXYGEN)
pappso::ABUNDANCEO18
const pappso_double ABUNDANCEO18(0.0020004)
pappso::ratioO18_cache
pappso_double ratioO18_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:91
pappso::ratioS33_cache
pappso_double ratioS33_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
Definition: peptidenaturalisotope.cpp:89
pappso::DIFFO16O18
const pappso_double DIFFO16O18(17.9991610 - MASSOXYGEN)
pappso::AtomIsotopeSurvey
AtomIsotopeSurvey
Definition: types.h:97
pappso::AtomIsotopeSurvey::O
@ O
pappso::Isotope::S33
@ S33