libpappsomspp
Library for mass spectrometry
peptidenaturalisotopelist.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/peptide/peptidenaturalisotopelist.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// make test ARGS="-V -I 7,7"
31#include <QDebug>
33#include "../pappsoexception.h"
34
35namespace pappso
36{
37
39 const PeptideInterfaceSp &peptide, pappso_double minimum_ratio_to_compute)
40 : msp_peptide(peptide)
41{
42 // qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__ << " "
43 // << minimum_ratio_to_compute;
44 int number_of_fixed_oxygen =
45 msp_peptide.get()->getNumberOfIsotope(Isotope::O18) +
46 msp_peptide.get()->getNumberOfIsotope(Isotope::O17);
47 int number_of_fixed_sulfur =
48 msp_peptide.get()->getNumberOfIsotope(Isotope::S33) +
49 msp_peptide.get()->getNumberOfIsotope(Isotope::S34) +
50 msp_peptide.get()->getNumberOfIsotope(Isotope::S36);
51 int number_of_fixed_nitrogen =
52 msp_peptide.get()->getNumberOfIsotope(Isotope::N15);
53 int number_of_fixed_hydrogen =
54 msp_peptide.get()->getNumberOfIsotope(Isotope::H2);
55
56 int total_carbon(msp_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::C) -
57 msp_peptide.get()->getNumberOfIsotope(Isotope::C13));
58
59 // qDebug() << "total_carbon " << total_carbon;
60 // qDebug() << "total_sulfur " << total_sulfur;
61 std::map<Isotope, int> map_isotope;
62 map_isotope.insert(std::pair<Isotope, int>(Isotope::C13, 0));
63 map_isotope.insert(std::pair<Isotope, int>(Isotope::H2, 0));
64 map_isotope.insert(std::pair<Isotope, int>(Isotope::N15, 0));
65 map_isotope.insert(std::pair<Isotope, int>(Isotope::O17, 0));
66 map_isotope.insert(std::pair<Isotope, int>(Isotope::O18, 0));
67 map_isotope.insert(std::pair<Isotope, int>(Isotope::S33, 0));
68 map_isotope.insert(std::pair<Isotope, int>(Isotope::S34, 0));
69 map_isotope.insert(std::pair<Isotope, int>(Isotope::S36, 0));
70
71 for(int nbc13 = 0; nbc13 <= total_carbon; nbc13++)
72 {
73 map_isotope[Isotope::C13] = nbc13;
74 PeptideNaturalIsotopeSp pepIsotope =
75 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
76 this->msp_peptide_natural_isotope_list.push_back(pepIsotope);
77 if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
78 {
79 break;
80 }
81 }
82 std::list<PeptideNaturalIsotopeSp> temp_list;
83
84 // ******************************************************************
85 // Sulfur isotope list
86 int total_sulfur(msp_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::S) -
87 number_of_fixed_sulfur);
88 std::list<PeptideNaturalIsotopeSp>::iterator it =
90 while(it != msp_peptide_natural_isotope_list.end())
91 {
92 map_isotope = it->get()->getIsotopeMap();
93 for(int nbS34 = 1; nbS34 <= total_sulfur; nbS34++)
94 {
95 map_isotope[Isotope::S34] = nbS34;
96 PeptideNaturalIsotopeSp pepIsotope =
97 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
98 temp_list.push_back(pepIsotope);
99 if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
100 {
101 // qDebug() << "peptide " << pepIsotope.get()->getFormula(1) << "
102 // " << pepIsotope.get()->getIntensityRatio(1); it++;
103 break;
104 }
105 }
106 it++;
107 }
109 it, temp_list.begin(), temp_list.end());
110
111 // compute S33 abundance
112 temp_list.resize(0);
114 while(it != msp_peptide_natural_isotope_list.end())
115 {
116 // qDebug() << "peptide S33 " << it->get()->getFormula(1) << " "
117 // <<it->get()->getIntensityRatio(1);
118 map_isotope = it->get()->getIsotopeMap();
119 for(int nbS33 = 1; nbS33 <= (total_sulfur - map_isotope[Isotope::S34]);
120 nbS33++)
121 {
122 map_isotope[Isotope::S33] = nbS33;
123 PeptideNaturalIsotopeSp pepIsotopeS33 =
124 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
125 temp_list.push_back(pepIsotopeS33);
126 if(pepIsotopeS33.get()->getIntensityRatio(1) <
127 minimum_ratio_to_compute)
128 {
129 // it++;
130 break;
131 }
132 }
133 it++;
134 }
136 it, temp_list.begin(), temp_list.end());
137
138 // compute S36 abundance
139 temp_list.resize(0);
141 while(it != msp_peptide_natural_isotope_list.end())
142 {
143 map_isotope = it->get()->getIsotopeMap();
144 for(int nbS36 = 1; nbS36 <= (total_sulfur - map_isotope[Isotope::S34] -
145 map_isotope[Isotope::S33]);
146 nbS36++)
147 {
148 map_isotope[Isotope::S36] = nbS36;
149 PeptideNaturalIsotopeSp pepIsotopeS36 =
150 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
151 temp_list.push_back(pepIsotopeS36);
152 if(pepIsotopeS36.get()->getIntensityRatio(1) <
153 minimum_ratio_to_compute)
154 {
155 // it++;
156 break;
157 }
158 }
159 it++;
160 }
162 it, temp_list.begin(), temp_list.end());
163
164 // ******************************************************************
165
166
167 // ******************************************************************
168 // Hydrogen isotope list
169 temp_list.resize(0);
170 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
171 // total_hydrogen";
172 int total_hydrogen(msp_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::H) -
173 number_of_fixed_hydrogen);
175 while(it != msp_peptide_natural_isotope_list.end())
176 {
177 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
178 // getIsotopeMap " << it->getFormula(1) << " " <<
179 // msp_peptide_natural_isotope_list.size();
180 map_isotope = it->get()->getIsotopeMap();
181 for(int nbH2 = 1; nbH2 <= total_hydrogen; nbH2++)
182 {
183 map_isotope[Isotope::H2] = nbH2;
184 PeptideNaturalIsotopeSp pepIsotope =
185 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
186 temp_list.push_back(pepIsotope);
187 if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
188 {
189 // it++;
190 break;
191 }
192 }
193 it++;
194 }
196 it, temp_list.begin(), temp_list.end());
197 // ******************************************************************
198
199
200 // ******************************************************************
201 // Oxygen isotope list
202 temp_list.resize(0);
203 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
204 // total_oxygen";
205 unsigned int total_oxygen(
206 msp_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::O) -
207 number_of_fixed_oxygen);
209 while(it != msp_peptide_natural_isotope_list.end())
210 {
211 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
212 // getIsotopeMap " << it->getFormula(1) << " " <<
213 // msp_peptide_natural_isotope_list.size();
214 map_isotope = it->get()->getIsotopeMap();
215 for(unsigned int nbO18 = 1; nbO18 <= total_oxygen; nbO18++)
216 {
217 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
218 // nbO18 " << nbO18;
219 map_isotope[Isotope::O18] = nbO18;
220 PeptideNaturalIsotopeSp pepIsotope =
221 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
222 temp_list.push_back(pepIsotope);
223 if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
224 {
225 // it++;
226 break;
227 }
228 }
229 it++;
230 }
232 it, temp_list.begin(), temp_list.end());
233 // ******************************************************************
234
235
236 // ******************************************************************
237 // Nitrogen isotope list
238 temp_list.resize(0);
239 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
240 // total_nitrogen";
241 unsigned int total_nitrogen(
242 msp_peptide.get()->getNumberOfAtom(AtomIsotopeSurvey::N) -
243 number_of_fixed_nitrogen);
245 while(it != msp_peptide_natural_isotope_list.end())
246 {
247 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
248 // getIsotopeMap " << it->getFormula(1) << " " <<
249 // msp_peptide_natural_isotope_list.size();
250 map_isotope = it->get()->getIsotopeMap();
251 for(unsigned int nbN15 = 1; nbN15 <= total_nitrogen; nbN15++)
252 {
253 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList
254 // nbN15 " << nbN15;
255 map_isotope[Isotope::N15] = nbN15;
256 PeptideNaturalIsotopeSp pepIsotope =
257 std::make_shared<PeptideNaturalIsotope>(msp_peptide, map_isotope);
258 temp_list.push_back(pepIsotope);
259 if(pepIsotope.get()->getIntensityRatio(1) < minimum_ratio_to_compute)
260 {
261 // it++;
262 break;
263 }
264 }
265 it++;
266 }
268 it, temp_list.begin(), temp_list.end());
269 // ******************************************************************
270
271 // qDebug() << "PeptideNaturalIsotopeList::PeptideNaturalIsotopeList end
272 // size="<<msp_peptide_natural_isotope_list.size();
273}
274
277{
278 return std::make_shared<PeptideNaturalIsotopeList>(*this);
279}
280
282 const PeptideNaturalIsotopeList &other)
283 : msp_peptide(other.msp_peptide),
284 msp_peptide_natural_isotope_list(other.msp_peptide_natural_isotope_list)
285{
286}
287
289{
290}
291
292const std::map<unsigned int, pappso_double>
294{
295 std::list<PeptideNaturalIsotopeSp>::const_iterator it =
297 std::map<unsigned int, pappso_double> map_isotope_number;
298
299 while(it != msp_peptide_natural_isotope_list.end())
300 {
301 unsigned int number = it->get()->getIsotopeNumber();
302 std::pair<std::map<unsigned int, pappso_double>::iterator, bool> mapnew =
303 map_isotope_number.insert(
304 std::pair<unsigned int, pappso_double>(number, 0));
305 if(mapnew.second == false)
306 {
307 // mapit = map_isotope_number.insert(std::pair<unsigned
308 // int,pappso_double>(number, 0));
309 }
310 mapnew.first->second += it->get()->getIntensityRatio(1);
311 it++;
312 }
313 return map_isotope_number;
314}
315
316
317/** /brief get a sorted (by expected intensity) vector of isotopes of the same
318 * level
319 *
320 * */
321
322std::vector<PeptideNaturalIsotopeSp>
324 unsigned int charge) const
325{
326 std::vector<PeptideNaturalIsotopeSp> v_isotope_list;
327
328 for(auto &&isotopeSp : msp_peptide_natural_isotope_list)
329 {
330 if(isotopeSp.get()->getIsotopeNumber() == isotope_number)
331 {
332 v_isotope_list.push_back(isotopeSp);
333 }
334 }
335 std::sort(v_isotope_list.begin(),
336 v_isotope_list.end(),
337 [charge](const PeptideNaturalIsotopeSp &m,
338 const PeptideNaturalIsotopeSp &n) {
339 return (m.get()->getIntensityRatio(charge) >
340 n.get()->getIntensityRatio(charge));
341 });
342 return v_isotope_list;
343}
344
345
346/** /brief get a sorted (by expected intensity) vector of natural isotope
347 * average by isotope number
348 *
349 * */
350std::vector<PeptideNaturalIsotopeAverageSp>
352 unsigned int charge,
353 PrecisionPtr precision,
354 unsigned int isotopeNumber,
355 pappso_double minimumIntensity)
356{
357
358 // qDebug() << "getByIntensityRatioByIsotopeNumber begin";
359 unsigned int askedIsotopeRank;
360 unsigned int maxAskedIsotopeRank = 10;
361 pappso_double cumulativeRatio = 0;
362 std::vector<PeptideNaturalIsotopeAverageSp> v_isotopeAverageList;
363 std::vector<PeptideNaturalIsotopeAverageSp> v_isotopeAverageListResult;
364
365 std::vector<unsigned int> previousIsotopeRank;
366 bool isEmpty = false;
367 for(askedIsotopeRank = 1;
368 (askedIsotopeRank < maxAskedIsotopeRank) && (!isEmpty);
369 askedIsotopeRank++)
370 {
371 PeptideNaturalIsotopeAverage isotopeAverage(
372 peptide, askedIsotopeRank, isotopeNumber, charge, precision);
373 isEmpty = isotopeAverage.isEmpty();
374 if(isEmpty)
375 {
376 }
377 else
378 {
379 if(std::find(previousIsotopeRank.begin(),
380 previousIsotopeRank.end(),
381 isotopeAverage.getIsotopeRank()) ==
382 previousIsotopeRank.end())
383 { // not Found
384 previousIsotopeRank.push_back(isotopeAverage.getIsotopeRank());
385 v_isotopeAverageList.push_back(
386 isotopeAverage.makePeptideNaturalIsotopeAverageSp());
387 }
388 }
389 }
390 if(v_isotopeAverageList.size() == 0)
391 return v_isotopeAverageListResult;
392
393 // qDebug() << "getByIntensityRatioByIsotopeNumber comp";
394 std::sort(v_isotopeAverageList.begin(),
395 v_isotopeAverageList.end(),
398 return (m.get()->getIntensityRatio() >
399 n.get()->getIntensityRatio());
400 });
401
402 cumulativeRatio = 0;
403 auto it = v_isotopeAverageList.begin();
404 v_isotopeAverageListResult.clear();
405 // qDebug() << "getByIntensityRatioByIsotopeNumber cumul";
406 while((it != v_isotopeAverageList.end()) &&
407 (cumulativeRatio < minimumIntensity))
408 {
409 cumulativeRatio += it->get()->getIntensityRatio();
410 v_isotopeAverageListResult.push_back(*it);
411 it++;
412 }
413
414 // qDebug() << "getByIntensityRatioByIsotopeNumber end";
415
416 return v_isotopeAverageListResult;
417}
418
419
420/** /brief get a sorted (by expected intensity) vector of natural isotope
421 * average
422 *
423 * */
424std::vector<PeptideNaturalIsotopeAverageSp>
426 unsigned int charge,
427 PrecisionPtr precision,
428 pappso_double minimumIntensityRatio) const
429{
430
431 // qDebug() << "PeptideNaturalIsotopeList::getByIntensityRatio begin";
432 std::vector<PeptideNaturalIsotopeAverageSp>
433 peptide_natural_isotope_average_list;
434
435 std::map<unsigned int, pappso::pappso_double> map_isotope_number =
437 std::vector<std::pair<unsigned int, pappso::pappso_double>>
438 sorted_number_ratio;
439
440 for(unsigned int i = 0; i < map_isotope_number.size(); i++)
441 {
442 sorted_number_ratio.push_back(
443 std::pair<unsigned int, pappso::pappso_double>(i,
444 map_isotope_number[i]));
445 unsigned int asked_rank = 0;
446 unsigned int given_rank = 0;
447 bool more_rank = true;
448 while(more_rank)
449 {
450 asked_rank++;
451 pappso::PeptideNaturalIsotopeAverage isotopeAverageMono(
452 *this, asked_rank, i, charge, precision);
453 given_rank = isotopeAverageMono.getIsotopeRank();
454 if(given_rank < asked_rank)
455 {
456 more_rank = false;
457 }
458 else if(isotopeAverageMono.getIntensityRatio() == 0)
459 {
460 more_rank = false;
461 }
462 else
463 {
464 // isotopeAverageMono.makePeptideNaturalIsotopeAverageSp();
465 peptide_natural_isotope_average_list.push_back(
466 isotopeAverageMono.makePeptideNaturalIsotopeAverageSp());
467 }
468 }
469 }
470
471
472 // sort by intensity ratio
473 std::sort(sorted_number_ratio.begin(),
474 sorted_number_ratio.end(),
475 [](const std::pair<unsigned int, pappso::pappso_double> &m,
476 const std::pair<unsigned int, pappso::pappso_double> &n) {
477 return (m.second > n.second);
478 });
479
480
481 double cumulativeRatio = 0;
482 std::vector<unsigned int> selected_isotope_number_list;
483 for(auto &pair_isotope_number : sorted_number_ratio)
484 {
485 if(cumulativeRatio <= minimumIntensityRatio)
486 {
487 selected_isotope_number_list.push_back(pair_isotope_number.first);
488 }
489 else
490 {
491 break;
492 }
493 cumulativeRatio += pair_isotope_number.second;
494 }
495
496 auto it_remove =
497 std::remove_if(peptide_natural_isotope_average_list.begin(),
498 peptide_natural_isotope_average_list.end(),
499 [selected_isotope_number_list](
500 const PeptideNaturalIsotopeAverageSp &average) {
501 auto it = std::find(selected_isotope_number_list.begin(),
502 selected_isotope_number_list.end(),
503 average.get()->getIsotopeNumber());
504 return (it == selected_isotope_number_list.end());
505 });
506 peptide_natural_isotope_average_list.erase(
507 it_remove, peptide_natural_isotope_average_list.end());
508 return peptide_natural_isotope_average_list;
509}
510
511
514{
516}
517
520{
522}
523
524unsigned int
526{
528}
529
530const PeptideInterfaceSp &
532{
533 return msp_peptide;
534}
535} // namespace pappso
PeptideNaturalIsotopeAverageSp makePeptideNaturalIsotopeAverageSp() const
const PeptideInterfaceSp & getPeptideInterfaceSp() const
PeptideNaturalIsotopeList(const PeptideInterfaceSp &peptide, pappso_double minimum_ratio_to_compute=0.001)
compute the list of possible isotopes for a peptide
std::list< PeptideNaturalIsotopeSp >::const_iterator const_iterator
PeptideNaturalIsotopeListSp makePeptideNaturalIsotopeListSp() const
std::vector< PeptideNaturalIsotopeAverageSp > getByIntensityRatio(unsigned int charge, PrecisionPtr precision, pappso_double minimum_isotope_pattern_ratio) const
get the list of natural isotopes representing at least a minimum ratio of the whole isotope pattern
std::list< PeptideNaturalIsotopeSp > msp_peptide_natural_isotope_list
const std::map< unsigned int, pappso_double > getIntensityRatioPerIsotopeNumber() const
std::vector< PeptideNaturalIsotopeSp > getByIsotopeNumber(unsigned int isotopeLevel, unsigned int charge) const
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
std::shared_ptr< const PeptideNaturalIsotope > PeptideNaturalIsotopeSp
std::vector< PeptideNaturalIsotopeAverageSp > getByIntensityRatioByIsotopeNumber(const PeptideInterfaceSp &peptide, unsigned int charge, PrecisionPtr precision, unsigned int isotopeNumber, pappso_double minimumIntensity)
std::shared_ptr< const PeptideInterface > PeptideInterfaceSp
double pappso_double
A type definition for doubles.
Definition: types.h:49
std::shared_ptr< const PeptideNaturalIsotopeAverage > PeptideNaturalIsotopeAverageSp
std::shared_ptr< const PeptideNaturalIsotopeList > PeptideNaturalIsotopeListSp
peptide natural isotope model