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(Enums::PeptideIon::y);
  m_ionList.push_back(Enums::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(Enums::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::Enums::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::Enums::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::Enums::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(Enums::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;
}