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
32#include "../pappsoexception.h"
33
34#include <cmath>
35#include <QDebug>
36
37using namespace std;
38
39namespace pappso
40{
41
42#define CACHE_ARRAY_SIZE 500
43
45
46uint64_t
47Combinations(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
69enum class AtomIsotope
70{
71 C,
72 H,
73 O,
74 N,
75 S
76};
77
79isotopem_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
95isotopem_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
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
228 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::C) -
229 m_peptide.get()->getNumberOfIsotope(Isotope::C13),
233 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::N) -
234 m_peptide.get()->getNumberOfIsotope(Isotope::N15),
238 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
239 number_of_fixed_oxygen,
243 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
244 number_of_fixed_oxygen,
248 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
249 number_of_fixed_sulfur,
253 m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
254 number_of_fixed_sulfur,
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 *
304 (m_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::H) + charge) -
305 m_peptide.get()->getNumberOfIsotope(Isotope::H2),
307}
308
309
310int
312{
313 return m_peptide.get()->getNumberOfAtom(atom);
314}
315
316int
318{
319 return m_mapIsotope.at(isotope) +
320 m_peptide.get()->getNumberOfIsotope(isotope);
321}
322
323const std::map<Isotope, int> &
325{
326 return m_mapIsotope;
327}
328
329
330bool
332{
333 return m_peptide.get()->isPalindrome();
334}
335
336
337unsigned int
339{
340 return m_peptide.get()->size();
341}
342
343const QString
345{
346 return m_peptide.get()->getSequence();
347}
348
349unsigned 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:89
double pappso_double
A type definition for doubles.
Definition: types.h:50
Isotope
Definition: types.h:104
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