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 
35 namespace 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 
292 const 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 
322 std::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  * */
350 std::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(),
396  [](const PeptideNaturalIsotopeAverageSp &m,
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  * */
424 std::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 {
515  return msp_peptide_natural_isotope_list.begin();
516 }
517 
520 {
522 }
523 
524 unsigned int
526 {
527  return msp_peptide_natural_isotope_list.size();
528 }
529 
530 const 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::shared_ptr< const PeptideInterface > PeptideInterfaceSp
double pappso_double
A type definition for doubles.
Definition: types.h:49
std::shared_ptr< const PeptideNaturalIsotopeAverage > PeptideNaturalIsotopeAverageSp
std::vector< PeptideNaturalIsotopeAverageSp > getByIntensityRatioByIsotopeNumber(const PeptideInterfaceSp &peptide, unsigned int charge, PrecisionPtr precision, unsigned int isotopeNumber, pappso_double minimumIntensity)
std::shared_ptr< const PeptideNaturalIsotopeList > PeptideNaturalIsotopeListSp
peptide natural isotope model