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
virtual int getNumberOfAtom(AtomIsotopeSurvey atom) const override
get the number of atom C, O, N, H in the molecule
virtual const QString getSequence() const override
amino acid sequence without modification
virtual int getNumberOfIsotope(Isotope isotope) const override
get the number of isotopes C13, H2, O17, O18, N15, S33, S34, S36 in the molecule
const std::map< Isotope, int > & getIsotopeMap() const
const PeptideInterfaceSp m_peptide
PeptideNaturalIsotope(const PeptideInterfaceSp &peptide, const std::map< Isotope, int > &map_isotope)
virtual bool isPalindrome() const override
tells if the peptide sequence is a palindrome
virtual unsigned int size() const override
virtual unsigned int getIsotopeNumber() const
const std::map< Isotope, int > m_mapIsotope
pappso_double getIntensityRatio(unsigned int charge) const
pappso_double getMass() const override
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
const pappso_double DIFFS32S33(32.9714589101 - MASSSULFUR)
uint64_t Combinations(unsigned int n, unsigned int k)
const pappso_double DIFFS32S34(33.9678670300 - MASSSULFUR)
pappso_double ratioO17_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double DIFFO16O17(16.99913150 - MASSOXYGEN)
const pappso_double ABUNDANCES36(0.00010999120070394368536836893213148869108408689498901367187500)
std::shared_ptr< const PeptideInterface > PeptideInterfaceSp
pappso_double ratioS34_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double ABUNDANCEN15(0.00364198543205827118818262988497735932469367980957031250000000)
const pappso_double DIFFS32S36(35.9670812000 - MASSSULFUR)
pappso_double ratioC13_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
AtomIsotopeSurvey
Definition: types.h:77
double pappso_double
A type definition for doubles.
Definition: types.h:49
Isotope
Definition: types.h:92
pappso_double isotopem_ratio(pappso_double abundance, unsigned int total, unsigned int heavy)
pappso_double ratioH2_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double ABUNDANCEC13(0.01078805814953308406245469086570665240287780761718750000000000)
const pappso_double ABUNDANCEO17(0.00038099847600609595965615028489992255344986915588378906250000)
pappso_double isotopem_ratio_cache(Isotope isotope, unsigned int total, unsigned int heavy)
const pappso_double ABUNDANCEH2(0.00011570983569203332000374651045149221317842602729797363281250)
pappso_double ratioO18_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double ABUNDANCES34(0.04252059835213182203972337447339668869972229003906250000000000)
pappso_double ratioS36_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double ABUNDANCEO18(0.00205139179443282221315669744399201590567827224731445312500000)
uint64_t combinations_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double DIFFO16O18(17.9991610 - MASSOXYGEN)
const pappso_double ABUNDANCES33(0.00751939844812414937003097747947322204709053039550781250000000)
const pappso_double DIFFN14N15(15.0001088982 - MASSNITROGEN)
pappso_double ratioN15_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
pappso_double ratioS33_cache[CACHE_ARRAY_SIZE][CACHE_ARRAY_SIZE]
const pappso_double DIFFC12C13(1.0033548378)
const pappso_double DIFFH1H2(2.0141017778 - MPROTIUM)
#define CACHE_ARRAY_SIZE
peptide natural isotope model