Loading...
Searching...
No Matches
spomsspectrum.cpp
Go to the documentation of this file.
/**
* \file pappsomspp/processing/specpeptidoms/spomsspectrum.cpp
* \date 24/03/2025
* \author Aurélien Berthier
* \brief SpecPeptidOMS Spectrum
*
* C++ implementation of the SpecPeptidOMS algorithm described in :
* (1) Benoist, É.; Jean, G.; Rogniaux, H.; Fertin, G.; Tessier, D. SpecPeptidOMS Directly and
* Rapidly Aligns Mass Spectra on Whole Proteomes and Identifies Peptides That Are Not Necessarily
* Tryptic: Implications for Peptidomics. J. Proteome Res. 2025.
* https://doi.org/10.1021/acs.jproteome.4c00870.
*/
/*
* Copyright (c) 2025 Aurélien Berthier
* <aurelien.berthier@ls2n.fr>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <algorithm>
#include <unordered_set>
#include "spomsspectrum.h"
#include "pappsomspp/core/processing/specpeptidoms/types.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "pappsomspp/core/amino_acid/aacode.h"
namespace pappso
{
namespace specpeptidoms
{
// SpOMSSpectrum::SpOMSSpectrum(const specglob::ExperimentalSpectrum &exp_spectrum)
pappso::specpeptidoms::SpOMSSpectrum::SpOMSSpectrum(pappso::QualifiedMassSpectrum &qmass_spectrum,
pappso::PrecisionPtr precision_ptr,
const pappso::AaCode &aaCode)
: std::vector<pappso::specglob::ExperimentalSpectrumDataPoint>(
specglob::ExperimentalSpectrum(qmass_spectrum, precision_ptr)),
m_qualifiedMassSpectrum(qmass_spectrum),
m_precision_ptr(precision_ptr),
m_aaCode(aaCode),
m_precursor_mass_error(0)
{
m_aapositions.reserve(m_aaCode.getSize());
for(std::size_t iter = 0; iter < m_aaCode.getSize(); iter++)
{
m_aapositions.push_back(std::make_shared<std::vector<AaPosition>>());
m_aapositions.back()->reserve(this->size() - 1);
}
m_supported_peaks.reserve(this->size());
m_supported_peaks.push_back(std::make_shared<std::vector<uint8_t>>());
m_reindexed_peaks.push_back(0);
for(std::size_t iter = 1; iter < this->size(); iter++)
{
m_supported_peaks.push_back(std::make_shared<std::vector<uint8_t>>());
m_reindexed_peaks.push_back(-1);
}
this->at(0).peak_mz = pappso::MHPLUS + 2 * pappso::MPROTIUM + pappso::MASSOXYGEN;
this->back().peak_mz = m_qualifiedMassSpectrum.getPrecursorMass() + pappso::MHPLUS;
preprocessSpectrum();
}
pappso::specpeptidoms::SpOMSSpectrum::SpOMSSpectrum(const SpOMSSpectrum &other)
: std::vector<pappso::specglob::ExperimentalSpectrumDataPoint>(
pappso::specglob::ExperimentalSpectrum(other.m_qualifiedMassSpectrum, other.m_precision_ptr)),
m_qualifiedMassSpectrum(other.m_qualifiedMassSpectrum),
m_aapositions(other.m_aapositions),
m_precision_ptr(other.m_precision_ptr),
m_supported_peaks(other.m_supported_peaks),
m_reindexed_peaks(other.m_reindexed_peaks),
m_aaCode(other.m_aaCode),
m_complementary_peak_indexes(other.m_complementary_peak_indexes),
m_precursor_mass_error(other.m_precursor_mass_error)
{
}
pappso::specpeptidoms::SpOMSSpectrum::SpOMSSpectrum(const SpOMSSpectrum &other,
double precursor_mass_error)
: std::vector<pappso::specglob::ExperimentalSpectrumDataPoint>(
pappso::specglob::ExperimentalSpectrum(
other.m_qualifiedMassSpectrum, other.m_precision_ptr, precursor_mass_error)),
m_qualifiedMassSpectrum(other.m_qualifiedMassSpectrum),
m_precision_ptr(other.m_precision_ptr),
m_aaCode(other.m_aaCode),
m_precursor_mass_error(precursor_mass_error)
{
m_aapositions.reserve(m_aaCode.getSize());
for(std::size_t iter = 0; iter < m_aaCode.getSize(); iter++)
{
m_aapositions.push_back(std::make_shared<std::vector<AaPosition>>());
m_aapositions.back()->reserve(this->size() - 1);
}
m_supported_peaks.reserve(this->size());
m_supported_peaks.push_back(std::make_shared<std::vector<uint8_t>>());
m_reindexed_peaks.push_back(0);
for(std::size_t iter = 1; iter < this->size(); iter++)
{
m_supported_peaks.push_back(std::make_shared<std::vector<uint8_t>>());
m_reindexed_peaks.push_back(-1);
}
this->at(0).peak_mz = pappso::MHPLUS + 2 * pappso::MPROTIUM + pappso::MASSOXYGEN;
this->back().peak_mz =
m_qualifiedMassSpectrum.getPrecursorMass() + pappso::MHPLUS + precursor_mass_error;
preprocessSpectrum();
}
pappso::specpeptidoms::SpOMSSpectrum::~SpOMSSpectrum()
{
}
// Add comments !!
void
pappso::specpeptidoms::SpOMSSpectrum::preprocessSpectrum()
{
// bool found;
uint8_t aa;
std::vector<double>::iterator iter1, iter2;
std::size_t peak1, peak2, next_l_peak;
std::vector<double> mass_list = getMassList();
peak1 = -1;
for(iter1 = mass_list.begin(); iter1 != mass_list.end(); iter1++)
{
peak1++;
peak2 = peak1;
for(iter2 = iter1 + 1; iter2 != mass_list.end(); iter2++)
{
peak2++;
aa = m_aaCode.getAaCodeByMass(*(iter2) - *(iter1), m_precision_ptr);
if(aa != 0)
{
next_l_peak = 0;
for(std::size_t iter = 1; iter < peak1;
iter++) // Search of the closer supported left peak.
// Possible optimization => search from the right
{
if(m_reindexed_peaks.at(iter) >= 0)
{
next_l_peak = iter;
}
}
if(m_reindexed_peaks.at(peak2) == -1)
{
addSupportedPeak(peak2);
m_supported_peaks.at(peak2)->push_back(aa);
}
else
{
auto it = std::find(
m_supported_peaks.at(peak2)->begin(), m_supported_peaks.at(peak2)->end(), aa);
if(it == m_supported_peaks.at(peak2)->end())
{
m_supported_peaks.at(peak2)->push_back(aa);
}
}
if(m_reindexed_peaks.at(peak1) >= 0)
{
addAaPosition(aa, peak2, peak1, *(iter1), true);
}
else
{
addAaPosition(aa, peak2, next_l_peak, *(iter1), false);
}
}
}
}
removeUnsupportedMasses();
correctPeakIndexes();
for(uint8_t aa = 1; aa < m_aaCode.getSize() + 1; ++aa)
{
qDebug() << m_aaCode.getAa(aa).getLetter();
for(auto iter = m_aapositions.at(aa - 1)->begin(); iter != m_aapositions.at(aa - 1)->end();
++iter)
{
qDebug() << iter->l_peak << this->at(iter->l_peak).peak_mz << iter->r_peak
<< this->at(iter->r_peak).peak_mz << iter->l_support << iter->condition;
}
}
// std::size_t i = 0;
// for(auto &data_point : *this)
// {
// data_point.indice = i;
// i++;
// }
fillComplementaryPeakIndexes();
}
// pappso::Aa const *
// SpOMSSpectrum::findAAMass(double mass, bool *found) const
// {
// bool ok;
// // auto charge = m_qualifiedMassSpectrum.getPrecursorCharge(&ok);
// if(!ok)
// {
// throw pappso::PappsoException(
// QObject::tr("precursor charge is not defined in spectrum %1")
// .arg(m_qualifiedMassSpectrum.getMassSpectrumId().getNativeId()));
// }
// pappso::MzRange mz_range(mass / m_qualifiedMassSpectrum.getPrecursorCharge(),
// m_precision_ptr);
// for(std::unordered_map<const Aa, double>::const_iterator aa = aaMasses.begin();
// aa != aaMasses.end();
// aa++)
// {
// if(mz_range.contains(aa->second))
// {
// if(found != nullptr)
// {
// *found = true;
// }
// return &(aa->first);
// }
// }
// if(found != nullptr)
// {
// *found = false;
// }
// return nullptr;
// }
// Not sure if optimal
void
pappso::specpeptidoms::SpOMSSpectrum::removeUnsupportedMasses()
{
std::vector<specglob::ExperimentalSpectrumDataPoint> kept_peaks;
for(std::vector<specglob::ExperimentalSpectrumDataPoint>::iterator iter = this->begin();
iter != this->end();
iter++)
{
if(m_reindexed_peaks.at(iter->indice) >= 0)
{
kept_peaks.push_back(*iter);
}
}
this->clear();
this->assign(kept_peaks.begin(), kept_peaks.end());
}
void
pappso::specpeptidoms::SpOMSSpectrum::addAaPosition(
uint8_t aa, const std::size_t r_peak, const std::size_t l_peak, double l_mass, bool l_support)
{
// aa=0 corresponds to no amino acid identified, thus aa is always >=1. We substract 1 to aa to
// avoid keeping an empty, useless vector.
m_aapositions.at(aa - 1)->push_back(
{r_peak, l_peak, l_mass, computeCondition(l_peak, l_support), l_support});
qDebug() << "l_peak" << l_peak << "r_peak" << r_peak << "l_mass" << l_mass << "l_support"
<< l_support << "condition" << computeCondition(l_peak, l_support);
}
uint32_t
pappso::specpeptidoms::SpOMSSpectrum::computeCondition(const std::size_t l_peak,
bool l_support) const
{
uint32_t condition;
if(!l_support)
{
condition = 1;
}
else if(l_peak == 0)
{
condition = 2;
}
else
{
condition = 0;
for(std::vector<uint8_t>::iterator aa = m_supported_peaks.at(l_peak)->begin();
aa != m_supported_peaks.at(l_peak)->end();
aa++)
{
condition += 2 << *(aa);
}
}
return condition;
}
const std::vector<pappso::specpeptidoms::AaPosition> &
pappso::specpeptidoms::SpOMSSpectrum::getAaPositions(std::uint8_t aa_code) const
{
return *m_aapositions.at(aa_code - 1);
}
std::vector<pappso::specpeptidoms::AaPosition>
pappso::specpeptidoms::SpOMSSpectrum::getAaPositions(
std::uint8_t aa_code, std::vector<std::size_t> &peaks_to_remove) const
{
std::vector<AaPosition> aa_positions;
for(auto aap : *m_aapositions.at(aa_code - 1))
{
if(std::find(peaks_to_remove.begin(), peaks_to_remove.end(), aap.r_peak) ==
peaks_to_remove.end())
{
aa_positions.push_back(aap);
}
}
return aa_positions;
}
std::vector<double>
pappso::specpeptidoms::SpOMSSpectrum::getMassList() const
{
std::vector<double> mass_list;
for(const specglob::ExperimentalSpectrumDataPoint &n : *this)
{
mass_list.push_back(n.peak_mz);
}
return mass_list;
}
pappso::specglob::ExperimentalSpectrumDataPointType
pappso::specpeptidoms::SpOMSSpectrum::peakType(std::size_t indice) const
{
return this->at(indice).type;
}
uint
pappso::specpeptidoms::SpOMSSpectrum::getPrecursorCharge() const
{
return m_qualifiedMassSpectrum.getPrecursorCharge();
}
double
pappso::specpeptidoms::SpOMSSpectrum::getPrecursorMass() const
{
return m_qualifiedMassSpectrum.getPrecursorMass();
}
double
pappso::specpeptidoms::SpOMSSpectrum::getMZShift(std::size_t l_peak, std::size_t r_peak) const
{
if(std::max(r_peak, l_peak) > size())
{
throw pappso::ExceptionOutOfRange(
QObject::tr("getMZShift : l_peak %1 or r_peak %2 greater than size %3")
.arg(l_peak)
.arg(r_peak)
.arg(size()));
}
return this->at(r_peak).peak_mz - this->at(l_peak).peak_mz;
}
double
pappso::specpeptidoms::SpOMSSpectrum::getMissingMass(std::size_t peak) const
{
if(peak > size())
{
throw pappso::ExceptionOutOfRange(
QObject::tr("getMissingMass : peak %1 greater than size %2").arg(peak).arg(size()));
}
return this->m_qualifiedMassSpectrum.getPrecursorMass() - m_precursor_mass_error -
this->at(peak).peak_mz + MHPLUS;
}
void
pappso::specpeptidoms::SpOMSSpectrum::addSupportedPeak(std::size_t peak)
{
std::size_t counter = 0;
for(std::size_t iter = 0; iter < peak; iter++)
{
if(m_reindexed_peaks.at(iter) >= 0)
{
counter++;
}
}
m_reindexed_peaks.at(peak) = counter;
for(std::size_t iter = peak + 1; iter < m_reindexed_peaks.size(); iter++)
{
if(m_reindexed_peaks.at(iter) >= 0)
{
m_reindexed_peaks.at(iter)++;
}
}
}
void
pappso::specpeptidoms::SpOMSSpectrum::correctPeakIndexes()
{
for(auto aa = m_aapositions.begin(); aa != m_aapositions.end(); aa++)
{
for(auto aap = aa->get()->begin(); aap != aa->get()->end(); aap++)
{
aap->l_peak = m_reindexed_peaks.at(aap->l_peak);
aap->r_peak = m_reindexed_peaks.at(aap->r_peak);
}
}
}
void
pappso::specpeptidoms::SpOMSSpectrum::fillComplementaryPeakIndexes()
{
std::size_t left_index, right_index;
m_complementary_peak_indexes.reserve(this->size());
while(m_complementary_peak_indexes.size() < this->size())
{
m_complementary_peak_indexes.push_back(0);
}
left_index = 0;
right_index = this->size() - 1;
double comp_mass = m_qualifiedMassSpectrum.getPrecursorMass() + 2 * MHPLUS;
while(left_index < right_index)
{
pappso::MzRange mz_range(comp_mass - this->at(left_index).peak_mz, m_precision_ptr);
if(mz_range.contains(this->at(right_index).peak_mz))
{
m_complementary_peak_indexes.at(left_index) = right_index;
m_complementary_peak_indexes.at(right_index) = left_index;
qDebug() << left_index << right_index;
}
if(comp_mass - this->at(left_index).peak_mz - this->at(right_index).peak_mz >= 0)
{
left_index++;
}
else
{
right_index--;
}
}
}
std::size_t
pappso::specpeptidoms::SpOMSSpectrum::getComplementaryPeak(std::size_t peak) const
{
return m_complementary_peak_indexes.at(peak);
}
} // namespace specpeptidoms
} // namespace pappso
aacode.h
pappso::AaCode
collection of integer code for each amino acid 0 => null 1 to 20 => amino acid sorted by there mass (...
Definition aacode.h:44
pappso::ExceptionOutOfRange
Definition exceptionoutofrange.h:32
pappso::MzRange
Definition mzrange.h:46
pappso::MzRange::contains
bool contains(pappso_double) const
Definition mzrange.cpp:115
pappso::QualifiedMassSpectrum
Class representing a fully specified mass spectrum.
Definition qualifiedmassspectrum.h:99
pappso::specpeptidoms::SpOMSSpectrum::preprocessSpectrum
void preprocessSpectrum()
Preprocess the spectrum.
Definition spomsspectrum.cpp:124
pappso::specpeptidoms::SpOMSSpectrum::getMZShift
double getMZShift(std::size_t l_peak, std::size_t r_peak) const
Returns the mz difference between two peaks.
Definition spomsspectrum.cpp:348
pappso::specpeptidoms::SpOMSSpectrum::getPrecursorCharge
uint getPrecursorCharge() const
Returns the spectrum's precursor's charge.
Definition spomsspectrum.cpp:336
pappso::specpeptidoms::SpOMSSpectrum::SpOMSSpectrum
SpOMSSpectrum(pappso::QualifiedMassSpectrum &qmass_spectrum, pappso::PrecisionPtr precision_ptr, const pappso::AaCode &aaCode)
Definition spomsspectrum.cpp:45
pappso::specpeptidoms::SpOMSSpectrum::getMissingMass
double getMissingMass(std::size_t peak) const
Returns the missing mass between a peak and the precursor's mass (shift at the end).
Definition spomsspectrum.cpp:362
pappso::specpeptidoms::SpOMSSpectrum::m_complementary_peak_indexes
std::vector< std::size_t > m_complementary_peak_indexes
Definition spomsspectrum.h:166
pappso::specpeptidoms::SpOMSSpectrum::m_supported_peaks
std::vector< std::shared_ptr< std::vector< uint8_t > > > m_supported_peaks
Definition spomsspectrum.h:161
pappso::specpeptidoms::SpOMSSpectrum::computeCondition
uint32_t computeCondition(const std::size_t l_peak, bool l_support) const
Computes the "condition" integer, used to apply the three peaks rule.
Definition spomsspectrum.cpp:269
pappso::specpeptidoms::SpOMSSpectrum::removeUnsupportedMasses
void removeUnsupportedMasses()
Removes the unsupported peaks (without an amino acid to the left) from the spectrum.
Definition spomsspectrum.cpp:240
pappso::specpeptidoms::SpOMSSpectrum::m_precursor_mass_error
double m_precursor_mass_error
Definition spomsspectrum.h:167
pappso::specpeptidoms::SpOMSSpectrum::addAaPosition
void addAaPosition(uint8_t aa, const std::size_t r_peak, const std::size_t l_peak, double l_mass, bool l_support)
Adds an amino acid position to the data structure.
Definition spomsspectrum.cpp:257
pappso::specpeptidoms::SpOMSSpectrum::m_precision_ptr
pappso::PrecisionPtr m_precision_ptr
Definition spomsspectrum.h:159
pappso::specpeptidoms::SpOMSSpectrum::m_qualifiedMassSpectrum
pappso::QualifiedMassSpectrum m_qualifiedMassSpectrum
Definition spomsspectrum.h:157
pappso::specpeptidoms::SpOMSSpectrum::m_aapositions
std::vector< std::shared_ptr< std::vector< AaPosition > > > m_aapositions
Definition spomsspectrum.h:158
pappso::specpeptidoms::SpOMSSpectrum::correctPeakIndexes
void correctPeakIndexes()
Reindexes the peaks after removal of the unsupported peaks.
Definition spomsspectrum.cpp:395
pappso::specpeptidoms::SpOMSSpectrum::addSupportedPeak
void addSupportedPeak(std::size_t peak)
Add a peak to the supported peaks list.
Definition spomsspectrum.cpp:374
pappso::specpeptidoms::SpOMSSpectrum::fillComplementaryPeakIndexes
void fillComplementaryPeakIndexes()
For each point of the spectrum, indicate the index of its complementary peak;.
Definition spomsspectrum.cpp:408
pappso::specpeptidoms::SpOMSSpectrum::m_reindexed_peaks
std::vector< int > m_reindexed_peaks
Definition spomsspectrum.h:163
pappso::specpeptidoms::SpOMSSpectrum::getComplementaryPeak
std::size_t getComplementaryPeak(std::size_t peak) const
Definition spomsspectrum.cpp:442
pappso::specpeptidoms::SpOMSSpectrum::getAaPositions
const std::vector< AaPosition > & getAaPositions(std::uint8_t aa_code) const
Returns the list of aa_positions for a given amino acid code.
Definition spomsspectrum.cpp:296
pappso::specpeptidoms::SpOMSSpectrum::peakType
specglob::ExperimentalSpectrumDataPointType peakType(std::size_t indice) const
Returns the type of one of the spectrum's peaks.
Definition spomsspectrum.cpp:330
pappso::specpeptidoms::SpOMSSpectrum::getMassList
std::vector< double > getMassList() const
Returns the spectrum's list of masses.
Definition spomsspectrum.cpp:319
pappso::specpeptidoms::SpOMSSpectrum::~SpOMSSpectrum
virtual ~SpOMSSpectrum()
Definition spomsspectrum.cpp:118
pappso::specpeptidoms::SpOMSSpectrum::getPrecursorMass
double getPrecursorMass() const
Definition spomsspectrum.cpp:341
pappso::specpeptidoms::SpOMSSpectrum::m_aaCode
const pappso::AaCode & m_aaCode
Definition spomsspectrum.h:165
exceptionoutofrange.h
pappso::specglob
Definition experimentalspectrum.cpp:44
pappso::specglob::ExperimentalSpectrumDataPointType
ExperimentalSpectrumDataPointType
Definition types.h:78
pappso::specpeptidoms
Definition correctiontree.h:40
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappso::MHPLUS
const pappso_double MHPLUS(1.007276466879)
pappso::MPROTIUM
const pappso_double MPROTIUM(1.007825032241)
pappso::uint
unsigned int uint
Definition types.h:67
pappso::MASSOXYGEN
const pappso_double MASSOXYGEN(15.99491461956)
pappso::PrecisionPtr
const PrecisionBase * PrecisionPtr
Definition precision.h:122
pappsoexception.h
types.h
spomsspectrum.h
pappso::specglob::ExperimentalSpectrumDataPoint
Definition experimentalspectrum.h:51
pappso::specglob::ExperimentalSpectrumDataPoint::peak_mz
double peak_mz
Definition experimentalspectrum.h:53