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
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
54{
55}
56
57void
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,
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();
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";
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
135double
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}
148double
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
159double
161{
163}
164
165std::size_t
167{
168 return m_peakIonPairs.size();
169}
170
171const std::vector<
172 std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>> &
174{
175 return m_peakIonPairs;
176}
177
178double
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
191double
193{
194 return m_matchedMzDiffSd;
195}
196
197double
199{
200 return m_matchedMzDiffMean;
201}
202
203
204std::size_t
206{
207 return msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().size();
208}
209
210std::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
268std::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
284std::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
302double
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
324double
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
345double
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
365void
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
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:886
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:385
@ 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:48
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