Loading...
Searching...
No Matches
psmfeatures.cpp
Go to the documentation of this file.
/**
* \file pappsomspp/psm/features/psmfeatures.cpp
* \date 19/07/2022
* \author Olivier Langella
* \brief comutes various PSM (Peptide Spectrum Match) features
*/
/*******************************************************************************
* Copyright (c) 2022 Olivier Langella <Olivier.Langella@u-psud.fr>.
*
* This file is part of PAPPSOms-tools.
*
* PAPPSOms-tools 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.
*
* PAPPSOms-tools 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 PAPPSOms-tools. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "psmfeatures.h"
#include <memory>
#include <cmath>
using namespace pappso;
PsmFeatures::PsmFeatures(PrecisionPtr ms2precision, double minimumMz)
{
m_ms2precision = ms2precision;
m_ionList.push_back(PeptideIon::y);
m_ionList.push_back(PeptideIon::b);
msp_filterKeepGreater =
std::make_shared<FilterResampleKeepGreater>(minimumMz);
msp_filterChargeDeconvolution =
std::make_shared<FilterChargeDeconvolution>(m_ms2precision);
msp_filterMzExclusion = std::make_shared<FilterMzExclusion>(
PrecisionFactory::getPrecisionPtrFractionInstance(m_ms2precision, 0.5));
}
PsmFeatures::~PsmFeatures()
{
}
void
PsmFeatures::setPeptideSpectrumCharge(const pappso::PeptideSp peptideSp,
const MassSpectrum *p_spectrum,
unsigned int parent_charge,
unsigned int max_isotope_number)
{
m_peakIonPairs.clear();
msp_peptide = peptideSp;
MassSpectrum spectrum(*p_spectrum);
msp_filterKeepGreater.get()->filter(spectrum);
// msp_filterChargeDeconvolution.get()->filter(spectrum);
// msp_filterMzExclusion.get()->filter(spectrum);
msp_peptideSpectrumMatch =
std::make_shared<PeptideIsotopeSpectrumMatch>(spectrum,
peptideSp,
parent_charge,
m_ms2precision,
m_ionList,
max_isotope_number,
1);
msp_peptideSpectrumMatch.get()->dropPeaksLackingMonoisotope();
m_spectrumSumIntensity = spectrum.sumY();
qDebug() << " accumulate";
std::vector<double> delta_list;
// TODO compute number of matched complementary peaks having m/z compatible
// with the precursor
m_precursorTheoreticalMz = peptideSp.get()->getMz(parent_charge);
m_precursorTheoreticalMass = peptideSp.get()->getMass();
m_parentCharge = parent_charge;
findComplementIonPairs(peptideSp);
for(const pappso::PeakIonIsotopeMatch &peak_ion :
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
{
delta_list.push_back(
peak_ion.getPeptideFragmentIonSp().get()->getMz(peak_ion.getCharge()) -
peak_ion.getPeak().x);
}
pappso::pappso_double sum =
std::accumulate(delta_list.begin(), delta_list.end(), 0);
qDebug() << " delta_list.size()=" << delta_list.size();
m_matchedMzDiffMean = 0;
m_matchedMzDiffMedian = 0;
m_matchedMzDiffSd = 0;
if(delta_list.size() > 0)
{
m_matchedMzDiffMean = sum / ((pappso::pappso_double)delta_list.size());
std::sort(delta_list.begin(), delta_list.end());
m_matchedMzDiffMedian = delta_list[(delta_list.size() / 2)];
qDebug() << " sd";
m_matchedMzDiffSd = 0;
for(pappso::pappso_double val : delta_list)
{
// sd = sd + ((val - mean) * (val - mean));
m_matchedMzDiffSd += std::pow((val - m_matchedMzDiffMean), 2);
}
m_matchedMzDiffSd = m_matchedMzDiffSd / delta_list.size();
m_matchedMzDiffSd = std::sqrt(m_matchedMzDiffSd);
}
else
{
}
}
double
PsmFeatures::getIntensityOfMatchedIon(PeptideIon ion_type)
{
double sum = 0;
for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
{
if(peak_ion_match.getPeptideIonType() == ion_type)
{
sum += peak_ion_match.getPeak().y;
}
}
return sum;
}
double
PsmFeatures::getTotalIntensityOfMatchedIons() const
{
double sum = 0;
for(const PeakIonMatch &peak_ion_match : *msp_peptideSpectrumMatch.get())
{
sum += peak_ion_match.getPeak().y;
}
return sum;
}
double
PsmFeatures::getTotalIntensity() const
{
return m_spectrumSumIntensity;
}
std::size_t
pappso::PsmFeatures::countMatchedIonComplementPairs() const
{
return m_peakIonPairs.size();
}
const std::vector<
std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>> &
pappso::PsmFeatures::getPeakIonPairs() const
{
return m_peakIonPairs;
}
double
pappso::PsmFeatures::getTotalIntensityOfMatchedIonComplementPairs() const
{
double sum = 0;
for(auto &peak_pairs : m_peakIonPairs)
{
sum += peak_pairs.first.getPeak().y;
sum += peak_pairs.second.getPeak().y;
}
return sum;
}
double
pappso::PsmFeatures::getMatchedMzDiffSd() const
{
return m_matchedMzDiffSd;
}
double
pappso::PsmFeatures::getMatchedMzDiffMean() const
{
return m_matchedMzDiffMean;
}
std::size_t
pappso::PsmFeatures::getNumberOfMatchedIons() const
{
return msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().size();
}
std::size_t
pappso::PsmFeatures::getMaxConsecutiveIon(pappso::PeptideIon ion_type)
{
std::size_t max = 0;
auto peak_ion_match_list =
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
peak_ion_match_list.erase(
std::remove_if(
peak_ion_match_list.begin(),
peak_ion_match_list.end(),
[ion_type](const PeakIonIsotopeMatch &a) {
if(a.getPeptideIonType() != ion_type)
return true;
if(a.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() > 0)
return true;
return false;
}),
peak_ion_match_list.end());
peak_ion_match_list.sort(
[](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
if(a.getCharge() < b.getCharge())
return true;
if(a.getPeptideIonType() < b.getPeptideIonType())
return true;
if(a.getPeptideFragmentIonSp().get()->size() <
b.getPeptideFragmentIonSp().get()->size())
return true;
return false;
});
unsigned int charge = 0;
std::size_t size = 0;
std::size_t count = 0;
for(std::list<PeakIonIsotopeMatch>::iterator it = peak_ion_match_list.begin();
it != peak_ion_match_list.end();
it++)
{
qDebug()
<< it->toString() << max << " " << it->getPeak().x << " "
<< it->getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber();
count++;
if((charge != it->getCharge()) ||
(size != (it->getPeptideFragmentIonSp().get()->size() - 1)))
{
count = 1;
charge = it->getCharge();
}
if(max < count)
max = count;
size = it->getPeptideFragmentIonSp().get()->size();
}
return max;
}
std::size_t
pappso::PsmFeatures::getAaSequenceCoverage(pappso::PeptideIon ion_type)
{
std::vector<bool> covered;
covered.resize(msp_peptide.get()->size(), false);
for(auto &peak : msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList())
{
if(peak.getPeptideNaturalIsotopeAverageSp().get()->getIsotopeNumber() ==
0)
{ // only count with monotisotopic
if(peak.getPeptideIonType() == ion_type)
{
covered[peak.getPeptideFragmentIonSp().get()->size() - 1] = true;
}
}
}
return std::count(covered.begin(), covered.end(), true);
}
std::size_t
pappso::PsmFeatures::getComplementPairsAaSequenceCoverage()
{
std::vector<bool> covered;
covered.resize(msp_peptide.get()->size(), false);
for(auto &peak_pair : m_peakIonPairs)
{
std::size_t pos =
peak_pair.first.getPeptideFragmentIonSp().get()->size() - 1;
covered[pos] = true;
covered[pos + 1] = true;
}
return std::count(covered.begin(), covered.end(), true);
}
double
pappso::PsmFeatures::getMaxIntensityPeakIonMatch(
pappso::PeptideIon ion_type) const
{
std::list<pappso::PeakIonIsotopeMatch> peak_ion_type =
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList();
peak_ion_type.remove_if([ion_type](const PeakIonIsotopeMatch &a) {
return (a.getPeptideIonType() != ion_type);
});
auto peak_it = std::max_element(
peak_ion_type.begin(),
peak_ion_type.end(),
[](const PeakIonIsotopeMatch &a, const PeakIonIsotopeMatch &b) {
return (a.getPeak().y < b.getPeak().y);
});
if(peak_it == peak_ion_type.end())
return 0;
return peak_it->getPeak().y;
}
double
pappso::PsmFeatures::getMaxIntensityMatchedIonComplementPairPrecursorMassDelta()
const
{
auto peak_it = std::max_element(
m_peakIonPairs.begin(),
m_peakIonPairs.end(),
[](const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
&a,
const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
&b) {
return ((a.first.getPeak().y + a.second.getPeak().y) <
(b.first.getPeak().y + b.second.getPeak().y));
});
if(peak_it == m_peakIonPairs.end())
return 0;
return getIonPairPrecursorMassDelta(*peak_it);
}
double
pappso::PsmFeatures::getIonPairPrecursorMassDelta(
const std::pair<pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch>
&ion_pair) const
{
qDebug() << m_precursorTheoreticalMz << " " << ion_pair.first.getPeak().x
<< " " << ion_pair.second.getPeak().x << " "
<< ion_pair.second.getCharge() << " " << ion_pair.first.getCharge()
<< " " << m_parentCharge;
double diff =
(m_precursorTheoreticalMass + (MHPLUS * ion_pair.first.getCharge())) /
ion_pair.first.getCharge();
return (diff - (ion_pair.first.getPeak().x + ion_pair.second.getPeak().x -
((MHPLUS * ion_pair.first.getCharge())) /
ion_pair.first.getCharge()));
}
void
pappso::PsmFeatures::findComplementIonPairs(const pappso::PeptideSp &peptideSp)
{
std::size_t peptide_size = peptideSp.get()->size();
std::vector<PeakIonIsotopeMatch> ion_isotope_list(
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().begin(),
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().end());
for(const pappso::PeakIonIsotopeMatch &peak_ion_ext : ion_isotope_list)
{
if(peptideIonIsNter(peak_ion_ext.getPeptideIonType()))
{
auto it = findComplementIonType(ion_isotope_list.begin(),
ion_isotope_list.end(),
peak_ion_ext,
peptide_size);
if(it != ion_isotope_list.end())
{ // contains the complementary ion
m_peakIonPairs.push_back({peak_ion_ext, *it});
}
}
}
}
LinearRegression
pappso::PsmFeatures::getIonIsotopeLinearRegression() const
{
std::size_t peptide_size = msp_peptide.get()->size();
std::vector<PeakIonIsotopeMatch> ion_isotope_list(
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().begin(),
msp_peptideSpectrumMatch.get()->getPeakIonIsotopeMatchList().end());
Trace scaterplot;
for(PeptideIon ion_type : m_ionList)
{
std::vector<double> mono_th_intensities(peptide_size, 0);
std::vector<double> isotope_th_intensities(peptide_size, 0);
std::vector<double> mono_exp_intensities(peptide_size, 0);
std::vector<double> isotope_exp_intensities(peptide_size, 0);
for(const PeakIonIsotopeMatch &peak_ion_match : ion_isotope_list)
{
if(peak_ion_match.getPeptideIonType() == ion_type)
{
std::size_t vector_position =
peak_ion_match.getPeptideFragmentIonSp().get()->size() - 1;
PeptideNaturalIsotopeAverageSp iso_average_sp =
peak_ion_match.getPeptideNaturalIsotopeAverageSp();
if(iso_average_sp.get()->getIsotopeNumber() == 0)
{
mono_th_intensities[vector_position] =
iso_average_sp.get()->getIntensityRatio();
mono_exp_intensities[vector_position] =
peak_ion_match.getPeak().y;
}
else if(iso_average_sp.get()->getIsotopeNumber() == 1)
{
isotope_th_intensities[vector_position] =
iso_average_sp.get()->getIntensityRatio();
isotope_exp_intensities[vector_position] =
peak_ion_match.getPeak().y;
}
}
}
for(std::size_t i = 0; i < mono_th_intensities.size(); i++)
{
if((mono_th_intensities[i] != 0) && (isotope_th_intensities[i] != 0))
{
DataPoint xy(mono_th_intensities[i] / isotope_th_intensities[i],
mono_exp_intensities[i] /
isotope_exp_intensities[i]);
scaterplot.push_back(xy);
}
}
}
scaterplot.sortX();
LinearRegression linear_regression(scaterplot);
return linear_regression;
}
pappso::LinearRegression
Definition linearregression.h:33
pappso::MassSpectrum
Class to represent a mass spectrum.
Definition massspectrum.h:73
pappso::PeakIonIsotopeMatch
Definition peakionisotopematch.h:47
pappso::PeakIonMatch
Definition peakionmatch.h:40
pappso::PrecisionBase
Definition precision.h:44
pappso::PrecisionFactory::getPrecisionPtrFractionInstance
static PrecisionPtr getPrecisionPtrFractionInstance(PrecisionPtr origin, double fraction)
get the fraction of an existing precision pointer
Definition precision.cpp:206
pappso::PsmFeatures::getNumberOfMatchedIons
std::size_t getNumberOfMatchedIons() const
number of matched ions (peaks)
Definition psmfeatures.cpp:206
pappso::PsmFeatures::findComplementIonPairs
void findComplementIonPairs(const pappso::PeptideSp &peptideSp)
Definition psmfeatures.cpp:371
pappso::PsmFeatures::getTotalIntensity
double getTotalIntensity() const
sum of all peak intensities (matched or not)
Definition psmfeatures.cpp:161
pappso::PsmFeatures::getMatchedMzDiffMean
double getMatchedMzDiffMean() const
get mean deviation of matched peak mass delta
Definition psmfeatures.cpp:199
pappso::PsmFeatures::getTotalIntensityOfMatchedIonComplementPairs
double getTotalIntensityOfMatchedIonComplementPairs() const
intensity of matched ion complement
Definition psmfeatures.cpp:180
pappso::PsmFeatures::m_peakIonPairs
std::vector< std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > > m_peakIonPairs
Definition psmfeatures.h:161
pappso::PsmFeatures::msp_filterChargeDeconvolution
std::shared_ptr< FilterChargeDeconvolution > msp_filterChargeDeconvolution
Definition psmfeatures.h:143
pappso::PsmFeatures::getPeakIonPairs
const std::vector< std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > > & getPeakIonPairs() const
Definition psmfeatures.cpp:174
pappso::PsmFeatures::msp_filterKeepGreater
std::shared_ptr< FilterResampleKeepGreater > msp_filterKeepGreater
Definition psmfeatures.h:145
pappso::PsmFeatures::msp_filterMzExclusion
std::shared_ptr< FilterMzExclusion > msp_filterMzExclusion
Definition psmfeatures.h:144
pappso::PsmFeatures::getIonPairPrecursorMassDelta
double getIonPairPrecursorMassDelta(const std::pair< pappso::PeakIonIsotopeMatch, pappso::PeakIonIsotopeMatch > &ion_pair) const
Definition psmfeatures.cpp:351
pappso::PsmFeatures::m_precursorTheoreticalMz
double m_precursorTheoreticalMz
Definition psmfeatures.h:155
pappso::PsmFeatures::getMaxIntensityMatchedIonComplementPairPrecursorMassDelta
double getMaxIntensityMatchedIonComplementPairPrecursorMassDelta() const
get the precursor mass delta of the maximum intensity pair of complement ions
Definition psmfeatures.cpp:330
pappso::PsmFeatures::countMatchedIonComplementPairs
std::size_t countMatchedIonComplementPairs() const
count the number of matched ion complement
Definition psmfeatures.cpp:167
pappso::PsmFeatures::getComplementPairsAaSequenceCoverage
std::size_t getComplementPairsAaSequenceCoverage()
number of amino acid covered by matched complement pairs of ions
Definition psmfeatures.cpp:290
pappso::PsmFeatures::m_matchedMzDiffMean
double m_matchedMzDiffMean
Definition psmfeatures.h:163
pappso::PsmFeatures::getTotalIntensityOfMatchedIons
double getTotalIntensityOfMatchedIons() const
sum of matched peak intensities
Definition psmfeatures.cpp:150
pappso::PsmFeatures::getIonIsotopeLinearRegression
LinearRegression getIonIsotopeLinearRegression() const
Definition psmfeatures.cpp:395
pappso::PsmFeatures::m_parentCharge
unsigned int m_parentCharge
Definition psmfeatures.h:157
pappso::PsmFeatures::getMaxIntensityPeakIonMatch
double getMaxIntensityPeakIonMatch(PeptideIon ion_type) const
Definition psmfeatures.cpp:308
pappso::PsmFeatures::m_ionList
std::list< PeptideIon > m_ionList
Definition psmfeatures.h:151
pappso::PsmFeatures::m_matchedMzDiffMedian
double m_matchedMzDiffMedian
Definition psmfeatures.h:164
pappso::PsmFeatures::msp_peptideSpectrumMatch
std::shared_ptr< PeptideIsotopeSpectrumMatch > msp_peptideSpectrumMatch
Definition psmfeatures.h:147
pappso::PsmFeatures::getMaxConsecutiveIon
std::size_t getMaxConsecutiveIon(PeptideIon ion_type)
get the maximum consecutive fragments of one ion type
Definition psmfeatures.cpp:212
pappso::PsmFeatures::setPeptideSpectrumCharge
void setPeptideSpectrumCharge(const pappso::PeptideSp peptideSp, const MassSpectrum *p_spectrum, unsigned int parent_charge, unsigned int max_isotope_number)
Definition psmfeatures.cpp:58
pappso::PsmFeatures::m_ms2precision
PrecisionPtr m_ms2precision
Definition psmfeatures.h:150
pappso::PsmFeatures::getIntensityOfMatchedIon
double getIntensityOfMatchedIon(PeptideIon ion_type)
get the sum of intensity of a specific ion
Definition psmfeatures.cpp:137
pappso::PsmFeatures::msp_peptide
pappso::PeptideSp msp_peptide
Definition psmfeatures.h:148
pappso::PsmFeatures::m_precursorTheoreticalMass
double m_precursorTheoreticalMass
Definition psmfeatures.h:156
pappso::PsmFeatures::PsmFeatures
PsmFeatures(PrecisionPtr ms2precision, double minimumMz)
compute psm features
Definition psmfeatures.cpp:35
pappso::PsmFeatures::getMatchedMzDiffSd
double getMatchedMzDiffSd() const
get standard deviation of matched peak mass delta
Definition psmfeatures.cpp:193
pappso::PsmFeatures::m_spectrumSumIntensity
double m_spectrumSumIntensity
Definition psmfeatures.h:153
pappso::PsmFeatures::m_matchedMzDiffSd
double m_matchedMzDiffSd
Definition psmfeatures.h:165
pappso::PsmFeatures::getAaSequenceCoverage
std::size_t getAaSequenceCoverage(PeptideIon ion_type)
number of amino acid covered by matched ions
Definition psmfeatures.cpp:270
pappso::PsmFeatures::~PsmFeatures
~PsmFeatures()
Definition psmfeatures.cpp:53
pappso::Trace
A simple container of DataPoint instances.
Definition trace.h:148
pappso::Trace::sortX
void sortX(SortOrder sort_order=SortOrder::ascending)
Definition trace.cpp:1086
pappso::Trace::sumY
pappso_double sumY() const
Definition trace.cpp:1023
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappso::PeptideIon
PeptideIon
PeptideIon enum defines all types of ions (Nter or Cter)
Definition types.h:425
pappso::PeptideIon::y
@ y
Cter amino ions.
pappso::PeptideIon::b
@ b
Nter acylium ions.
pappso::MzFormat::xy
@ xy
(x,y) format
pappso::PeptideIonNter::a
@ a
pappso::PeptideIonNter::b
@ b
pappso::PeptideSp
std::shared_ptr< const Peptide > PeptideSp
Definition aamodification.h:48
pappso::MHPLUS
const pappso_double MHPLUS(1.007276466879)
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition types.h:50
pappso::PeptideNaturalIsotopeAverageSp
std::shared_ptr< const PeptideNaturalIsotopeAverage > PeptideNaturalIsotopeAverageSp
Definition peptidenaturalisotopeaverage.h:40
pappso::findComplementIonType
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...
Definition peakionisotopematch.cpp:30
pappso::peptideIonIsNter
bool peptideIonIsNter(PeptideIon ion_type)
tells if an ion is Nter
Definition peptide.cpp:60
pappso::XicExtractMethod::sum
@ sum
sum of intensities
pappso::XicExtractMethod::max
@ max
maximum of intensities
psmfeatures.h
comutes various PSM (Peptide Spectrum Match) features
pappso::DataPoint
Definition datapoint.h:22