experimentalspectrum.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specglob/experimentalspectrum.cpp
 * \date 08/11/2023
 * \author Olivier Langella
 * \brief transform a spectrum to SpecGlob spectra
 *
 * C++ implementation of the SpecGlob algorithm described in :
 * 1. Prunier, G. et al. Fast alignment of mass spectra in large proteomics
 * datasets, capturing dissimilarities arising from multiple complex
 * modifications of peptides. BMC Bioinformatics 24, 421 (2023).
 *
 *     HAL Id : hal-04296170 , version 1
 *     Mot de passe : hxo20cl
 *     DOI : 10.1186/s12859-023-05555-y
 */
 
/*
 * SpecGlobTool, Spectra to peptide alignment tool
 * Copyright (C) 2023  Olivier Langella
 * <olivier.langella@universite-paris-saclay.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 "experimentalspectrum.h"
#include <QDebug>
#include "pappsomspp/core/pappsoexception.h"
#include <QObject>
 
namespace pappso
{
namespace specglob
{
ExperimentalSpectrum::ExperimentalSpectrum(const pappso::QualifiedMassSpectrum &qmass_spectrum,
                                           pappso::PrecisionPtr precision_ptr)
  : std::vector<ExperimentalSpectrumDataPoint>()
{
  m_qualifiedMassSpectrum = qmass_spectrum;
  m_precisionPtr          = precision_ptr;
  createSymetricPeakList();
}
 
ExperimentalSpectrum::ExperimentalSpectrum(const pappso::QualifiedMassSpectrum &qmass_spectrum,
                                           pappso::PrecisionPtr precision_ptr,
                                           double precursor_mass_error)
  : std::vector<ExperimentalSpectrumDataPoint>()
{
  m_qualifiedMassSpectrum = qmass_spectrum;
  m_precisionPtr          = precision_ptr;
  createSymetricPeakList(precursor_mass_error);
}
 
ExperimentalSpectrum::ExperimentalSpectrum(const ExperimentalSpectrum &other)
  : std::vector<ExperimentalSpectrumDataPoint>(other)
{
  m_targetMzSum           = other.m_targetMzSum;
  m_qualifiedMassSpectrum = other.m_qualifiedMassSpectrum;
  m_precisionPtr          = other.m_precisionPtr;
}
 
ExperimentalSpectrum::~ExperimentalSpectrum()
{
}
 
void
ExperimentalSpectrum::createSymetricPeakList()
{
  createSymetricPeakList(0);
}
 
void
ExperimentalSpectrum::createSymetricPeakList(double precursor_mass_error)
{
 
  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()));
    }
 
  double mz_prec = m_qualifiedMassSpectrum.getPrecursorMz(&ok);
  if(!ok)
    {
      throw pappso::PappsoException(
        QObject::tr("precursor m/z is not defined in spectrum %1")
          .arg(m_qualifiedMassSpectrum.getMassSpectrumId().getNativeId()));
    }
 
  // compute precursor mass given the charge state
  m_precursorMass = mz_prec * (double)charge - precursor_mass_error;
  m_precursorMass -= pappso::MHPLUS * (double)charge;
  m_targetMzSum = m_precursorMass + pappso::MHPLUS + pappso::MHPLUS;
 
  std::vector<double> mz_list = m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->xValues();
 
  std::size_t mz_current_indice = 0;
  auto itend                    = m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->end();
  for(std::vector<pappso::DataPoint>::const_iterator it =
        m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->begin();
      it != itend;
      it++)
    {
      double current_mz = it->x;
 
      double symmetric_current_mz = getSymetricMz(current_mz);
 
      auto itpair_symmetric = findMz(symmetric_current_mz);
      if(itpair_symmetric != itend)
        {
          // there is a counterpart
          push_back({ExperimentalSpectrumDataPointType::both, it->x, 0});
          qDebug() << "current_mz=" << current_mz << " both";
        }
      else
        {
          // there is no counterpart
          push_back({ExperimentalSpectrumDataPointType::native, it->x, 0});
          push_back({ExperimentalSpectrumDataPointType::symmetric, symmetric_current_mz, 0});
          qDebug() << "current_mz=" << current_mz << " symmetrics";
        }
 
      mz_current_indice++;
    }
 
  // we add a peak with NT mass (1.0078) if it is not detected, to give better
  // chance to align the first amino acid in b if it is present
  if(findMz(pappso::MPROTIUM) == itend)
    {
      push_back({ExperimentalSpectrumDataPointType::synthetic, pappso::MPROTIUM, 0});
    }
 
  // We add the B peak corresponding to the precursor (complete peptide)
  double precusorBion = m_targetMzSum - pappso::MASSH2O - pappso::MPROTIUM;
 
  if(findMz(precusorBion) == itend)
    {
      push_back({ExperimentalSpectrumDataPointType::synthetic, precusorBion, 0});
    }
 
  std::sort(begin(),
            end(),
            [](const ExperimentalSpectrumDataPoint &a, const ExperimentalSpectrumDataPoint &b) {
              return (a.peak_mz < b.peak_mz);
            });
 
  std::size_t i = 0;
  for(auto &data_point : *this)
    {
      data_point.indice = i;
      i++;
    }
}
 
std::vector<pappso::DataPoint>::const_iterator
ExperimentalSpectrum::findMz(double mz)
{
  pappso::MzRange mz_range(mz, m_precisionPtr);
  auto itend = m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->end();
  auto it =
    findFirstEqualOrGreaterX(m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->begin(),
                             m_qualifiedMassSpectrum.getMassSpectrumCstSPtr().get()->end(),
                             mz_range.lower());
 
  if(it != itend)
    {
      if(it->x <= mz_range.upper())
        {
          return it;
        }
    }
  return itend;
}
 
double
ExperimentalSpectrum::getTargetMzSum() const
{
  return m_targetMzSum;
}
 
std::vector<double>
ExperimentalSpectrum::getMassList() const
{
  std::vector<double> mass_list;
  for(const ExperimentalSpectrumDataPoint &n : *this)
    {
      mass_list.push_back(n.peak_mz);
    }
 
  return mass_list;
}
 
double
ExperimentalSpectrum::getSymetricMz(double mz) const
{
  return m_targetMzSum - mz;
}
 
std::vector<double>
ExperimentalSpectrum::getMassList(ExperimentalSpectrumDataPointType type) const
{
  std::vector<double> mass_list;
  for(const ExperimentalSpectrumDataPoint &n : *this)
    {
      if(n.type == type)
        mass_list.push_back(n.peak_mz);
    }
 
  return mass_list;
}
 
QString
ExperimentalSpectrum::toString() const
{
  QStringList all_element;
  for(const ExperimentalSpectrumDataPoint &n : *this)
    {
      all_element << QString("%1 %2").arg(n.peak_mz).arg((std::uint8_t)n.type);
    }
  return QString("[%1]").arg(all_element.join("] ["));
}
 
double
ExperimentalSpectrum::getPrecursorMass() const
{
  return m_precursorMass;
}
 
std::vector<ExperimentalSpectrumDataPoint>::const_reverse_iterator
ExperimentalSpectrum::reverseFindDiffMz(std::size_t start_position,
                                        const pappso::MzRange &aaTheoMzRange) const
{
 
  qDebug() << "start_position" << start_position << " lookfor=" << aaTheoMzRange.getMz();
  std::vector<ExperimentalSpectrumDataPoint>::const_reverse_iterator itrbegin =
    rbegin() + (size() - 1 - start_position);
 
  qDebug() << itrbegin->indice << "mz=" << itrbegin->peak_mz;
  double eperimentalMzReference = itrbegin->peak_mz;
  auto itrend                   = this->rend();
  // We check all row from j to 0
  qDebug();
  for(auto itr = itrbegin + 1; itr != itrend; ++itr)
    {
      qDebug() << itr->indice;
      double experimentalMzDifference = eperimentalMzReference - itr->peak_mz;
 
      if(experimentalMzDifference > aaTheoMzRange.upper())
        {
          // if we pass the mass of the theoretical amino acid, we stop to not
          // over
          // calculate
          qDebug() << experimentalMzDifference << ">" << aaTheoMzRange.upper();
          return itrend;
        }
      else if(experimentalMzDifference < aaTheoMzRange.lower())
        {
          continue;
        }
 
      qDebug() << itr->indice << " diff=" << experimentalMzDifference;
      // if we found that j-k give an amino acid, we keep k value
      return itr;
    }
  qDebug() << "rend";
  return itrend; // a value of -1 to show that
                 // there is no k value
}
 
const pappso::QualifiedMassSpectrum &
ExperimentalSpectrum::getQualifiedMassSpectrum() const
{
  return m_qualifiedMassSpectrum;
}
} // namespace specglob
} // namespace pappso