spectralalignment.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specglob/spectraalignment.cpp
 * \date 08/11/2023
 * \author Olivier Langella
 * \brief petide to spectrum alignment
 *
 * 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 "spectralalignment.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/utils.h"
 
 
namespace pappso
{
namespace specglob
{
SpectralAlignment::SpectralAlignment(ScoreValues score_values, pappso::PrecisionPtr precision_ptr)
{
  m_precisionPtr = precision_ptr;
  m_scoreValues  = score_values;
}
 
SpectralAlignment::~SpectralAlignment()
{
}
 
 
ExperimentalSpectrumCsp
SpectralAlignment::getExperimentalSpectrumCsp() const
{
  return mcsp_experimentalSpectrum;
}
 
PeptideSpectraCsp
SpectralAlignment::getPeptideSpectraCsp() const
{
  return mcsp_peptideSpectrum;
}
void
SpectralAlignment::align(PeptideSpectraCsp peptide_spectrum,
                         ExperimentalSpectrumCsp experimental_spectrum)
{
  mcsp_experimentalSpectrum = experimental_spectrum;
  mcsp_peptideSpectrum      = peptide_spectrum;
  m_maxScore                = 0;
  m_matrix.resize(mcsp_peptideSpectrum.get()->size(), mcsp_experimentalSpectrum.get()->size());
 
  // set the value of difference between Precursor mass and theoretical mass of
  // the proposed peptide
  m_precursorMassDelta = mcsp_experimentalSpectrum.get()->getPrecursorMass() -
                         mcsp_peptideSpectrum.get()->getPeptideSp().get()->getMass();
 
  fillMassDelta(*mcsp_peptideSpectrum.get(), *mcsp_experimentalSpectrum.get());
 
  m_itPosMax = m_matrix.end2();
  // Make the alignment in filling the matrices
 
  auto last_theoretical_peptide = std::prev(m_matrix.end1());
  for(auto itmi = ++m_matrix.begin1(); itmi != m_matrix.end1(); ++itmi)
    {
      for(auto itmj = ++itmi.begin(); itmj != itmi.end(); itmj++)
        {
          // qDebug() << "i=" << itmj.index1() << " j=" << itmj.index2();
          fillMatricesWithScores(
            itmj, *mcsp_peptideSpectrum.get(), *mcsp_experimentalSpectrum.get());
          // keep the Max Score during matrixes filling
          // the condition is to keep the max score in the last line (last
          // theoretical amino acid) This is for try to keep all amino acids of
          // the initial Theoretical sequence (but can modify alignment choice)
          if(itmi == last_theoretical_peptide && (*itmj).score > m_maxScore)
            {
              m_maxScore = (*itmj).score;
              m_itPosMax = itmj;
            }
        }
    }
}
 
PeptideModel
SpectralAlignment::rtrim(PrecisionPtr precision_ptr)
{
  pappso::specglob::PeptideModel model;
  pappso::specglob::PeptideSpectraCsp ref_peptide = mcsp_peptideSpectrum;
 
  double precursor_mass = mcsp_experimentalSpectrum.get()->getPrecursorMass();
  std::size_t start     = 1;
  int max_score         = m_maxScore;
  pappso::specglob::PeptideSpectraCsp ref_best_peptide = ref_peptide;
 
  while(start < ref_peptide.get()->getPeptideSp().get()->size() - 5)
    {
      pappso::PeptideSp peptide_sp = std::make_shared<pappso::Peptide>(
        ref_peptide.get()->getPeptideSp().get()->getSequence().mid(start));
      std::size_t length = peptide_sp.get()->size();
      while(peptide_sp.get()->getMass() > (precursor_mass * 1.1))
        {
          peptide_sp =
            std::make_shared<pappso::Peptide>(peptide_sp.get()->getSequence().mid(0, length - 1));
          if(length > 1)
            length--;
        }
      pappso::specglob::PeptideSpectraCsp peptide_spectra =
        std::make_shared<pappso::specglob::PeptideSpectrum>(peptide_sp);
      align(peptide_spectra, mcsp_experimentalSpectrum);
 
      if(max_score < m_maxScore)
        {
          if(m_maxScore > 0)
            {
              model = buildPeptideModel();
              model.eliminateComplementaryDelta(precision_ptr);
              model.ltrimOnRemoval();
              max_score        = m_maxScore;
              ref_best_peptide = mcsp_peptideSpectrum;
            }
        }
      start++;
    }
  align(ref_best_peptide, mcsp_experimentalSpectrum);
  return model;
}
 
void
SpectralAlignment::fillMatricesWithScores(
  const boost::numeric::ublas::matrix<SpectralAlignmentDataPoint>::iterator2 &it_pos,
  const PeptideSpectrum &peptide_spectrum,
  const ExperimentalSpectrum &experimental_spectrum)
{
  // long theoIndiceI = it_pos.index1();
  // long expeIndiceJ = it_pos.index2();
  //  in first time, we set score corresponding to type of peak (initial/mirror
  //  or both) we set initial score to peak mirror or initial
  ExperimentalSpectrumDataPointType expePeakType = getExperimentalSpectrumDataPoint(it_pos).type;
  //  experimental_spectrum.at(expeIndiceJ).type;
 
  int alignScoreToAdd     = m_scoreValues.get(ScoreValueType::scoreAlignNative);
  int reAlignScoreToAdd   = m_scoreValues.get(ScoreValueType::scoreReAlignNative);
  int reAlignScoreNOToAdd = m_scoreValues.get(ScoreValueType::scoreReAlignNativeNO);
 
  switch(expePeakType)
    {
      case ExperimentalSpectrumDataPointType::symmetric:
        // this is a symmetric peak
        alignScoreToAdd     = m_scoreValues.get(ScoreValueType::scoreAlignSym);
        reAlignScoreToAdd   = m_scoreValues.get(ScoreValueType::scoreReAlignSym);
        reAlignScoreNOToAdd = m_scoreValues.get(ScoreValueType::scoreReAlignSymNO);
        if(m_BETTER_END_RA && it_pos.index1() == peptide_spectrum.getPeptideSp().get()->size())
          {
            reAlignScoreToAdd = m_scoreValues.get(ScoreValueType::scoreNonAlign) +
                                m_scoreValues.get(ScoreValueType::scoreAlignSym);
          }
        break;
      case ExperimentalSpectrumDataPointType::both:
        // this is a native peak with a symmetric corresponding
        alignScoreToAdd     = m_scoreValues.get(ScoreValueType::scoreAlignBoth);
        reAlignScoreToAdd   = m_scoreValues.get(ScoreValueType::scoreReAlignBoth);
        reAlignScoreNOToAdd = m_scoreValues.get(ScoreValueType::scoreReAlignBothNO);
        if(m_BETTER_END_RA && it_pos.index1() == peptide_spectrum.getPeptideSp().get()->size())
          {
            reAlignScoreToAdd = m_scoreValues.get(ScoreValueType::scoreNonAlign) +
                                m_scoreValues.get(ScoreValueType::scoreAlignBoth);
          }
        break;
      default:
        // this is a native or synthetic peak
        alignScoreToAdd     = m_scoreValues.get(ScoreValueType::scoreAlignNative);
        reAlignScoreToAdd   = m_scoreValues.get(ScoreValueType::scoreReAlignNative);
        reAlignScoreNOToAdd = m_scoreValues.get(ScoreValueType::scoreReAlignNativeNO);
        if(m_BETTER_END_RA && it_pos.index1() == peptide_spectrum.getPeptideSp().get()->size())
          {
            reAlignScoreToAdd = m_scoreValues.get(ScoreValueType::scoreNonAlign) +
                                m_scoreValues.get(ScoreValueType::scoreAlignNative);
          }
        break;
    }
  SpectralAlignmentDataPoint &matrix_data_point_i_j = *it_pos;
 
  /*
  long k = getkValue(it_pos,
                     peptide_spectrum.at(it_pos.index1()).diff_mz,
                     experimental_spectrum);
 
  if(it_pos.index1() == 15)
    qDebug() << "position1" << it_pos.index1() << "position2" << it_pos.index1()
  << " k=" << k << " " << peptide_spectrum.at(it_pos.index1()).diff_mz;
  */
  pappso::MzRange aaMassRange(peptide_spectrum.at(it_pos.index1()).diff_mz, m_precisionPtr);
  qDebug();
  // long k = -1;
  auto itKpeak = experimental_spectrum.reverseFindDiffMz(it_pos.index2(), aaMassRange);
  // if(itKpeak != experimental_spectrum.rend())
  // k = itKpeak->indice;
 
  if(itKpeak == experimental_spectrum.rend())
    {
      matrix_data_point_i_j.score = m_matrix(it_pos.index1() - 1, it_pos.index2()).score +
                                    (m_scoreValues.get(ScoreValueType::scoreNonAlign));
      matrix_data_point_i_j.origin_column_indices = it_pos.index2();
      matrix_data_point_i_j.alignment_type        = SpectralAlignmentType::nonAlign;
    }
  else
    {
      SpectralAlignmentDataPoint &matrix_data_point_previ_k =
        m_matrix(it_pos.index1() - 1, itKpeak->indice);
      int scoreAlignK = matrix_data_point_previ_k.score + alignScoreToAdd;
      // if it come from non align, we must verify that there is no offset from
      // previous align
      if(matrix_data_point_previ_k.alignment_type == SpectralAlignmentType::nonAlign)
        {
          int l;
          for(l = it_pos.index1() - 1; l > 0; l--)
            {
              if(m_matrix(l, itKpeak->indice).origin_column_indices != itKpeak->indice)
                break;
            }
          if(std::abs(m_matrix(l, itKpeak->indice).mass_difference - (*it_pos).mass_difference) >
             m_precisionPtr->getNominal())
            scoreAlignK = matrix_data_point_previ_k.score + reAlignScoreToAdd;
        }
 
      // int[0] = the j value m and int[1] = the score value
      SpectralAlignmentDataPoint reAlignBestScore =
        getBestRealignScore(it_pos, itKpeak->indice, reAlignScoreToAdd, reAlignScoreNOToAdd);
 
      // For debug to see value for any match
      // System.out.println("score k = " + scoreAlignK + " - score m = " +
      // reAlignBestScore[1] + " - origin m = "
      // + reAlignBestScore[0]);
 
      if(scoreAlignK >= reAlignBestScore.score)
        {
          // setMatricesData(theoIndiceI, expeIndiceJ, scoreAlignK, k, 2);
          matrix_data_point_i_j.score                 = scoreAlignK;
          matrix_data_point_i_j.origin_column_indices = itKpeak->indice;
          matrix_data_point_i_j.alignment_type        = SpectralAlignmentType::align;
        }
      else
        {
          /*setMatricesData(theoIndiceI,
                          expeIndiceJ,
                          reAlignBestScore[1],
                          reAlignBestScore[0],
                          reAlignBestScore[2]);*/
 
          matrix_data_point_i_j = reAlignBestScore;
        }
    }
}
 
 
SpectralAlignmentDataPoint
SpectralAlignment::getBestRealignScore(
  const boost::numeric::ublas::matrix<SpectralAlignmentDataPoint>::iterator2 &it_pos,
  std::size_t expeIndicesK,
  int reAlignScore,
  int alignScoreToAdd)
{
 
  std::size_t previous_peptide_row             = it_pos.index1() - 1;
  SpectralAlignmentDataPoint return_data_point = *it_pos;
  int bestScore                                = -10000;
  int origin                                   = -1;
 
  // find the best score column indice on previous row, walking back from
  // expeIndicesK m is a j value between 0 and k where a realign can be do if we
  // accept mass offset
  for(long m = expeIndicesK; m > -1; m--)
    {
      // the >= here is for keep the highest S value if there is multiple S with
      // the same best score
      if(m_matrix(previous_peptide_row, m).score > bestScore)
        {
          bestScore = m_matrix(previous_peptide_row, m).score;
          origin    = m;
        }
    }
 
  return_data_point.origin_column_indices = origin;
  return_data_point.score                 = bestScore + reAlignScore;
  return_data_point.alignment_type        = SpectralAlignmentType::reAlign; // [2] = 1;
 
  if(origin == -1)
    return return_data_point;
 
  // We check for the last alignment if we have chain of Non Alignment to
  // compare the last mass offset found
  std::size_t lastAlignIndiceI = previous_peptide_row;
  for(long l = previous_peptide_row; l > 0; l--)
    {
      if(m_matrix(l, origin).origin_column_indices != 0)
        {
          lastAlignIndiceI = l;
          break;
        }
    }
 
  // if the difference of mass offset between actual state and last align (or
  // realign) is null, we consider that to an alignment
  if((lastAlignIndiceI != (previous_peptide_row)) && (it_pos.index1() > 1) &&
     (std::abs(m_matrix(previous_peptide_row, expeIndicesK).mass_difference -
               m_matrix(lastAlignIndiceI, origin).mass_difference) < m_precisionPtr->getNominal()))
    {
      return_data_point.score          = bestScore + alignScoreToAdd;
      return_data_point.alignment_type = SpectralAlignmentType::align; //[2] = 2;
    }
 
  // we return the origin (value of m) and the associate calculated score and
  // the type of Alignment
  return return_data_point;
}
 
 
void
SpectralAlignment::fillMassDelta(const PeptideSpectrum &peptide_spectrum,
                                 const ExperimentalSpectrum &experimental_spectrum)
{
 
  auto it_peptide  = peptide_spectrum.begin();
  auto it_spectrum = experimental_spectrum.begin();
  for(auto itr1 = m_matrix.begin1(); itr1 != m_matrix.end1(); ++itr1, it_peptide++)
    {
      it_spectrum = experimental_spectrum.begin();
      for(auto itr2 = itr1.begin(); itr2 != itr1.end(); itr2++, it_spectrum++)
        {
          (*itr2).alignment_type = SpectralAlignmentType::nonAlign;
          if(it_peptide == peptide_spectrum.begin())
            (*itr2).alignment_type = SpectralAlignmentType::nonAlign;
          (*itr2).origin_column_indices = 0;
          (*itr2).score                 = 0;
          (*itr2).mass_difference       = it_spectrum->peak_mz - it_peptide->mz;
          // qDebug() << " i=" << itr2.index1() << " j=" << itr2.index2()
          //         << " mass_difference=" << (*itr2).mass_difference;
        }
    }
}
 
const matrix<SpectralAlignmentDataPoint> &
SpectralAlignment::getMatrix() const
{
  return m_matrix;
}
 
std::vector<int>
SpectralAlignment::getScoreRow(std::size_t row_indice) const
{
  std::vector<int> score;
 
  auto itr1 = m_matrix.begin1() + row_indice;
  for(auto itr2 = itr1.begin(); itr2 != itr1.end(); itr2++)
    {
      score.push_back((*itr2).score);
    }
 
  return score;
}
 
int
SpectralAlignment::getMaxScore() const
{
  return m_maxScore;
}
 
boost::numeric::ublas::matrix<SpectralAlignmentDataPoint>::iterator2
SpectralAlignment::getMaxPosIterator() const
{
  return m_itPosMax;
}
 
const ExperimentalSpectrumDataPoint &
SpectralAlignment::getExperimentalSpectrumDataPoint(
  const boost::numeric::ublas::matrix<SpectralAlignmentDataPoint>::iterator2 &itpos) const
{
  return mcsp_experimentalSpectrum.get()->at(itpos.index2());
}
 
double
SpectralAlignment::getPrecursorMzDelta() const
{
  return m_precursorMassDelta;
}
 
 
QString
SpectralAlignment::backTrack() const
{
  if(m_maxScore < 1)
    {
      throw pappso::PappsoException("no backtrack");
    }
  QString pepModified  = "";
  QString theoSequence = mcsp_peptideSpectrum.get()->getPeptideSp().get()->getSequence();
  int actualI          = m_itPosMax.index1();
  int prevI;
  int actualJ = m_itPosMax.index2();
  int prevJ;
  double actualDelMass;
  double prevDelMass;
  int modifCount        = 0;
  double totExplainMass = 0.0;
 
  // System.out.println(theoSequence);
  // System.out.println("actualI is : " + actualI + " and actualJ is : " +
  // actualJ);
 
  while(actualI > 0)
    {
      // System.out.println(actualI);
      // define the actual mass delta to the actual i and j
      QString tempPepSeq = "";
      QString tempAA     = "";
      QString aminoAcid  = "";
      actualDelMass      = m_matrix(actualI, actualJ).mass_difference;
      prevI              = actualI - 1;
      prevJ              = m_matrix(actualI, actualJ).origin_column_indices;
      prevDelMass        = m_matrix(prevI, prevJ).mass_difference;
 
      // We checking first if last Amino acid are not aligned
 
      if(m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::nonAlign)
        {
          // System.out.println("I'm NON ALIGN");
          tempPepSeq = QString("[%1]").arg(theoSequence.at(actualI - 1));
 
          while(m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign &&
                prevI > 0)
            {
              tempPepSeq = QString("[%1]").arg(theoSequence.at(prevI - 1)) + tempPepSeq;
              prevI--;
            }
          // modifCount++;
          actualI = prevI;
          actualJ = prevJ;
 
          pepModified = tempPepSeq + pepModified;
 
          qDebug() << "a1 pepModified=" << pepModified;
          // if not, we are in the case where there is an alignment or a
          // realignment if there is Align or a re-align, put the letter of
          // Amino Acid and check what we get before to see if we have a mass
          // change
        }
      else
        {
          // we put the actual amino acid because he is found, and in function
          // of alignment type, we can have a mass offset
          tempPepSeq = QString("%1").arg(theoSequence.at(actualI - 1));
          aminoAcid  = tempPepSeq;
 
          // we check if there is a deletion before the founded aminoAcid
          if(prevI > 0 && m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign)
            {
 
              // modifCount++;
 
              // we continue to check if this is a deletion bloc to keep the
              // mass difference due to the deletion of the block
              while(prevI > 0 &&
                    m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign)
                {
 
                  tempPepSeq = QString("[%1]").arg(theoSequence.at(prevI - 1)) + tempPepSeq;
                  tempAA     = QString("%1").arg(theoSequence.at(prevI - 1)) + tempAA;
 
                  prevI--;
                  prevDelMass = m_matrix(prevI, prevJ).mass_difference;
                  if(prevI == 0)
                    {
                      prevDelMass = 0.0;
                    }
                }
 
              // check the mass delta to avoid showing a null mass delta
              if(std::abs(actualDelMass - prevDelMass) > m_precisionPtr->getNominal())
                {
 
                  tempPepSeq = tempPepSeq.mid(0, tempPepSeq.length() - 1);
 
                  tempPepSeq += QString("[%1]").arg(actualDelMass - prevDelMass) + aminoAcid;
                  totExplainMass += (actualDelMass - prevDelMass);
                  modifCount++;
                  qDebug() << "a2a1 tempPepSeq=" << tempPepSeq;
                }
 
              // if there this is just a re-align, we need to indicate the mass
              // offset
            }
          else if(m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::reAlign)
            {
 
              tempPepSeq = QString("[%1]").arg(actualDelMass - prevDelMass) + tempPepSeq;
              modifCount++;
              totExplainMass += (actualDelMass - prevDelMass);
 
              // the fact when you align the first amino acid, but there is a
              // leak of Amino acid in OMS solution before
            }
          else if(actualI == 1 &&
                  m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::align &&
                  (std::abs(actualDelMass) > m_precisionPtr->getNominal()))
            {
              tempPepSeq = QString("[%1]").arg(actualDelMass) + tempPepSeq;
              totExplainMass += actualDelMass;
              modifCount++;
            }
 
          pepModified = tempPepSeq + pepModified;
          qDebug() << "a2 pepModified=" << pepModified;
          // System.out.println(pepModified);
 
          actualI = prevI;
          actualJ = prevJ;
        }
    }
 
  // setModificationNumber(modifCount);
 
  return QString("%1_[%2]").arg(pepModified).arg(m_precursorMassDelta - totExplainMass);
}
 
 
PeptideModel
SpectralAlignment::buildPeptideModel() const
{
  if(m_maxScore < 1)
    {
      throw pappso::PappsoException(
        QObject::tr("building peptide model failed m_maxScore == %1").arg(m_maxScore));
    }
  PeptideModel sg_peptide_model(mcsp_experimentalSpectrum.get()->getQualifiedMassSpectrum(),
                                *(mcsp_peptideSpectrum.get()->getPeptideSp().get()));
  int actualI = m_itPosMax.index1();
  int prevI;
  int actualJ = m_itPosMax.index2();
  int prevJ;
  double actualDelMass;
  double prevDelMass;
 
  // System.out.println(theoSequence);
  // System.out.println("actualI is : " + actualI + " and actualJ is : " +
  // actualJ);
 
  while(actualI > 0)
    {
      // System.out.println(actualI);
      // define the actual mass delta to the actual i and j
      QString aminoAcid = "";
      actualDelMass     = m_matrix(actualI, actualJ).mass_difference;
      prevI             = actualI - 1;
      prevJ             = m_matrix(actualI, actualJ).origin_column_indices;
      prevDelMass       = m_matrix(prevI, prevJ).mass_difference;
 
      // We checking first if last Amino acid are not aligned
 
      if(m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::nonAlign)
        {
          // System.out.println("I'm NON ALIGN");
 
          while(m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign &&
                prevI > 0)
            {
              if(prevI > 0)
                sg_peptide_model[prevI - 1].bracket = true;
              prevI--;
            }
          // modifCount++;
          actualI = prevI;
          actualJ = prevJ;
 
          // if not, we are in the case where there is an alignment or a
          // realignment if there is Align or a re-align, put the letter of
          // Amino Acid and check what we get before to see if we have a mass
          // change
        }
      else
        {
          // we put the actual amino acid because he is found, and in function
          // of alignment type, we can have a mass offset
 
          // we check if there is a deletion before the founded aminoAcid
          if(prevI > 0 && m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign)
            {
 
              // modifCount++;
 
              // we continue to check if this is a deletion bloc to keep the
              // mass difference due to the deletion of the block
              while(prevI > 0 &&
                    m_matrix(prevI, prevJ).alignment_type == SpectralAlignmentType::nonAlign)
                {
 
 
                  if(prevI > 0)
                    sg_peptide_model[prevI - 1].bracket = true;
                  prevI--;
                  prevDelMass = m_matrix(prevI, prevJ).mass_difference;
                  if(prevI == 0)
                    {
                      prevDelMass = 0.0;
                    }
                }
 
              // check the mass delta to avoid showing a null mass delta
              if(std::abs(actualDelMass - prevDelMass) > m_precisionPtr->getNominal())
                {
 
                  int mass_i = actualI - 1;
                  if(mass_i == 0)
                    {
                      sg_peptide_model.setBeginMassDelta(actualDelMass - prevDelMass);
                    }
                  else
                    {
                      sg_peptide_model[mass_i - 1].mass_difference = actualDelMass - prevDelMass;
                    }
                }
 
              // if there this is just a re-align, we need to indicate the mass
              // offset
            }
          else if(m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::reAlign)
            {
 
              int mass_i = actualI - 1;
              if(mass_i == 0)
                {
                  sg_peptide_model.setBeginMassDelta(actualDelMass - prevDelMass);
                }
              else
                {
                  sg_peptide_model[mass_i - 1].mass_difference = actualDelMass - prevDelMass;
                }
              // the fact when you align the first amino acid, but there is a
              // leak of Amino acid in OMS solution before
            }
          else if(actualI == 1 &&
                  m_matrix(actualI, actualJ).alignment_type == SpectralAlignmentType::align &&
                  (std::abs(actualDelMass) > m_precisionPtr->getNominal()))
            {
              sg_peptide_model.setBeginMassDelta(actualDelMass);
            }
 
          actualI = prevI;
          actualJ = prevJ;
        }
    }
 
  // setModificationNumber(modifCount);
 
  // return QString("%1_[%2]")
  //  .arg(pepModified)
  //   .arg(m_precursorMassDelta - totExplainMass);
  return sg_peptide_model;
}
} // namespace specglob
} // namespace pappso