libpappsomspp
Library for mass spectrometry
psmfeatures.cpp
Go to the documentation of this file.
1 /**
2  * \file pappsomspp/psm/features/psmfeatures.cpp
3  * \date 19/07/2022
4  * \author Olivier Langella
5  * \brief comutes various PSM (Peptide Spectrum Match) features
6  */
7 
8 /*******************************************************************************
9  * Copyright (c) 2022 Olivier Langella <Olivier.Langella@u-psud.fr>.
10  *
11  * This file is part of PAPPSOms-tools.
12  *
13  * PAPPSOms-tools 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-tools 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-tools. If not, see <http://www.gnu.org/licenses/>.
25  *
26  ******************************************************************************/
27 
28 
29 #include "psmfeatures.h"
30 #include <memory>
31 #include <cmath>
32 
33 using namespace pappso;
34 
35 PsmFeatures::PsmFeatures(PrecisionPtr ms2precision, double minimumMz)
36 {
37 
38  m_ms2precision = ms2precision;
39 
40  m_ionList.push_back(PeptideIon::y);
41  m_ionList.push_back(PeptideIon::b);
42 
43 
45  std::make_shared<FilterResampleKeepGreater>(minimumMz);
46 
48  std::make_shared<FilterChargeDeconvolution>(m_ms2precision);
49  msp_filterMzExclusion = std::make_shared<FilterMzExclusion>(
51 }
52 
54 {
55 }
56 
57 void
59  const MassSpectrum *p_spectrum,
60  unsigned int parent_charge)
61 {
62  m_peakIonPairs.clear();
63  msp_peptide = peptideSp;
64  MassSpectrum spectrum(*p_spectrum);
65  msp_filterKeepGreater.get()->filter(spectrum);
66  // msp_filterChargeDeconvolution.get()->filter(spectrum);
67  // msp_filterMzExclusion.get()->filter(spectrum);
68 
70  std::make_shared<PeptideIsotopeSpectrumMatch>(spectrum,
71  peptideSp,
72  parent_charge,
74  m_ionList,
75  (unsigned int)1,
76  0);
77 
78  msp_peptideSpectrumMatch.get()->dropPeaksLackingMonoisotope();
79  m_spectrumSumIntensity = spectrum.sumY();
80 
81 
82  qDebug() << " accumulate";
83  std::vector<double> delta_list;
84 
85 
86  // TODO compute number of matched complementary peaks having m/z compatible
87  // with the precursor
88 
89  m_precursorTheoreticalMz = peptideSp.get()->getMz(parent_charge);
90  m_precursorTheoreticalMass = peptideSp.get()->getMass();
91  m_parentCharge = parent_charge;
92 
93 
94  findComplementIonPairs(peptideSp);
95 
96 
97  for(const pappso::PeakIonIsotopeMatch &peak_ion :
98  msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
99  {
100  delta_list.push_back(
101  peak_ion.getPeptideFragmentIonSp().get()->getMz(peak_ion.getCharge()) -
102  peak_ion.getPeak().x);
103  }
105  std::accumulate(delta_list.begin(), delta_list.end(), 0);
106 
107  qDebug() << " delta_list.size()=" << delta_list.size();
110  m_matchedMzDiffSd = 0;
111  if(delta_list.size() > 0)
112  {
113  m_matchedMzDiffMean = sum / ((pappso::pappso_double)delta_list.size());
114 
115  std::sort(delta_list.begin(), delta_list.end());
116  m_matchedMzDiffMedian = delta_list[(delta_list.size() / 2)];
117 
118 
119  qDebug() << " sd";
120  m_matchedMzDiffSd = 0;
121  for(pappso::pappso_double val : delta_list)
122  {
123  // sd = sd + ((val - mean) * (val - mean));
124  m_matchedMzDiffSd += std::pow((val - m_matchedMzDiffMean), 2);
125  }
126  m_matchedMzDiffSd = m_matchedMzDiffSd / delta_list.size();
128  }
129  else
130  {
131  }
132 }
133 
134 
135 double
137 {
138  double sum = 0;
139  for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
140  {
141  if(peak_ion_match.getPeptideIonType() == ion_type)
142  {
143  sum += peak_ion_match.getPeak().y;
144  }
145  }
146  return sum;
147 }
148 double
150 {
151  double sum = 0;
152  for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
153  {
154  sum += peak_ion_match.getPeak().y;
155  }
156  return sum;
157 }
158 
159 double
161 {
162  return m_spectrumSumIntensity;
163 }
164 
165 std::size_t
167 {
168  return m_peakIonPairs.size();
169 }
170 
171 const std::vector<
172  std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>> &
174 {
175  return m_peakIonPairs;
176 }
177 
178 double
180 {
181 
182  double sum = 0;
183  for(auto &peak_pairs : m_peakIonPairs)
184  {
185  sum += peak_pairs.first.getPeak().y;
186  sum += peak_pairs.second.getPeak().y;
187  }
188  return sum;
189 }
190 
191 double
193 {
194  return m_matchedMzDiffSd;
195 }
196 
197 double
199 {
200  return m_matchedMzDiffMean;
201 }
202 
203 
204 std::size_t
206 {
207  return msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().size();
208 }
209 
210 std::size_t
212 {
213  std::size_t max = 0;
214  auto peak_ion_match_list =
215  msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
216 
217  peak_ion_match_list.erase(
218  std::remove_if(
219  peak_ion_match_list.begin(),
220  peak_ion_match_list.end(),
221  [ion_type](const PeakIonIsotopeMatch &a) {
222  if(a.getPeptideIonType() != ion_type)
223  return true;
224  if(a.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() > 0)
225  return true;
226  return false;
227  }),
228  peak_ion_match_list.end());
229 
230  peak_ion_match_list.sort(
231  [](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
232  if(a.getCharge() < b.getCharge())
233  return true;
234  if(a.getPeptideIonType() < b.getPeptideIonType())
235  return true;
236  if(a.getPeptideFragmentIonSp().get()->size() <
237  b.getPeptideFragmentIonSp().get()->size())
238  return true;
239  return false;
240  });
241 
242  unsigned int charge = 0;
243  std::size_t size = 0;
244  std::size_t count = 0;
245  for(std::list<PeakIonIsotopeMatch>::iterator it = peak_ion_match_list.begin();
246  it != peak_ion_match_list.end();
247  it++)
248  {
249  qDebug()
250  << it->toString() << max << " " << it->getPeak().x << " "
251  << it->getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber();
252  count++;
253  if((charge != it->getCharge()) ||
254  (size != (it->getPeptideFragmentIonSp().get()->size() - 1)))
255  {
256  count = 1;
257  charge = it->getCharge();
258  }
259  if(max < count)
260  max = count;
261 
262  size = it->getPeptideFragmentIonSp().get()->size();
263  }
264 
265  return max;
266 }
267 
268 std::size_t
270 {
271  std::vector<bool> covered;
272  covered.resize(msp_peptide.get()->size(), false);
273 
274  for(auto &peak : msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
275  {
276  if(peak.getPeptideIonType() == ion_type)
277  {
278  covered[peak.getPeptideFragmentIonSp().get()->size() - 1] = true;
279  }
280  }
281  return std::count(covered.begin(), covered.end(), true);
282 }
283 
284 std::size_t
286 {
287 
288  std::vector<bool> covered;
289  covered.resize(msp_peptide.get()->size(), false);
290 
291  for(auto &peak_pair : m_peakIonPairs)
292  {
293  std::size_t pos =
294  peak_pair.first.getPeptideFragmentIonSp().get()->size() - 1;
295  covered[pos] = true;
296  covered[pos + 1] = true;
297  }
298  return std::count(covered.begin(), covered.end(), true);
299 }
300 
301 
302 double
304  pappso::PeptideIon ion_type) const
305 {
306  std::list<pappso::PeakIonIsotopeMatch> peak_ion_type =
307  msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
308 
309  peak_ion_type.remove_if([ion_type](const PeakIonIsotopeMatch &a) {
310  return (a.getPeptideIonType() != ion_type);
311  });
312  auto peak_it = std::max_element(
313  peak_ion_type.begin(),
314  peak_ion_type.end(),
315  [](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
316  return (a.getPeak().y < b.getPeak().y);
317  });
318 
319  if(peak_it == peak_ion_type.end())
320  return 0;
321  return peak_it->getPeak().y;
322 }
323 
324 double
326  const
327 {
328  auto peak_it = std::max_element(
329  m_peakIonPairs.begin(),
330  m_peakIonPairs.end(),
331  [](const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
332  &a,
333  const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
334  &b) {
335  return ((a.first.getPeak().y + a.second.getPeak().y) <
336  (b.first.getPeak().y + b.second.getPeak().y));
337  });
338 
339  if(peak_it == m_peakIonPairs.end())
340  return 0;
341 
342  return getIonPairPrecursorMassDelta(*peak_it);
343 }
344 
345 double
347  const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
348  &ion_pair) const
349 {
350  qDebug() << m_precursorTheoreticalMz << " " << ion_pair.first.getPeak().x
351  << " " << ion_pair.second.getPeak().x << " "
352  << ion_pair.second.getCharge() << " " << ion_pair.first.getCharge()
353  << " " << m_parentCharge;
354  double diff =
355  (m_precursorTheoreticalMass + (MHPLUS * ion_pair.first.getCharge())) /
356  ion_pair.first.getCharge();
357 
358 
359  return (diff - (ion_pair.first.getPeak().x + ion_pair.second.getPeak().x -
360  ((MHPLUS * ion_pair.first.getCharge())) /
361  ion_pair.first.getCharge()));
362 }
363 
364 
365 void
367 {
368  std::size_t peptide_size = peptideSp.get()->size();
369  std::vector<PeakIonIsotopeMatch> ion_isotope_list(
370  msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().begin(),
371  msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().end());
372  for(const pappso::PeakIonIsotopeMatch &peak_ion_ext : ion_isotope_list)
373  {
374  if(peptideIonIsNter(peak_ion_ext.getPeptideIonType()))
375  {
376  auto it = findComplementIonType(ion_isotope_list.begin(),
377  ion_isotope_list.end(),
378  peak_ion_ext,
379  peptide_size);
380  if(it != ion_isotope_list.end())
381  { // contains the complementary ion
382  m_peakIonPairs.push_back({peak_ion_ext, *it});
383  }
384  }
385  }
386 }
Class to represent a mass spectrum.
Definition: massspectrum.h:71
static PrecisionPtr getPrecisionPtrFractionInstance(PrecisionPtr origin, double fraction)
get the fraction of an existing precision pointer
Definition: precision.cpp:203
std::size_t getNumberOfMatchedIons() const
number of matched ions (peaks)
void findComplementIonPairs(const pappso::PeptideSp &peptideSp)
double getTotalIntensity() const
sum of all peak intensities (matched or not)
double getMatchedMzDiffMean() const
get mean deviation of matched peak mass delta
double getTotalIntensityOfMatchedIonComplementPairs() const
intensity of matched ion complement
std::vector< std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > > m_peakIonPairs
Definition: psmfeatures.h:155
std::shared_ptr< FilterChargeDeconvolution > msp_filterChargeDeconvolution
Definition: psmfeatures.h:137
void setPeptideSpectrumCharge(const pappso::PeptideSp peptideSp, const MassSpectrum *p_spectrum, unsigned int parent_charge)
Definition: psmfeatures.cpp:58
const std::vector< std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > > & getPeakIonPairs() const
std::shared_ptr< FilterResampleKeepGreater > msp_filterKeepGreater
Definition: psmfeatures.h:139
std::shared_ptr< FilterMzExclusion > msp_filterMzExclusion
Definition: psmfeatures.h:138
double getIonPairPrecursorMassDelta(const std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > &ion_pair) const
double m_precursorTheoreticalMz
Definition: psmfeatures.h:149
double getMaxIntensityMatchedIonComplementPairPrecursorMassDelta() const
get the precursor mass delta of the maximum intensity pair of complement ions
std::size_t countMatchedIonComplementPairs() const
count the number of matched ion complement
std::size_t getComplementPairsAaSequenceCoverage()
number of amino acid covered by matched complement pairs of ions
double m_matchedMzDiffMean
Definition: psmfeatures.h:157
double getTotalIntensityOfMatchedIons() const
sum of matched peak intensities
unsigned int m_parentCharge
Definition: psmfeatures.h:151
double getMaxIntensityPeakIonMatch(PeptideIon ion_type) const
std::list< PeptideIon > m_ionList
Definition: psmfeatures.h:145
double m_matchedMzDiffMedian
Definition: psmfeatures.h:158
std::shared_ptr< PeptideIsotopeSpectrumMatch > msp_peptideSpectrumMatch
Definition: psmfeatures.h:141
std::size_t getMaxConsecutiveIon(PeptideIon ion_type)
get the maximum consecutive fragments of one ion type
PrecisionPtr m_ms2precision
Definition: psmfeatures.h:144
double getIntensityOfMatchedIon(PeptideIon ion_type)
get the sum of intensity of a specific ion
pappso::PeptideSp msp_peptide
Definition: psmfeatures.h:142
double m_precursorTheoreticalMass
Definition: psmfeatures.h:150
PsmFeatures(PrecisionPtr ms2precision, double minimumMz)
compute psm features
Definition: psmfeatures.cpp:35
double getMatchedMzDiffSd() const
get standard deviation of matched peak mass delta
double m_spectrumSumIntensity
Definition: psmfeatures.h:147
std::size_t getAaSequenceCoverage(PeptideIon ion_type)
number of amino acid covered by matched ions
pappso_double sumY() const
Definition: trace.cpp:906
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
PeptideIon
PeptideIon enum defines all types of ions (Nter or Cter)
Definition: types.h:386
@ y
Cter amino ions.
@ b
Nter acylium ions.
std::shared_ptr< const Peptide > PeptideSp
const pappso_double MHPLUS(1.007276466879)
double pappso_double
A type definition for doubles.
Definition: types.h:49
std::vector< PeakIonIsotopeMatch >::iterator findComplementIonType(std::vector< PeakIonIsotopeMatch >::iterator begin, std::vector< PeakIonIsotopeMatch >::iterator end, const PeakIonIsotopeMatch &peak_ion, std::size_t peptide_size)
find the first element containing the complementary ion complementary ion of y1 is b(n-1) for instanc...
bool peptideIonIsNter(PeptideIon ion_type)
tells if an ion is Nter
Definition: peptide.cpp:60
@ sum
sum of intensities
@ max
maximum of intensities
comutes various PSM (Peptide Spectrum Match) features