posttreatment.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specglob/sgposttreatment.cpp
 * \date 15/11/2023
 * \author Olivier Langella
 * \brief SpecGlobTool peptide model post treatment
 *
 * 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 "posttreatment.h"
 
 
namespace pappso
{
namespace specglob
{
PostTreatment::PostTreatment(pappso::PrecisionPtr precision, const PeptideModel &peptide_model)
  : m_originalPeptideModel(peptide_model), m_betterPeptideModel(peptide_model)
{
  m_precision = precision;
  m_originalPeptideModel.matchExperimentalPeaks(m_precision);
 
  findBetterPeptideModel();
}
 
PostTreatment::~PostTreatment()
{
}
 
const PeptideModel &
PostTreatment::getBetterPeptideModel() const
{
  return m_betterPeptideModel;
}
 
void
PostTreatment::findBetterPeptideModel()
{
  if(std::abs(m_originalPeptideModel.getMassDelta()) > m_precision->getNominal())
    {
      // calculate the number of shared peaks with the deltaM on the last aa
      m_betterPeptideModel.assignResidualMass2Cter();
      m_betterPeptideModel.matchExperimentalPeaks(m_precision);
      if(m_betterPeptideModel.getCountSharedPeaks() > m_originalPeptideModel.getCountSharedPeaks())
        {
        }
      else
        {
          // not better, get back to original peptide
          m_betterPeptideModel.copyDeep(m_originalPeptideModel);
        }
    }
  else
    {
      m_betterPeptideModel.matchExperimentalPeaks(m_precision);
    }
 
 
  // Try to remove complementary mass offset
  // Not done if the number of peaks is increased by more than two peaks
  // If the number of peaks > +2, it is considered as an information that must
  // be kept in the alignment Could be modified...the information is in the
  // preAligned column
  PeptideModel test_peptide_model_remove_complementary(m_betterPeptideModel);
  if(test_peptide_model_remove_complementary.eliminateComplementaryDelta(m_precision))
    {
 
      test_peptide_model_remove_complementary.matchExperimentalPeaks(m_precision);
 
      qDebug() << "test_peptide_model_remove_complementary.getCountSharedPeaks()="
               << test_peptide_model_remove_complementary.getCountSharedPeaks()
               << " m_betterPeptideModel.getCountSharedPeaks()="
               << m_betterPeptideModel.getCountSharedPeaks();
      if(test_peptide_model_remove_complementary.getCountSharedPeaks() <
         m_betterPeptideModel.getCountSharedPeaks() + 2)
        {
          qDebug() << "replace with " << test_peptide_model_remove_complementary.toString();
          m_betterPeptideModel.copyDeep(test_peptide_model_remove_complementary);
        }
    }
 
  // if bestModified contains negative offSet, try to remove amino acids to
  // explain it
  PeptideModel test_peptide_model_remove_negative = m_betterPeptideModel;
  if(test_peptide_model_remove_negative.eliminateNegativeOffset(m_precision))
    {
      test_peptide_model_remove_negative.matchExperimentalPeaks(m_precision);
      if(test_peptide_model_remove_negative.getCountSharedPeaks() >=
         m_betterPeptideModel.getCountSharedPeaks())
        {
          qDebug() << "replace with " << test_peptide_model_remove_negative.toString();
          m_betterPeptideModel.copyDeep(test_peptide_model_remove_negative);
        }
    }
 
  // findBetterMassDifferencePosition
  tryToRemoveOffsets();
  tryBetterPositionOffsets();
 
 
  if(m_betterPeptideModel.getMassDelta() < -m_precision->getNominal())
    {
      tryToCumulateOffSets();
    }
 
  m_betterPeptideModel.removeBracketsForAlignedAA();
}
 
void
PostTreatment::tryBetterPositionOffsets()
{
 
  PeptideModel better_peptide_model_move_offset = m_betterPeptideModel;
  std::size_t best_count_shared_peaks           = m_betterPeptideModel.getCountSharedPeaks();
  qDebug() << best_count_shared_peaks;
 
  qDebug() << m_betterPeptideModel.toString();
  bool move = false;
  for(std::size_t i = 0; i < m_betterPeptideModel.size(); i++)
    {
 
      qDebug() << "i =" << i;
      if(m_betterPeptideModel.at(i).remove == true)
        continue;
      if(m_betterPeptideModel.at(i).mass_difference == 0.0)
        continue;
      if(m_betterPeptideModel.at(i).bracket == true)
        {
          // try to move if
          double mass_diff                            = m_betterPeptideModel.at(i).mass_difference;
          PeptideModel test_peptide_model_move_offset = better_peptide_model_move_offset;
          test_peptide_model_move_offset.at(i).mass_difference = 0;
          qDebug() << test_peptide_model_move_offset.toString();
 
          std::size_t j = i;
          while(j > 0)
            {
              j -= 1;
              if(test_peptide_model_move_offset.at(j).bracket == false)
                break;
              if(test_peptide_model_move_offset.at(j).mass_difference != 0.0)
                break;
 
              test_peptide_model_move_offset.at(j).mass_difference = mass_diff;
              test_peptide_model_move_offset.matchExperimentalPeaks(m_precision);
 
              qDebug() << test_peptide_model_move_offset.toString() << " "
                       << test_peptide_model_move_offset.getCountSharedPeaks() << " "
                       << best_count_shared_peaks;
              if(test_peptide_model_move_offset.getCountSharedPeaks() > best_count_shared_peaks)
                {
                  best_count_shared_peaks = test_peptide_model_move_offset.getCountSharedPeaks();
                  better_peptide_model_move_offset.copyDeep(test_peptide_model_move_offset);
                  move = true;
                }
              else
                {
                  test_peptide_model_move_offset = better_peptide_model_move_offset;
                }
            }
        }
    }
 
  if(move)
    {
      m_betterPeptideModel.copyDeep(better_peptide_model_move_offset);
    }
  qDebug() << m_betterPeptideModel.toString();
}
 
 
void
PostTreatment::tryToRemoveOffsets()
{
 
  // We'll try to evaluate for each deltaM if it can be a neutral loss or not
  // For this, we compare the number of shared peaks with and without the shift
 
  PeptideModel better_peptide_model_for_neutral_loss;
  std::size_t best_count_shared_peaks = m_betterPeptideModel.getCountSharedPeaks();
  qDebug() << best_count_shared_peaks;
  for(std::size_t i = 0; i < m_betterPeptideModel.size(); i++)
    {
      if(m_betterPeptideModel.at(i).remove == true)
        continue;
      if(m_betterPeptideModel.at(i).mass_difference == 0.0)
        continue;
 
 
      // try in the case we remove the offset
      PeptideModel test_peptide_model_for_neutral_loss          = m_betterPeptideModel;
      test_peptide_model_for_neutral_loss.at(i).mass_difference = 0;
      test_peptide_model_for_neutral_loss.matchExperimentalPeaks(m_precision);
 
      qDebug() << test_peptide_model_for_neutral_loss.toString() << " "
               << test_peptide_model_for_neutral_loss.getCountSharedPeaks() << " "
               << best_count_shared_peaks;
      if(test_peptide_model_for_neutral_loss.getCountSharedPeaks() >= best_count_shared_peaks)
        {
          best_count_shared_peaks = test_peptide_model_for_neutral_loss.getCountSharedPeaks();
          better_peptide_model_for_neutral_loss.copyDeep(test_peptide_model_for_neutral_loss);
        }
    }
 
  if(better_peptide_model_for_neutral_loss.size() > 0)
    {
      m_betterPeptideModel.copyDeep(better_peptide_model_for_neutral_loss);
    }
}
 
void
PostTreatment::tryToCumulateOffSets()
{
  double mass_delta                   = m_betterPeptideModel.getMassDelta();
  std::size_t best_count_shared_peaks = m_betterPeptideModel.getCountSharedPeaks();
  PeptideModel better_peptide_model_cumulating_residual_mass_delta;
  for(std::size_t i = 0; i < m_betterPeptideModel.size(); i++)
    {
      if(m_betterPeptideModel.at(i).remove == true)
        continue;
      if(m_betterPeptideModel.at(i).mass_difference == 0.0)
        continue;
 
      PeptideModel test_peptide_model_cumulating_residual_mass_delta = m_betterPeptideModel;
      // try to cumulate the current mass difference with unexplained mass delta
      // :
      test_peptide_model_cumulating_residual_mass_delta.at(i).mass_difference += mass_delta;
      test_peptide_model_cumulating_residual_mass_delta.matchExperimentalPeaks(m_precision);
      if(test_peptide_model_cumulating_residual_mass_delta.getCountSharedPeaks() >=
         best_count_shared_peaks)
        {
          best_count_shared_peaks =
            test_peptide_model_cumulating_residual_mass_delta.getCountSharedPeaks();
          better_peptide_model_cumulating_residual_mass_delta.copyDeep(
            test_peptide_model_cumulating_residual_mass_delta);
        }
    }
 
  if(better_peptide_model_cumulating_residual_mass_delta.size() > 0)
    {
      m_betterPeptideModel.copyDeep(better_peptide_model_cumulating_residual_mass_delta);
    }
}
 
 
const PeptideModel &
PostTreatment::getOriginalPeptideModel() const
{
  return m_originalPeptideModel;
}
 
bool
PostTreatment::findReplaceMutations()
{
  auto peptide_model = m_betterPeptideModel;
 
  bool modif                                = false;
  std::vector<Enums::AminoAcidChar> aa_list = {Enums::AminoAcidChar::alanine,
                                               Enums::AminoAcidChar::aspartic_acid,
                                               Enums::AminoAcidChar::glutamic_acid,
                                               Enums::AminoAcidChar::phenylalanine,
                                               Enums::AminoAcidChar::glycine,
                                               Enums::AminoAcidChar::histidine,
                                               Enums::AminoAcidChar::isoleucine,
                                               Enums::AminoAcidChar::lysine,
                                               Enums::AminoAcidChar::leucine,
                                               Enums::AminoAcidChar::methionine,
                                               Enums::AminoAcidChar::asparagine,
                                               Enums::AminoAcidChar::proline,
                                               Enums::AminoAcidChar::glutamine,
                                               Enums::AminoAcidChar::arginine,
                                               Enums::AminoAcidChar::serine,
                                               Enums::AminoAcidChar::threonine,
                                               Enums::AminoAcidChar::valine,
                                               Enums::AminoAcidChar::tryptophan,
                                               Enums::AminoAcidChar::tyrosine};
  if(peptide_model.checkForMutations(aa_list, m_precision))
    {
      modif = true;
    }
  pappso::Aa cysteine('C');
  cysteine.addAaModification(pappso::AaModification::getInstance("MOD:00397")); // carbamido
  if(peptide_model.checkForMutation(cysteine, m_precision))
    {
      modif = true;
    }
  pappso::Aa methionine('M');
  methionine.addAaModification(pappso::AaModification::getInstance("MOD:00719")); // oxydation
  if(peptide_model.checkForMutation(methionine, m_precision))
    {
      modif = true;
    }
  if(modif)
    {
      m_betterPeptideModel = peptide_model;
    }
  return modif;
}
} // namespace specglob
} // namespace pappso