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);
                }
              if(m_reindexed_peaks.at(peak1) >= 0)
                {
                  addAaPosition(aa, peak2, peak1, true);
                }
              else
                {
                  addAaPosition(aa, peak2, next_l_peak, 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;
        }
    }
 
  // 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,
                                                    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.
  if(l_support)
    {
      m_aapositions.at(aa - 1)->push_back(
        {r_peak, l_peak, computeCondition(l_peak, l_support), l_support});
    }
  else
    {
      m_aapositions.at(aa - 1)->push_back(
        {r_peak, l_peak, computeCondition(l_peak, l_support), l_support});
    }
}
 
uint32_t
pappso::specpeptidoms::SpOMSSpectrum::computeCondition(const std::size_t l_peak,
                                                       bool l_support) const
{
  uint32_t condition;
  if(l_peak == 0)
    {
      condition = 2;
    }
  else if(!l_support)
    {
      condition = 1;
    }
  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