peptidemodel.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specglob/peptidemodel.h
 * \date 14/11/2023
 * \author Olivier Langella
 * \brief SpecGlobTool peptide model
 *
 * 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 ipetide to spectrum
 * alignmentt 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 "peptidemodel.h"
#include "../../peptide/peptiderawfragmentmasses.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/utils.h"
 
 
namespace pappso
{
namespace specglob
{
PeptideModel::PeptideModel()
{
  // not valid
}
 
PeptideModel::PeptideModel(const pappso::QualifiedMassSpectrum &qmass_spectrum,
                           const pappso::Peptide &peptide)
  : std::vector<AminoAcidModel>()
{
  m_experimentalPrecursorMass = qmass_spectrum.getPrecursorMass();
  mcsp_peakList               = qmass_spectrum.getMassSpectrumCstSPtr();
  m_nterModification          = peptide.getCleavageNterModification();
  m_cterModification          = peptide.getCleavageCterModification();
  for(auto aa_pep : peptide)
    {
 
      push_back({aa_pep, SpectralAlignmentType::nonAlign, false, 0.0});
    }
 
  // internal:Nter_hydrolytic_cleavage_H
 
  // internal:Cter_hydrolytic_cleavage_HO
  //
}
 
PeptideModel::PeptideModel(const PeptideModel &other) : std::vector<AminoAcidModel>(other)
{
  m_experimentalPrecursorMass = other.m_experimentalPrecursorMass;
  mcsp_peakList               = other.mcsp_peakList;
 
  m_experimentalPrecursorMass = other.m_experimentalPrecursorMass;
  m_beginMassDelta            = other.m_beginMassDelta;
 
  m_nterModification           = other.m_nterModification;
  m_cterModification           = other.m_cterModification;
  m_countSharedPeaks           = other.m_countSharedPeaks;
  m_intensitySharedPeaks       = other.m_intensitySharedPeaks;
  m_intensityExperimentalPeaks = other.m_intensityExperimentalPeaks;
}
 
PeptideModel &
PeptideModel::operator=(const PeptideModel &other)
{
 
  assign(other.begin(), other.end());
  m_experimentalPrecursorMass = other.m_experimentalPrecursorMass;
  mcsp_peakList               = other.mcsp_peakList;
 
  m_experimentalPrecursorMass = other.m_experimentalPrecursorMass;
  m_beginMassDelta            = other.m_beginMassDelta;
 
  m_nterModification           = other.m_nterModification;
  m_cterModification           = other.m_cterModification;
  m_countSharedPeaks           = other.m_countSharedPeaks;
  m_intensitySharedPeaks       = other.m_intensitySharedPeaks;
  m_intensityExperimentalPeaks = other.m_intensityExperimentalPeaks;
  m_theoreticalPeakList.clear();
  return *this;
}
 
PeptideModel &
PeptideModel::copyDeep(const PeptideModel &other)
{
  operator=(other);
 
  m_theoreticalPeakList = other.m_theoreticalPeakList;
  return *this;
}
 
PeptideModel::~PeptideModel()
{
}
 
double
PeptideModel::getMz(unsigned int charge) const
{
  return ((getMass() + (pappso::MHPLUS * charge)) / (double)charge);
}
 
QString
PeptideModel::toString() const
{
  QString peptide_str;
 
  if(m_beginMassDelta != 0.0)
    peptide_str.append(QString("[%1]").arg(QString::number(m_beginMassDelta, 'f', 2)));
 
  for(auto aa_model : *this)
    {
      qDebug() << " aa_model.amino_acid.toString()=" << aa_model.amino_acid.toAbsoluteString()
               << " aa_model.alignment_type=" << Utils::toString(aa_model.alignment_type)
               << " aa_model.mass_difference=" << aa_model.mass_difference
               << " aa_model.bracket=" << aa_model.bracket;
 
      if(aa_model.remove)
        peptide_str.append("-(");
      if(aa_model.bracket)
        {
          peptide_str.append(QString("[%1]").arg(aa_model.amino_acid.toString()));
        }
      else
        {
          peptide_str.append(aa_model.amino_acid.toString());
        }
 
 
      if(aa_model.mass_difference != 0.0)
        {
          peptide_str.append(
            QString("[%1]").arg(QString::number(aa_model.mass_difference, 'f', 2)));
        }
 
      if(aa_model.remove)
        peptide_str.append(")");
    }
 
  peptide_str.replace(")-(", "");
  peptide_str.append(QString("_[%1]").arg(QString::number(getMassDelta(), 'f', 2)));
  peptide_str.replace("_[0.00]", "");
  peptide_str.replace("_[-0.00]", "");
  return peptide_str;
}
 
 
bool
PeptideModel::hasRemoval() const
{
  for(auto aa_model : *this)
    {
      if(std::abs(aa_model.amino_acid.getMass() + aa_model.mass_difference) < 0.001)
        return true;
    }
  return false;
}
 
bool
PeptideModel::ltrimOnRemoval()
{
  AminoAcidModel *amino_acid_to_remove = nullptr;
  int pos                              = -1;
  for(auto aa_model : *this)
    {
      pos++;
      if(std::abs(aa_model.amino_acid.getMass() + aa_model.mass_difference) < 0.001)
        {
          amino_acid_to_remove = &aa_model;
          break;
        }
    }
  if(amino_acid_to_remove != nullptr)
    {
      if(std::abs(front().amino_acid.getMass() - amino_acid_to_remove->amino_acid.getMass()) <
         0.001)
        {
          erase(begin());
          if(pos > 0)
            at(pos - 1).mass_difference = 0;
          return true;
        }
    }
  return false;
}
 
 
QString
PeptideModel::toProForma() const
{
  QString peptide_str("[");
 
  // if(m_beginMassDelta != 0.0)
  //   peptide_str.append(
  //    QString("[%1]").arg(QString::number(m_beginMassDelta, 'f', 2)));
  // ProForma version 2 unexplained mass delta :
  //  peptide_str.append(
  //  QString("_[%1]").arg(QString::number(getMassDelta(), 'f', 2)));
 
  if(getMassDelta() > 0)
    {
      peptide_str.append("+");
    }
  peptide_str.append(QString("%1]?").arg(QString::number(getMassDelta(), 'f', 4)));
  peptide_str.replace("[+0.0000]?", "");
  peptide_str.replace("[-0.0000]?", "");
 
 
  double mass_diff = m_beginMassDelta;
  for(auto aa_model : *this)
    {
      qDebug() << " aa_model.amino_acid.toString()=" << aa_model.amino_acid.toAbsoluteString()
               << " aa_model.alignment_type=" << Utils::toString(aa_model.alignment_type)
               << " aa_model.mass_difference=" << aa_model.mass_difference
               << " aa_model.bracket=" << aa_model.bracket;
 
      if(aa_model.bracket)
        {
        }
 
      if(aa_model.remove)
        {
          peptide_str.append(aa_model.amino_acid.getLetter());
 
          peptide_str.append(
            "[" +
            AaModification::getInstanceRemovalAccessionByAaLetter(aa_model.amino_acid.getLetter())
              ->getAccession() +
            "]");
        }
      else
        {
          peptide_str.append(aa_model.amino_acid.toProForma());
 
          mass_diff += aa_model.mass_difference;
          if(mass_diff != 0.0)
            {
              if(mass_diff > 0)
                {
                  peptide_str.append(QString("[+%1]").arg(QString::number(mass_diff, 'f', 4)));
                }
              else
                {
                  peptide_str.append(QString("[%1]").arg(QString::number(mass_diff, 'f', 4)));
                }
            }
        }
      mass_diff = 0.0;
    }
  return peptide_str;
}
double
PeptideModel::getMass() const
{
  double mass = m_beginMassDelta;
  if(m_cterModification != nullptr)
    mass += m_cterModification->getMass();
  if(m_nterModification != nullptr)
    mass += m_nterModification->getMass();
 
  for(auto aa_model : *this)
    {
      mass += aa_model.amino_acid.getMass() + aa_model.mass_difference;
    }
  return mass;
}
 
double
PeptideModel::getMassDelta() const
{
  return m_experimentalPrecursorMass - getMass();
}
 
std::size_t
PeptideModel::modifCount() const
{
  std::size_t count = 0;
  if(m_beginMassDelta != 0.0)
    count++;
  for(auto aa_model : *this)
    {
      if(aa_model.mass_difference != 0.0)
        count++;
    }
  return count;
}
 
void
PeptideModel::setBeginMassDelta(double delta)
{
  m_beginMassDelta = delta;
}
 
void
PeptideModel::matchExperimentalPeaks(pappso::PrecisionPtr precision)
{
  if(m_theoreticalPeakList.size() == 0)
    {
      m_intensityExperimentalPeaks = mcsp_peakList.get()->totalIonCurrent();
      generateTheoreticalPeaks(precision);
    }
  auto it_expe_peaks_end = mcsp_peakList.get()->end();
  auto it_expe_peaks     = mcsp_peakList.get()->begin();
  auto it_theo_peaks     = m_theoreticalPeakList.begin();
  while((it_expe_peaks != mcsp_peakList.get()->end()) &&
        (it_theo_peaks != m_theoreticalPeakList.end()))
    {
      while((it_expe_peaks != mcsp_peakList.get()->end()) &&
            (it_expe_peaks->x < it_theo_peaks->mz_range.lower()))
        { // find the first possible ion hit
          it_expe_peaks++;
        }
      if(it_expe_peaks == mcsp_peakList.get()->end())
        continue;
 
      if(it_expe_peaks->x < it_theo_peaks->mz_range.upper())
        {
          // ok, we've got a hit
          auto it_theo_peaks_loop = it_theo_peaks;
          while((it_theo_peaks_loop != m_theoreticalPeakList.end()) &&
                (it_theo_peaks_loop->mz_range.contains(it_expe_peaks->x)))
            {
              if((it_theo_peaks_loop->it_experimental_peak_match == it_expe_peaks_end) ||
                 (it_theo_peaks_loop->it_experimental_peak_match->y < it_expe_peaks->y))
                {
                  it_theo_peaks_loop->it_experimental_peak_match = it_expe_peaks;
                  qDebug() << "match: mz=" << it_expe_peaks->x << " intensity=" << it_expe_peaks->y;
                  break;
                }
              it_theo_peaks_loop++;
            }
          it_expe_peaks++;
        }
      else
        {
          it_theo_peaks++;
        }
      // it_expe_peaks++;
    }
 
  // get number of shared peaks
  m_countSharedPeaks     = 0;
  m_intensitySharedPeaks = 0;
  std::vector<std::vector<pappso::DataPoint>::const_iterator> matched_experimental_peaks;
  for(auto &theo_peak : m_theoreticalPeakList)
    {
      if(theo_peak.it_experimental_peak_match != it_expe_peaks_end)
        {
          matched_experimental_peaks.push_back(theo_peak.it_experimental_peak_match);
          // m_countSharedPeaks++;
        }
    }
  auto it_end = std::unique(matched_experimental_peaks.begin(), matched_experimental_peaks.end());
  m_countSharedPeaks     = std::distance(matched_experimental_peaks.begin(), it_end);
  m_intensitySharedPeaks = std::accumulate(
    matched_experimental_peaks.begin(),
    it_end,
    0.0,
    [](double k, const std::vector<pappso::DataPoint>::const_iterator &l) { return (k + l->y); });
}
 
void
PeptideModel::generateTheoreticalPeaks(pappso::PrecisionPtr precision)
{
  m_theoreticalPeakList.clear();
  double cumul_mass =
    m_beginMassDelta + m_nterModification->getMass() - pappso::MPROTIUM + pappso::MHPLUS;
 
 
  auto it_expe_peaks_end = mcsp_peakList.get()->end();
  std::size_t i          = 0;
  for(auto aa_model : *this)
    {
      // compute B ions
      cumul_mass += aa_model.amino_acid.getMass() + aa_model.mass_difference;
 
      qDebug() << cumul_mass;
      m_theoreticalPeakList.push_back(
        {pappso::MzRange(cumul_mass, precision), it_expe_peaks_end, i});
      i++;
    }
  // remove the complete ion B (entire peptide, no fragment):
  m_theoreticalPeakList.pop_back();
 
  std::size_t bmaxindice = m_theoreticalPeakList.size();
  double yion_delta      = pappso::PeptideRawFragmentMasses::getDeltaMass(Enums::PeptideIon::y);
  cumul_mass += yion_delta + pappso::MHPLUS;
  for(std::size_t i = 0; i < bmaxindice; i++)
    {
      m_theoreticalPeakList.push_back(
        {pappso::MzRange(cumul_mass - m_theoreticalPeakList[i].mz_range.getMz(), precision),
         it_expe_peaks_end,
         i + 1});
      qDebug() << m_theoreticalPeakList.back().mz_range.getMz();
    }
 
  std::sort(m_theoreticalPeakList.begin(),
            m_theoreticalPeakList.end(),
            [](const TheoreticalPeakDataPoint &a, const TheoreticalPeakDataPoint &b) {
              return a.mz_range.getMz() < b.mz_range.getMz();
            });
}
 
std::vector<double>
PeptideModel::getTheoreticalIonMassList() const
{
  std::vector<double> mass_list;
  for(auto &datapoint : m_theoreticalPeakList)
    {
      mass_list.push_back(datapoint.mz_range.getMz());
    }
  return mass_list;
}
 
std::size_t
PeptideModel::getCountSharedPeaks() const
{
  return m_countSharedPeaks;
}
 
double
PeptideModel::getIntensitySharedPeaks() const
{
  return m_intensitySharedPeaks;
}
 
double
PeptideModel::getIntensityExperimentalPeaks() const
{
  return m_intensityExperimentalPeaks;
}
 
void
PeptideModel::assignResidualMass2Cter()
{
  back().mass_difference = getMassDelta();
}
 
bool
PeptideModel::eliminateComplementaryDelta(pappso::PrecisionPtr precision)
{
  double precision_ref                                 = precision->getNominal();
  double offset                                        = 0;
  bool modified                                        = false;
  std::vector<AminoAcidModel>::iterator it_begin_block = end();
  if(std::abs(m_beginMassDelta) > precision_ref)
    {
      it_begin_block = begin();
      offset         = m_beginMassDelta;
 
      qDebug() << " new block offset=" << offset;
    }
  qDebug() << offset;
  for(auto it = begin(); it != end(); it++)
    {
      qDebug() << " it->mass_difference=" << it->mass_difference;
      if(it_begin_block == end())
        {
 
          if(std::abs(it->mass_difference) > precision_ref)
            {
              // new block
              it_begin_block = it;
              offset         = it->mass_difference;
            }
          qDebug() << " new block offset=" << offset;
        }
      if(it_begin_block != end())
        {
          // existing block
          qDebug() << " existing block offset=" << offset;
          if(it == begin())
            {
              offset += it->mass_difference;
            }
          else
            {
              // where is the end of the block ?
              qDebug() << " where is the end of the block ? offset=" << offset;
              auto it_end_block =
                std::find_if(it, end(), [offset, precision_ref](const AminoAcidModel &point) {
                  qDebug() << point.mass_difference << " " << offset + point.mass_difference;
                  if(std::abs(offset + point.mass_difference) <= precision_ref)
                    return true;
                  return false;
                });
 
              if(it_end_block != end())
                {
                  // eliminate mass mass_difference
                  modified = true;
                  if(it_begin_block == begin())
                    m_beginMassDelta = 0;
                  it_begin_block->mass_difference = 0.0;
                  it_end_block->mass_difference   = 0.0;
                  // start new block search
                  it_begin_block = end();
                  offset         = 0;
                  qDebug() << " existing block end block found";
                }
              else
                {
                  qDebug() << " existing block end block not found offset=" << offset;
                  // close block
                  it_begin_block = end();
                  offset         = 0;
                }
            }
        }
    }
  qDebug() << toString();
  return modified;
}
 
 
bool
PeptideModel::eliminateNegativeOffset(pappso::PrecisionPtr precision)
{ // System.out.println("Elimination des neg offsets");
  double precision_ref = precision->getNominal();
  auto it_reverse      = rbegin();
  bool modified        = false;
 
  while(it_reverse != rend())
    {
      if(it_reverse->remove)
        {
          // this is a part of already removed block
        }
      else
        {
          double mass_difference = it_reverse->mass_difference;
          if(mass_difference < 0)
            {
              // try to find an amino acid blocl of that mass
              double cumul_mass_aa = 0;
              auto it_block_end    = it_reverse;
              while(it_block_end != rend())
                {
                  cumul_mass_aa += it_block_end->amino_acid.getMass();
                  it_block_end++;
                  if(std::abs(mass_difference + cumul_mass_aa) < precision_ref)
                    {
                      modified = true;
                      // block found -> remove this range
                      for(auto it_block = it_reverse; it_block != it_block_end; it_block++)
                        {
                          it_block->remove = true;
                        }
                    }
                  else
                    {
                      if((it_block_end != rend()) && (it_block_end->mass_difference != 0.0))
                        {
                          // interrupted block
                          break;
                        }
                    }
                }
            }
        }
      it_reverse++;
    }
 
  return modified;
}
 
bool
PeptideModel::removeBracketsForAlignedAA()
{
  if(m_theoreticalPeakList.size() == 0)
    {
      throw pappso::PappsoException("missing peak matching");
    }
  bool modified          = false;
  auto it_expe_peaks_end = mcsp_peakList.get()->end();
  for(auto &theo_peak : m_theoreticalPeakList)
    {
      if(theo_peak.it_experimental_peak_match != it_expe_peaks_end)
        {
          if(at(theo_peak.aa_indice).bracket == true)
            {
              modified                        = true;
              at(theo_peak.aa_indice).bracket = false;
            }
        }
    }
  return modified;
}
 
const std::vector<TheoreticalPeakDataPoint> &
PeptideModel::getTheoreticalPeakDataPointList() const
{
  return m_theoreticalPeakList;
}
 
QString
PeptideModel::getTheoreticalPeakDataPointListToString() const
{
 
  auto it_expe_peaks_end = mcsp_peakList.get()->end();
  QString str;
  for(auto &theo_peak : m_theoreticalPeakList)
    {
      double intensity = 0;
      if(theo_peak.it_experimental_peak_match != it_expe_peaks_end)
        {
          intensity = theo_peak.it_experimental_peak_match->y;
        }
      str.append(QString("{%1,%2,%3,%4}")
                   .arg(theo_peak.mz_range.getMz())
                   .arg(intensity)
                   .arg(theo_peak.aa_indice)
                   .arg(at(theo_peak.aa_indice).amino_acid.toString()));
    }
  return str;
}
 
bool
PeptideModel::checkForMutation(const pappso::Aa &amino_acid_candidate,
                               pappso::PrecisionPtr precision)
{
  bool modif                  = false;
  double total_mass_to_insert = amino_acid_candidate.getMass();
  auto insert_aa_modification =
    AaModification::getInstanceInsertionAccessionByAaLetter(amino_acid_candidate.getLetter());
 
  double mass_difference = m_beginMassDelta;
  for(auto &aa_model : *this)
    {
      mass_difference += aa_model.mass_difference;
      auto aa_ref = aa_model.amino_acid;
      if(mass_difference != 0.0)
        {
          mass_difference += aa_ref.getTotalModificationMass();
          AaModificationP aa_ref_remove =
            AaModification::getInstanceRemovalAccessionByAaLetter(aa_model.amino_acid.getLetter());
          // check for possible mutation
          double current_diff = (aa_ref_remove->getMass() - mass_difference) + total_mass_to_insert;
 
 
          qDebug() << "current_diff=" << current_diff;
          if(std::abs(current_diff) < precision->getNominal())
            {
              qDebug() << "modif = true=";
              aa_model.mass_difference = 0;
              aa_model.amino_acid.removeAllButInternalModification();
              aa_model.amino_acid.addAaModification(aa_ref_remove);
              aa_model.amino_acid.addAaModification(insert_aa_modification);
              for(auto &mod : amino_acid_candidate.getModificationList())
                {
                  aa_model.amino_acid.addAaModification(mod);
                }
              modif = true;
            }
        }
      mass_difference = 0;
    }
  return modif;
}
 
bool
PeptideModel::checkForMutations(const std::vector<Enums::AminoAcidChar> &aa_list,
                                pappso::PrecisionPtr precision)
{
  bool modif = false;
  // AaModificationP carbamido = AaModification::getInstance("MOD:00397");
  // AaModificationP metox = AaModification::getInstance("MOD:00719");
 
  std::vector<AaModificationP> modificationInsertList;
  for(const auto aa : aa_list)
    {
      modificationInsertList.push_back(
        AaModification::getInstanceInsertionAccessionByAaLetter((char)aa));
    }
  double mass_difference = m_beginMassDelta;
  for(auto &aa_model : *this)
    {
      mass_difference += aa_model.mass_difference;
      auto aa_ref = aa_model.amino_acid;
      if(mass_difference != 0.0)
        {
          mass_difference += aa_ref.getTotalModificationMass();
          AaModificationP aa_ref_remove =
            AaModification::getInstanceRemovalAccessionByAaLetter(aa_model.amino_acid.getLetter());
          AaModificationP aa_ref_insert = nullptr;
          double better_diff            = 10;
          // check for possible mutation
          for(const auto aaInsertion : modificationInsertList)
            {
              double current_diff =
                (aa_ref_remove->getMass() - mass_difference) + aaInsertion->getMass();
              qDebug() << aa_model.amino_acid.getLetter() << " " << aaInsertion->getName() << " "
                       << aaInsertion->getMass() << " " << current_diff;
              if(std::abs(current_diff) < std::abs(better_diff))
                {
                  aa_ref_insert = aaInsertion;
                  better_diff   = current_diff;
                }
            }
          qDebug() << "better_diff=" << better_diff;
          if((aa_ref_insert != nullptr) && (std::abs(better_diff) < precision->getNominal()))
            {
 
              qDebug() << "modif = true=";
              aa_model.mass_difference = 0;
              aa_model.amino_acid.removeAllButInternalModification();
              aa_model.amino_acid.addAaModification(aa_ref_remove);
              aa_model.amino_acid.addAaModification(aa_ref_insert);
              modif = true;
            }
        }
      mass_difference = 0;
    }
  return modif;
}
 
bool
PeptideModel::checkForAaModificationP(pappso::AaModificationP modification,
                                      pappso::PrecisionPtr precision)
{
 
  bool modif = false;
 
  double mass_difference = m_beginMassDelta;
  for(auto &aa_model : *this)
    {
      mass_difference += aa_model.mass_difference;
      auto aa_ref = aa_model.amino_acid;
      if(mass_difference != 0.0)
        {
          if(std::abs(modification->getMass() - mass_difference) < precision->getNominal())
            {
              aa_model.mass_difference = 0;
              aa_model.amino_acid.addAaModification(modification);
              modif = true;
            }
        }
      mass_difference = 0;
    }
  return modif;
}
 
 
bool
PeptideModel::checkForAaModification(const pappso::Aa &aa_modified, pappso::PrecisionPtr precision)
{
 
  bool modif = false;
 
  double mass_difference = m_beginMassDelta;
  for(auto &aa_model : *this)
    {
      mass_difference += aa_model.mass_difference;
      auto aa_ref = aa_model.amino_acid;
      if(aa_ref.getLetter() == aa_modified.getLetter())
        {
          if(mass_difference != 0.0)
            {
              double mass_modification = aa_modified.getTotalModificationMass();
              if(std::abs(mass_modification - mass_difference) < precision->getNominal())
                {
                  aa_model.mass_difference = 0;
                  for(auto &modification : aa_modified.getModificationList())
                    {
                      aa_model.amino_acid.addAaModification(modification);
                    }
                  modif = true;
                }
            }
        }
      mass_difference = 0;
    }
  return modif;
}
 
} // namespace specglob
} // namespace pappso