libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
experimentalspectrum.cpp
Go to the documentation of this file.
1/**
2 * \file pappsomspp/processing/specglob/experimentalspectrum.cpp
3 * \date 08/11/2023
4 * \author Olivier Langella
5 * \brief transform a spectrum to SpecGlob spectra
6 *
7 * C++ implementation of the SpecGlob algorithm described in :
8 * 1. Prunier, G. et al. Fast alignment of mass spectra in large proteomics
9 * datasets, capturing dissimilarities arising from multiple complex
10 * modifications of peptides. BMC Bioinformatics 24, 421 (2023).
11 *
12 * HAL Id : hal-04296170 , version 1
13 * Mot de passe : hxo20cl
14 * DOI : 10.1186/s12859-023-05555-y
15 */
16
17/*
18 * SpecGlobTool, Spectra to peptide alignment tool
19 * Copyright (C) 2023 Olivier Langella
20 * <olivier.langella@universite-paris-saclay.fr>
21 *
22 * This program is free software: you can redistribute it and/or modify
23 * it under the terms of the GNU General Public License as published by
24 * the Free Software Foundation, either version 3 of the License, or
25 * (at your option) any later version.
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU General Public License for more details.
31 *
32 * You should have received a copy of the GNU General Public License
33 * along with this program. If not, see <http://www.gnu.org/licenses/>.
34 *
35 */
37#include <QDebug>
38#include "../../pappsoexception.h"
39#include <QObject>
40
41namespace pappso
42{
43namespace specglob
44{
46 const pappso::QualifiedMassSpectrum &qmass_spectrum,
47 pappso::PrecisionPtr precision_ptr)
48 : std::vector<ExperimentalSpectrumDataPoint>()
49{
50 m_qualifiedMassSpectrum = qmass_spectrum;
51 m_precisionPtr = precision_ptr;
53}
54
62
66
67void
69{
70
71 bool ok;
73 if(!ok)
74 {
76 QObject::tr("precursor charge is not defined in spectrum %1")
78 }
79
80 double mz_prec = m_qualifiedMassSpectrum.getPrecursorMz(&ok);
81 if(!ok)
82 {
84 QObject::tr("precursor m/z is not defined in spectrum %1")
86 }
87
88 // compute precursor mass given the charge state
89 m_precursorMass = mz_prec * (double)charge;
90 m_precursorMass -= pappso::MHPLUS * (double)charge;
92
93 std::vector<double> mz_list =
95
96 std::size_t mz_current_indice = 0;
97 auto itend = m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->end();
98 for(std::vector<pappso::DataPoint>::const_iterator it =
100 it != itend;
101 it++)
102 {
103 double current_mz = it->x;
104
105 double symetric_current_mz = getSymetricMz(current_mz);
106
107 auto itpair_symetric = findMz(symetric_current_mz);
108 if(itpair_symetric != itend)
109 {
110 // there is a counterpart
111 push_back({ExperimentalSpectrumDataPointType::both, it->x, 0});
112 qDebug() << "current_mz=" << current_mz << " both";
113 }
114 else
115 {
116 // there is no counterpart
117 push_back({ExperimentalSpectrumDataPointType::native, it->x, 0});
119 symetric_current_mz,
120 0});
121 qDebug() << "current_mz=" << current_mz << " symetrics";
122 }
123
124 mz_current_indice++;
125 }
126
127 // we add a peak with NT mass (1.0078) if it is not detected, to give better
128 // chance to align the first amino acid in b if it is present
129 if(findMz(pappso::MPROTIUM) == itend)
130 {
131 push_back(
133 }
134
135 // We add the B peak corresponding to the precursor (complete peptide)
136 double precusorBion = m_targetMzSum - pappso::MASSH2O - pappso::MPROTIUM;
137
138 if(findMz(precusorBion) == itend)
139 {
140 push_back(
142 }
143
144 std::sort(begin(),
145 end(),
148 return (a.symetric_mz < b.symetric_mz);
149 });
150
151 std::size_t i = 0;
152 for(auto &data_point : *this)
153 {
154 data_point.indice = i;
155 i++;
156 }
157}
158
159std::vector<pappso::DataPoint>::const_iterator
161{
163 auto itend = m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->end();
164 auto it = findFirstEqualOrGreaterX(
167 mz_range.lower());
168
169 if(it != itend)
170 {
171 if(it->x <= mz_range.upper())
172 {
173 return it;
174 }
175 }
176 return itend;
177}
178
179double
184
185std::vector<double>
187{
188 std::vector<double> mass_list;
189 for(const ExperimentalSpectrumDataPoint &n : *this)
190 {
191 mass_list.push_back(n.symetric_mz);
192 }
193
194 return mass_list;
195}
196
197double
199{
200 return m_targetMzSum - mz;
201}
202
203std::vector<double>
205{
206 std::vector<double> mass_list;
207 for(const ExperimentalSpectrumDataPoint &n : *this)
208 {
209 if(n.type == type)
210 mass_list.push_back(n.symetric_mz);
211 }
212
213 return mass_list;
214}
215
216QString
218{
219 QStringList all_element;
220 for(const ExperimentalSpectrumDataPoint &n : *this)
221 {
222 all_element
223 << QString("%1 %2").arg(n.symetric_mz).arg((std::uint8_t)n.type);
224 }
225 return QString("[%1]").arg(all_element.join("] ["));
226}
227
228double
233
234std::vector<ExperimentalSpectrumDataPoint>::const_reverse_iterator
236 std::size_t start_position, const pappso::MzRange &aaTheoMzRange) const
237{
238
239 qDebug() << "start_position" << start_position
240 << " lookfor=" << aaTheoMzRange.getMz();
241 std::vector<ExperimentalSpectrumDataPoint>::const_reverse_iterator itrbegin =
242 rbegin() + (size() - 1 - start_position);
243
244 qDebug() << itrbegin->indice << "mz=" << itrbegin->symetric_mz;
245 double eperimentalMzReference = itrbegin->symetric_mz;
246 auto itrend = this->rend();
247 // We check all row from j to 0
248 qDebug();
249 for(auto itr = itrbegin + 1; itr != itrend; ++itr)
250 {
251 qDebug() << itr->indice;
252 double experimentalMzDifference =
253 eperimentalMzReference - itr->symetric_mz;
254
255 if(experimentalMzDifference > aaTheoMzRange.upper())
256 {
257 // if we pass the mass of the theoretical amino acid, we stop to not
258 // over
259 // calculate
260 qDebug() << experimentalMzDifference << ">" << aaTheoMzRange.upper();
261 return itrend;
262 }
263 else if(experimentalMzDifference < aaTheoMzRange.lower())
264 {
265 continue;
266 }
267
268 qDebug() << itr->indice << " diff=" << experimentalMzDifference;
269 // if we found that j-k give an amino acid, we keep k value
270 return itr;
271 }
272 qDebug() << "rend";
273 return itrend; // a value of -1 to show that
274 // there is no k value
275}
276
282} // namespace specglob
283} // namespace pappso
const QString & getNativeId() const
pappso_double getMz() const
Definition mzrange.cpp:114
pappso_double lower() const
Definition mzrange.h:71
pappso_double upper() const
Definition mzrange.h:77
Class representing a fully specified mass spectrum.
MassSpectrumCstSPtr getMassSpectrumCstSPtr() const
Get the MassSpectrumCstSPtr.
uint getPrecursorCharge(bool *ok=nullptr) const
get precursor charge
const MassSpectrumId & getMassSpectrumId() const
Get the MassSpectrumId.
pappso_double getPrecursorMz(bool *ok=nullptr) const
get precursor mz
double getSymetricMz(double mz) const
compute the symetric mass for debuggin purpose
const pappso::QualifiedMassSpectrum & getQualifiedMassSpectrum() const
std::vector< ExperimentalSpectrumDataPoint >::const_reverse_iterator reverseFindDiffMz(std::size_t start_position, const pappso::MzRange &targeted_mass_range) const
find the peak for wich mass difference from rbegin corresponds to aaTheoMass Find if a peak back in t...
ExperimentalSpectrum(const pappso::QualifiedMassSpectrum &qmass_spectrum, pappso::PrecisionPtr precision_ptr)
void createSymetricPeakList()
add symetric peaks to the spectrum Create a SymetricPeakList that contain symmetric peaks and the inf...
std::vector< pappso::DataPoint >::const_iterator findMz(double mz)
find the correspondin mz in the mass spectrum (given the precision)
pappso::QualifiedMassSpectrum m_qualifiedMassSpectrum
transform a spectrum to SpecGlob spectra
ExperimentalSpectrumDataPointType
Definition types.h:78
@ symetric
new peak : computed symetric mass from a corresponding native peak
@ synthetic
does not correspond to existing peak, for computational purpose
@ both
both, the ion and the complement exists in the original spectrum
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
std::vector< DataPoint >::iterator findFirstEqualOrGreaterX(std::vector< DataPoint >::iterator begin, std::vector< DataPoint >::iterator end, const double &value)
find the first element in which X is equal or greater than the value searched important : it implies ...
Definition trace.cpp:71
const pappso_double MHPLUS(1.007276466879)
const pappso_double MPROTIUM(1.007825032241)
const pappso_double MASSH2O((MPROTIUM *2)+MASSOXYGEN)