libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
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
33using namespace pappso;
34
35PsmFeatures::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
56
57void
59 const MassSpectrum *p_spectrum,
60 unsigned int parent_charge,
61 unsigned int max_isotope_number)
62{
63 m_peakIonPairs.clear();
64 msp_peptide = peptideSp;
65 MassSpectrum spectrum(*p_spectrum);
66 msp_filterKeepGreater.get()->filter(spectrum);
67 // msp_filterChargeDeconvolution.get()->filter(spectrum);
68 // msp_filterMzExclusion.get()->filter(spectrum);
69
71 std::make_shared<PeptideIsotopeSpectrumMatch>(spectrum,
72 peptideSp,
73 parent_charge,
76 max_isotope_number,
77 1);
78
79 msp_peptideSpectrumMatch.get()->dropPeaksLackingMonoisotope();
80 m_spectrumSumIntensity = spectrum.sumY();
81
82
83 qDebug() << " accumulate";
84 std::vector<double> delta_list;
85
86
87 // TODO compute number of matched complementary peaks having m/z compatible
88 // with the precursor
89
90 m_precursorTheoreticalMz = peptideSp.get()->getMz(parent_charge);
91 m_precursorTheoreticalMass = peptideSp.get()->getMass();
92 m_parentCharge = parent_charge;
93
94
95 findComplementIonPairs(peptideSp);
96
97
98 for(const pappso::PeakIonIsotopeMatch &peak_ion :
99 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
100 {
101 delta_list.push_back(
102 peak_ion.getPeptideFragmentIonSp().get()->getMz(peak_ion.getCharge()) -
103 peak_ion.getPeak().x);
104 }
106 std::accumulate(delta_list.begin(), delta_list.end(), 0);
107
108 qDebug() << " delta_list.size()=" << delta_list.size();
112 if(delta_list.size() > 0)
113 {
114 m_matchedMzDiffMean = sum / ((pappso::pappso_double)delta_list.size());
115
116 std::sort(delta_list.begin(), delta_list.end());
117 m_matchedMzDiffMedian = delta_list[(delta_list.size() / 2)];
118
119
120 qDebug() << " sd";
122 for(pappso::pappso_double val : delta_list)
123 {
124 // sd = sd + ((val - mean) * (val - mean));
125 m_matchedMzDiffSd += std::pow((val - m_matchedMzDiffMean), 2);
126 }
127 m_matchedMzDiffSd = m_matchedMzDiffSd / delta_list.size();
129 }
130 else
131 {
132 }
133}
134
135
136double
138{
139 double sum = 0;
140 for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
141 {
142 if(peak_ion_match.getPeptideIonType() == ion_type)
143 {
144 sum += peak_ion_match.getPeak().y;
145 }
146 }
147 return sum;
148}
149double
151{
152 double sum = 0;
153 for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
154 {
155 sum += peak_ion_match.getPeak().y;
156 }
157 return sum;
158}
159
160double
165
166std::size_t
168{
169 return m_peakIonPairs.size();
170}
171
172const std::vector<
173 std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>> &
175{
176 return m_peakIonPairs;
177}
178
179double
181{
182
183 double sum = 0;
184 for(auto &peak_pairs : m_peakIonPairs)
185 {
186 sum += peak_pairs.first.getPeak().y;
187 sum += peak_pairs.second.getPeak().y;
188 }
189 return sum;
190}
191
192double
194{
195 return m_matchedMzDiffSd;
196}
197
198double
200{
201 return m_matchedMzDiffMean;
202}
203
204
205std::size_t
207{
208 return msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().size();
209}
210
211std::size_t
213{
214 std::size_t max = 0;
215 auto peak_ion_match_list =
216 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
217
218 peak_ion_match_list.erase(
219 std::remove_if(
220 peak_ion_match_list.begin(),
221 peak_ion_match_list.end(),
222 [ion_type](const PeakIonIsotopeMatch &a) {
223 if(a.getPeptideIonType() != ion_type)
224 return true;
225 if(a.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() > 0)
226 return true;
227 return false;
228 }),
229 peak_ion_match_list.end());
230
231 peak_ion_match_list.sort(
232 [](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
233 if(a.getCharge() < b.getCharge())
234 return true;
235 if(a.getPeptideIonType() < b.getPeptideIonType())
236 return true;
237 if(a.getPeptideFragmentIonSp().get()->size() <
238 b.getPeptideFragmentIonSp().get()->size())
239 return true;
240 return false;
241 });
242
243 unsigned int charge = 0;
244 std::size_t size = 0;
245 std::size_t count = 0;
246 for(std::list<PeakIonIsotopeMatch>::iterator it = peak_ion_match_list.begin();
247 it != peak_ion_match_list.end();
248 it++)
249 {
250 qDebug()
251 << it->toString() << max << " " << it->getPeak().x << " "
252 << it->getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber();
253 count++;
254 if((charge != it->getCharge()) ||
255 (size != (it->getPeptideFragmentIonSp().get()->size() - 1)))
256 {
257 count = 1;
258 charge = it->getCharge();
259 }
260 if(max < count)
261 max = count;
262
263 size = it->getPeptideFragmentIonSp().get()->size();
264 }
265
266 return max;
267}
268
269std::size_t
271{
272 std::vector<bool> covered;
273 covered.resize(msp_peptide.get()->size(), false);
274
275 for(auto &peak : msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
276 {
277 if(peak.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() ==
278 0)
279 { // only count with monotisotopic
280 if(peak.getPeptideIonType() == ion_type)
281 {
282 covered[peak.getPeptideFragmentIonSp().get()->size() - 1] = true;
283 }
284 }
285 }
286 return std::count(covered.begin(), covered.end(), true);
287}
288
289std::size_t
291{
292
293 std::vector<bool> covered;
294 covered.resize(msp_peptide.get()->size(), false);
295
296 for(auto &peak_pair : m_peakIonPairs)
297 {
298 std::size_t pos =
299 peak_pair.first.getPeptideFragmentIonSp().get()->size() - 1;
300 covered[pos] = true;
301 covered[pos + 1] = true;
302 }
303 return std::count(covered.begin(), covered.end(), true);
304}
305
306
307double
309 pappso::PeptideIon ion_type) const
310{
311 std::list<pappso::PeakIonIsotopeMatch> peak_ion_type =
312 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
313
314 peak_ion_type.remove_if([ion_type](const PeakIonIsotopeMatch &a) {
315 return (a.getPeptideIonType() != ion_type);
316 });
317 auto peak_it = std::max_element(
318 peak_ion_type.begin(),
319 peak_ion_type.end(),
320 [](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
321 return (a.getPeak().y < b.getPeak().y);
322 });
323
324 if(peak_it == peak_ion_type.end())
325 return 0;
326 return peak_it->getPeak().y;
327}
328
329double
331 const
332{
333 auto peak_it = std::max_element(
334 m_peakIonPairs.begin(),
335 m_peakIonPairs.end(),
336 [](const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
337 &a,
338 const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
339 &b) {
340 return ((a.first.getPeak().y + a.second.getPeak().y) <
341 (b.first.getPeak().y + b.second.getPeak().y));
342 });
343
344 if(peak_it == m_peakIonPairs.end())
345 return 0;
346
347 return getIonPairPrecursorMassDelta(*peak_it);
348}
349
350double
352 const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
353 &ion_pair) const
354{
355 qDebug() << m_precursorTheoreticalMz << " " << ion_pair.first.getPeak().x
356 << " " << ion_pair.second.getPeak().x << " "
357 << ion_pair.second.getCharge() << " " << ion_pair.first.getCharge()
358 << " " << m_parentCharge;
359 double diff =
360 (m_precursorTheoreticalMass + (MHPLUS * ion_pair.first.getCharge())) /
361 ion_pair.first.getCharge();
362
363
364 return (diff - (ion_pair.first.getPeak().x + ion_pair.second.getPeak().x -
365 ((MHPLUS * ion_pair.first.getCharge())) /
366 ion_pair.first.getCharge()));
367}
368
369
370void
372{
373 std::size_t peptide_size = peptideSp.get()->size();
374 std::vector<PeakIonIsotopeMatch> ion_isotope_list(
375 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().begin(),
376 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().end());
377 for(const pappso::PeakIonIsotopeMatch &peak_ion_ext : ion_isotope_list)
378 {
379 if(peptideIonIsNter(peak_ion_ext.getPeptideIonType()))
380 {
381 auto it = findComplementIonType(ion_isotope_list.begin(),
382 ion_isotope_list.end(),
383 peak_ion_ext,
384 peptide_size);
385 if(it != ion_isotope_list.end())
386 { // contains the complementary ion
387 m_peakIonPairs.push_back({peak_ion_ext, *it});
388 }
389 }
390 }
391}
392
393
396{
397 std::size_t peptide_size = msp_peptide.get()->size();
398 std::vector<PeakIonIsotopeMatch> ion_isotope_list(
399 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().begin(),
400 msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().end());
401
402 Trace scaterplot;
403
404 for(PeptideIon ion_type : m_ionList)
405 {
406 std::vector<double> mono_th_intensities(peptide_size, 0);
407 std::vector<double> isotope_th_intensities(peptide_size, 0);
408
409 std::vector<double> mono_exp_intensities(peptide_size, 0);
410 std::vector<double> isotope_exp_intensities(peptide_size, 0);
411 for(const PeakIonIsotopeMatch &peak_ion_match : ion_isotope_list)
412 {
413 if(peak_ion_match.getPeptideIonType() == ion_type)
414 {
415 std::size_t vector_position =
416 peak_ion_match.getPeptideFragmentIonSp().get()->size() - 1;
417 PeptideNaturalIsotopeAverageSp iso_average_sp =
418 peak_ion_match.getPeptideNaturalIsotopeAverageSp();
419 if(iso_average_sp.get()->getIsotopeNumber() == 0)
420 {
421 mono_th_intensities[vector_position] =
422 iso_average_sp.get()->getIntensityRatio();
423 mono_exp_intensities[vector_position] =
424 peak_ion_match.getPeak().y;
425 }
426 else if(iso_average_sp.get()->getIsotopeNumber() == 1)
427 {
428 isotope_th_intensities[vector_position] =
429 iso_average_sp.get()->getIntensityRatio();
430 isotope_exp_intensities[vector_position] =
431 peak_ion_match.getPeak().y;
432 }
433 }
434 }
435
436 for(std::size_t i = 0; i < mono_th_intensities.size(); i++)
437 {
438 if((mono_th_intensities[i] != 0) && (isotope_th_intensities[i] != 0))
439 {
440 DataPoint xy(mono_th_intensities[i] / isotope_th_intensities[i],
441 mono_exp_intensities[i] /
442 isotope_exp_intensities[i]);
443 scaterplot.push_back(xy);
444 }
445 }
446 }
447
448 scaterplot.sortX();
449
450 LinearRegression linear_regression(scaterplot);
451
452 return linear_regression;
453}
Class to represent a mass spectrum.
static PrecisionPtr getPrecisionPtrFractionInstance(PrecisionPtr origin, double fraction)
get the fraction of an existing precision pointer
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
std::shared_ptr< FilterChargeDeconvolution > msp_filterChargeDeconvolution
const std::vector< std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > > & getPeakIonPairs() const
std::shared_ptr< FilterResampleKeepGreater > msp_filterKeepGreater
std::shared_ptr< FilterMzExclusion > msp_filterMzExclusion
double getIonPairPrecursorMassDelta(const std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > &ion_pair) const
double m_precursorTheoreticalMz
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 getTotalIntensityOfMatchedIons() const
sum of matched peak intensities
LinearRegression getIonIsotopeLinearRegression() const
unsigned int m_parentCharge
double getMaxIntensityPeakIonMatch(PeptideIon ion_type) const
std::list< PeptideIon > m_ionList
std::shared_ptr< PeptideIsotopeSpectrumMatch > msp_peptideSpectrumMatch
std::size_t getMaxConsecutiveIon(PeptideIon ion_type)
get the maximum consecutive fragments of one ion type
void setPeptideSpectrumCharge(const pappso::PeptideSp peptideSp, const MassSpectrum *p_spectrum, unsigned int parent_charge, unsigned int max_isotope_number)
PrecisionPtr m_ms2precision
double getIntensityOfMatchedIon(PeptideIon ion_type)
get the sum of intensity of a specific ion
pappso::PeptideSp msp_peptide
double m_precursorTheoreticalMass
PsmFeatures(PrecisionPtr ms2precision, double minimumMz)
compute psm features
double getMatchedMzDiffSd() const
get standard deviation of matched peak mass delta
double m_spectrumSumIntensity
std::size_t getAaSequenceCoverage(PeptideIon ion_type)
number of amino acid covered by matched ions
A simple container of DataPoint instances.
Definition trace.h:148
void sortX(SortOrder sort_order=SortOrder::ascending)
Definition trace.cpp:1086
pappso_double sumY() const
Definition trace.cpp:1023
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:429
@ y
Cter amino ions.
@ b
Nter acylium ions.
@ xy
(x,y) format
std::shared_ptr< const Peptide > PeptideSp
const pappso_double MHPLUS(1.007276466879)
double pappso_double
A type definition for doubles.
Definition types.h:50
std::shared_ptr< const PeptideNaturalIsotopeAverage > PeptideNaturalIsotopeAverageSp
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