semiglobalalignment.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specpeptidoms/semiglobalalignment.cpp
 * \date 24/03/2025
 * \author Aurélien Berthier
 * \brief protein to spectrum alignment
 *
 * 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 <QString>
#include <QStringRef>
#include "semiglobalalignment.h"
#include "types.h"
#include "../../amino_acid/aacode.h"
#include "../../protein/protein.h"
#include "correctiontree.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
 
void
pappso::specpeptidoms::Alignment::reset()
{
  peaks.clear();
  interpretation = "";
  score          = 0;
  begin_shift    = 0.0;
  end_shift      = 0.0;
  shifts.clear();
  SPC       = 0;
  beginning = 0;
  end       = 0;
}
 
 
pappso::specpeptidoms::SemiGlobalAlignment::SemiGlobalAlignment(
  const ScoreValues &score_values,
  const pappso::PrecisionPtr precision_ptr,
  const pappso::AaCode &aaCode)
  : m_scorevalues(score_values), m_aaCode(aaCode)
{
  m_precision_ptr = precision_ptr;
  m_interest_cells.push_back({0, 0, 0, 0});
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::fastAlign(const SpOMSSpectrum &spectrum,
                                                      const QString &protein_seq_in,
                                                      const QString &protein_id)
{
  // TODO don't forget to reset any important variable
  m_best_alignment.reset();
  m_best_corrected_alignment.reset();
  m_best_post_processed_alignment.reset();
 
  std::size_t sequence_length = protein_seq_in.size();
  QString protein_seq(protein_seq_in);
  std::reverse(protein_seq.begin(), protein_seq.end());
  pappso::Enums::AminoAcidChar aa;
  std::vector<AaPosition> aa_positions;
  m_interest_cells.at(0).n_row     = 0;
  m_interest_cells.at(0).score     = 0;
  m_interest_cells.at(0).beginning = 0;
  m_interest_cells.at(0).tree_id   = 0;
 
  for(std::size_t i = 1; i < m_interest_cells.size(); i++)
    {
      m_interest_cells.at(i).n_row     = 0;
      m_interest_cells.at(i).score     = m_scorevalues.get(ScoreType::init);
      m_interest_cells.at(i).beginning = 0;
      m_interest_cells.at(i).tree_id   = 0;
    }
 
  for(std::size_t iter = m_interest_cells.size(); iter < spectrum.size(); iter++)
    {
      m_interest_cells.push_back({0, m_scorevalues.get(ScoreType::init), 0, 0});
    }
  m_location_saver.resetLocationSaver();
  for(std::size_t row_number = 1; row_number <= sequence_length; row_number++)
    {
      if(protein_seq[row_number - 1] == 'U' or protein_seq[row_number - 1] == 'X')
        {
          continue;
        }
      try
        {
          aa = pappso::Enums::AminoAcidChar(protein_seq[row_number - 1].toLatin1());
        }
      catch(const pappso::PappsoException &err)
        {
          throw err;
        }
      aa_positions = spectrum.getAaPositions(aa);
      updateAlignmentMatrix(protein_seq, row_number, aa_positions, spectrum, true, protein_id);
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::preciseAlign(const SpOMSSpectrum &spectrum,
                                                         const QString &protein_seq,
                                                         const QString &protein_id,
                                                         const std::size_t beginning,
                                                         const std::size_t length)
{
  qDebug();
  std::size_t length2;
  if((qsizetype)(beginning + length) <= protein_seq.size())
    {
      length2 = length;
    }
  else
    {
      length2 = protein_seq.size() - beginning;
    }
 
  qDebug();
  QString sequence = protein_seq.sliced(protein_seq.size() - beginning - length2, length2);
  std::reverse(sequence.begin(), sequence.end());
  std::vector<AaPosition> aa_positions;
  CorrectionTree correction_tree;
 
  qDebug();
  m_scenario.reserve(length2 + 1, spectrum.size());
  m_interest_cells.reserve(spectrum.size());
  m_interest_cells.at(0).n_row     = 0;
  m_interest_cells.at(0).score     = 0;
  m_interest_cells.at(0).beginning = 0;
  m_interest_cells.at(0).tree_id   = 0;
  for(std::size_t i = 1; i < m_interest_cells.size(); i++)
    {
      m_interest_cells.at(i).n_row     = 0;
      m_interest_cells.at(i).score     = m_scorevalues.get(ScoreType::init);
      m_interest_cells.at(i).beginning = 0;
      m_interest_cells.at(i).tree_id   = 0;
    }
  qDebug();
  for(std::size_t iter = m_interest_cells.size(); iter < spectrum.size(); iter++)
    {
      m_interest_cells.push_back({0, m_scorevalues.get(ScoreType::init), 0, 0});
    }
  qDebug();
  m_scenario.resetScenario();
  qDebug();
  for(std::size_t row_number = 1; row_number <= length2; row_number++)
    {
      qDebug() << "row_number - 1=" << row_number - 1 << " sequence.size()=" << sequence.size();
      if(sequence[row_number - 1] == 'U' or sequence[row_number - 1] == 'X')
        {
          continue;
        }
 
      qDebug() << "row_number - 1=" << row_number - 1 << " sequence.size()=" << sequence.size();
      // aa           = Aa(sequence[row_number - 1].unicode());
      aa_positions =
        spectrum.getAaPositions(pappso::Enums::AminoAcidChar(sequence[row_number - 1].toLatin1()));
      updateAlignmentMatrix(sequence, row_number, aa_positions, spectrum, false, protein_id);
    }
  qDebug();
  saveBestAlignment(sequence, spectrum, protein_seq.size() - beginning);
 
  qDebug() << m_scenario.getBestScore() << " " << MIN_ALIGNMENT_SCORE;
  // Correction : if complementary peaks are used, corrected spectra without one of the two peaks
  // are generated and aligned. The best alignment is kept.
  if(m_scenario.getBestScore() >
     MIN_ALIGNMENT_SCORE) // We only correct alignments with acceptable scores
    {
      qDebug();
      m_best_corrected_alignment.score = 0;
      for(std::size_t iter : m_best_alignment.peaks)
        {
          if(iter > spectrum.getComplementaryPeak(iter))
            {
              break;
            }
          else if(std::find(m_best_alignment.peaks.begin(),
                            m_best_alignment.peaks.end(),
                            spectrum.getComplementaryPeak(iter)) != m_best_alignment.peaks.end())
            {
              correction_tree.addPeaks(iter, spectrum.getComplementaryPeak(iter));
            }
        }
      qDebug();
      std::vector<std::vector<std::size_t>> corrections = correction_tree.getPeaks();
      if(corrections.size() > 0)
        {
          m_best_alignment.score = 0; //  Reset the best alignment score (we dont want to keep the
                                      //  original alignment if corrections are needed)
          qDebug();
          for(auto peaks_to_remove : corrections)
            {
              qDebug();
              correctAlign(
                1, sequence, protein_id, spectrum, peaks_to_remove, protein_seq.size() - beginning);
              qDebug();
            }
          qDebug();
          m_best_alignment = m_best_corrected_alignment;
        }
    }
  qDebug();
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::correctAlign(int recursive_call_count,
                                                         const QString &sequence,
                                                         const QString &protein_id,
                                                         const SpOMSSpectrum &spectrum,
                                                         std::vector<std::size_t> peaks_to_remove,
                                                         std::size_t offset)
{
  qDebug() << recursive_call_count;
  std::vector<AaPosition> aa_positions;
  CorrectionTree correction_tree;
 
  m_interest_cells.at(0).n_row     = 0;
  m_interest_cells.at(0).score     = 0;
  m_interest_cells.at(0).beginning = 0;
  m_interest_cells.at(0).tree_id   = 0;
  for(std::size_t i = 1; i < m_interest_cells.size(); i++)
    {
      m_interest_cells.at(i).n_row     = 0;
      m_interest_cells.at(i).score     = m_scorevalues.get(ScoreType::init);
      m_interest_cells.at(i).beginning = 0;
      m_interest_cells.at(i).tree_id   = 0;
    }
 
  m_scenario.resetScenario();
  qDebug();
  for(qsizetype row_number = 1; row_number <= sequence.size(); row_number++)
    {
      qDebug() << row_number - 1 << " " << sequence.size();
      if(sequence[row_number - 1] == 'U' or sequence[row_number - 1] == 'X')
        {
          continue;
        }
      pappso::Enums::AminoAcidChar aa =
        pappso::Enums::AminoAcidChar(sequence[row_number - 1].toLatin1());
      qDebug();
      aa_positions = spectrum.getAaPositions(aa, peaks_to_remove);
      qDebug();
      updateAlignmentMatrix(sequence, row_number, aa_positions, spectrum, false, protein_id);
      qDebug();
    }
 
  qDebug();
  // Correction : if complementary peaks are used, corrected spectra without one of the two peaks
  // are generated and aligned. The best alignment is kept.
  qDebug() << m_scenario.getBestScore();
  if(m_scenario.getBestScore() >
     MIN_ALIGNMENT_SCORE) // We only correct alignments with acceptable scores
    {
      qDebug();
      qDebug() << sequence << " " << recursive_call_count;
      qDebug() << offset;
      qDebug() << spectrum.getPrecursorCharge();
      saveBestAlignment(sequence, spectrum, offset);
      qDebug();
      for(std::size_t iter : m_best_alignment.peaks)
        {
          if(iter > spectrum.getComplementaryPeak(iter))
            {
              break;
            }
          else if(std::find(m_best_alignment.peaks.begin(),
                            m_best_alignment.peaks.end(),
                            spectrum.getComplementaryPeak(iter)) != m_best_alignment.peaks.end())
            {
              correction_tree.addPeaks(iter, spectrum.getComplementaryPeak(iter));
            }
        }
      std::vector<std::vector<std::size_t>> corrections = correction_tree.getPeaks();
      if(corrections.size() > 0)
        {
          for(auto peaks_to_remove : corrections)
            {
              correctAlign(
                recursive_call_count + 1, sequence, protein_id, spectrum, peaks_to_remove, offset);
            }
        }
      else if(m_scenario.getBestScore() > m_best_corrected_alignment.score)
        {
          m_best_corrected_alignment = m_best_alignment;
        }
    }
  qDebug();
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::postProcessingAlign(const SpOMSSpectrum &spectrum,
                                                                const QString &protein_seq,
                                                                const QString &protein_id,
                                                                std::size_t beginning,
                                                                std::size_t length,
                                                                const std::vector<double> &shifts)
{
  std::size_t current_SPC               = m_best_alignment.SPC;
  int current_best_score                = m_best_alignment.score;
  m_best_post_processed_alignment.score = 0;
  for(double precursor_mass_error : shifts)
    {
      SpOMSSpectrum corrected_spectrum(spectrum, precursor_mass_error);
      preciseAlign(corrected_spectrum, protein_seq, protein_id, beginning, length);
      if(m_best_alignment.score >= m_best_post_processed_alignment.score)
        {
          m_best_post_processed_alignment = m_best_alignment;
        }
    }
  if(m_best_post_processed_alignment.SPC > current_SPC &&
     m_best_post_processed_alignment.score >= current_best_score)
    {
      m_best_alignment = m_best_post_processed_alignment;
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::updateAlignmentMatrix(
  const QString &sequence,
  const std::size_t row_number,
  const std::vector<AaPosition> aa_positions,
  const SpOMSSpectrum &spectrum,
  const bool fast_align,
  const QString &protein)
{
  int score_found, score_shift, best_score, alt_score, tree_id;
  uint32_t condition; // FIXME : may be used uninitialised
  std::size_t best_column, shift, beginning, missing_aas, length, perfect_shift_origin;
  KeyCell *current_cell_ptr, *tested_cell_ptr;
  AlignType alignment_type, temp_align_type;
 
  m_updated_cells.reserve(aa_positions.size());
 
  // Computation of the threePeaks condition, see spomsspectrum.h for more details.
  if(fast_align)
    {
      condition = 3;
      if(row_number > 1)
        {
          condition += 2 << m_aaCode.getAaCode(sequence[row_number - 2].unicode());
          qDebug() << "condition" << condition << sequence[row_number - 2]
                   << m_aaCode.getAaCode(sequence[row_number - 2].unicode());
        }
    }
 
  for(std::vector<AaPosition>::const_iterator aa_position = aa_positions.begin();
      aa_position != aa_positions.end();
      aa_position++)
    {
      if(((condition & aa_position->condition) != 0) ||
         !fast_align) // Verification of the threePeaks condition (only during first alignment).
        {
          current_cell_ptr = &m_interest_cells.at(aa_position->r_peak);
          if(spectrum.peakType(aa_position->r_peak) ==
             specglob::ExperimentalSpectrumDataPointType::both)
            {
              score_found = m_scorevalues.get(ScoreType::foundDouble);
              score_shift = m_scorevalues.get(ScoreType::foundShiftDouble);
            }
          else
            {
              score_found = m_scorevalues.get(ScoreType::found);
              score_shift = m_scorevalues.get(ScoreType::foundShift);
            }
 
          // not found case (always computed)
          best_column = aa_position->r_peak;
          best_score  = current_cell_ptr->score + (row_number - current_cell_ptr->n_row) *
                                                   m_scorevalues.get(ScoreType::notFound);
          beginning      = current_cell_ptr->beginning;
          tree_id        = current_cell_ptr->tree_id;
          alignment_type = AlignType::notFound;
 
          // found case (Can only happen if the left peak is supported)
          if(aa_position->l_support)
            {
              tested_cell_ptr = &m_interest_cells.at(aa_position->l_peak);
              if(aa_position->l_peak == 0)
                {
                  alt_score = tested_cell_ptr->score + score_found;
                }
              else
                {
                  if(tested_cell_ptr->n_row == row_number - 1)
                    {
                      alt_score = tested_cell_ptr->score +
                                  (row_number - tested_cell_ptr->n_row - 1) *
                                    m_scorevalues.get(ScoreType::notFound) +
                                  score_found;
                    }
                  else
                    {
                      alt_score = tested_cell_ptr->score +
                                  (row_number - tested_cell_ptr->n_row - 1) *
                                    m_scorevalues.get(ScoreType::notFound) +
                                  score_shift;
                    }
                }
              if(alt_score >= best_score)
                {
                  alignment_type = AlignType::found;
                  best_score     = alt_score;
                  best_column    = aa_position->l_peak;
                  if(best_column == 0)
                    {
                      if(row_number < ALIGNMENT_SURPLUS)
                        {
                          beginning = 0;
                        }
                      else
                        {
                          beginning =
                            std::max((std::size_t)(row_number - ALIGNMENT_SURPLUS), (std::size_t)0);
                        }
                      if(fast_align)
                        {
                          tree_id = m_location_saver.getNextTree();
                        }
                    }
                  else
                    {
                      beginning = tested_cell_ptr->beginning;
                      tree_id   = tested_cell_ptr->tree_id;
                    }
                }
            }
 
          // generic shift case (all shifts are tested)
          shift = 0;
          while(shift < aa_position->next_l_peak)
            {
              tested_cell_ptr = &m_interest_cells.at(aa_position->next_l_peak - shift);
              // verification saut parfait
              if(perfectShiftPossible(sequence,
                                      spectrum,
                                      tested_cell_ptr->n_row,
                                      row_number,
                                      aa_position->next_l_peak - shift,
                                      aa_position->r_peak))
                {
                  alt_score = tested_cell_ptr->score +
                              (row_number - tested_cell_ptr->n_row - 1) *
                                m_scorevalues.get(ScoreType::notFound) +
                              score_found;
                  temp_align_type = AlignType::perfectShift;
                }
              else
                {
                  alt_score = tested_cell_ptr->score +
                              (row_number - tested_cell_ptr->n_row - 1) *
                                m_scorevalues.get(ScoreType::notFound) +
                              score_shift;
                  temp_align_type = AlignType::shift;
                }
              if(alt_score > best_score)
                {
                  alignment_type = temp_align_type;
                  best_score     = alt_score;
                  best_column    = aa_position->next_l_peak - shift;
                  beginning      = tested_cell_ptr->beginning;
                  tree_id        = tested_cell_ptr->tree_id;
                }
              shift++;
            }
 
          // case shift from column 0 (no penalties if all precedent amino acids are missed)
          tested_cell_ptr = &m_interest_cells.at(0);
          // verification saut parfait
          perfect_shift_origin =
            perfectShiftPossibleFrom0(sequence, spectrum, row_number, aa_position->r_peak);
          if(perfect_shift_origin != row_number)
            {
              alt_score       = tested_cell_ptr->score + score_found;
              temp_align_type = AlignType::perfectShift;
            }
          else
            {
              alt_score       = tested_cell_ptr->score + score_shift;
              temp_align_type = AlignType::shift;
            }
          if(alt_score > best_score)
            {
              alignment_type = temp_align_type;
              best_score     = alt_score;
              best_column    = 0;
              missing_aas =
                std::floor(spectrum.getMZShift(0, aa_position->l_peak) / m_aaCode.getMass('G'));
              if(row_number < ALIGNMENT_SURPLUS + missing_aas)
                {
                  beginning = 0;
                }
              else
                {
                  beginning = std::max((std::size_t)(row_number - missing_aas - ALIGNMENT_SURPLUS),
                                       (std::size_t)0);
                }
              if(fast_align)
                {
                  tree_id = m_location_saver.getNextTree();
                }
            }
          if(best_column != aa_position->r_peak)
            {
              m_updated_cells.push_back(
                {aa_position->r_peak, {row_number, best_score, beginning, tree_id}});
            }
          if(best_score > m_location_saver.getMinScore(tree_id) && fast_align)
            {
              length =
                row_number - beginning + 1 +
                std::ceil(spectrum.getMissingMass(aa_position->r_peak) / m_aaCode.getMass('G')) +
                ALIGNMENT_SURPLUS;
              m_location_saver.addLocation(beginning, length, tree_id, best_score, protein);
            }
          else if(!fast_align)
            {
              if(alignment_type == AlignType::perfectShift && best_column == 0)
                {
                  m_scenario.saveOrigin(row_number,
                                        aa_position->r_peak,
                                        perfect_shift_origin,
                                        0,
                                        best_score,
                                        AlignType::perfectShift);
                }
              else
                {
                  m_scenario.saveOrigin(row_number,
                                        aa_position->r_peak,
                                        m_interest_cells.at(best_column).n_row,
                                        best_column,
                                        best_score,
                                        alignment_type);
                }
            }
        }
    }
 
  //  Update row number in column 0
  m_updated_cells.push_back({0, {row_number, 0, 0, 0}});
 
  //  Save updated key cells in the matrix
  while(m_updated_cells.size() > 0)
    {
      m_interest_cells.at(m_updated_cells.back().first) = m_updated_cells.back().second;
      m_updated_cells.pop_back();
    }
}
 
bool
pappso::specpeptidoms::SemiGlobalAlignment::perfectShiftPossible(const QString &sequence,
                                                                 const SpOMSSpectrum &spectrum,
                                                                 const std::size_t origin_row,
                                                                 const std::size_t current_row,
                                                                 const std::size_t l_peak,
                                                                 const std::size_t r_peak) const
{
  double missing_mass = 0;
  for(QString::const_iterator iter = sequence.begin() + origin_row;
      iter != sequence.begin() + current_row;
      iter++)
    {
      missing_mass += m_aaCode.getMass(iter->toLatin1()); // Aa(iter->unicode()).getMass();
    }
  pappso::MzRange mz_range(missing_mass, m_precision_ptr);
  return mz_range.contains(spectrum.getMZShift(l_peak, r_peak));
}
 
std::size_t
pappso::specpeptidoms::SemiGlobalAlignment::perfectShiftPossibleFrom0(
  const QString &sequence,
  const SpOMSSpectrum &spectrum,
  const std::size_t current_row,
  const std::size_t r_peak) const
{
  std::size_t perfect_shift_origin = current_row;
  double missing_mass              = spectrum.getMZShift(0, r_peak);
  pappso::MzRange mz_range(missing_mass, m_precision_ptr);
  double aa_mass = 0;
  while(aa_mass < missing_mass && perfect_shift_origin > 0 && !mz_range.contains(aa_mass))
    {
      aa_mass += m_aaCode.getMass(
        sequence.at(perfect_shift_origin - 1)
          .toLatin1()); // Aa(sequence.at(perfect_shift_origin - 1).unicode()).getMass();
      perfect_shift_origin--;
    }
  if(mz_range.contains(aa_mass))
    {
      return perfect_shift_origin;
    }
  else
    {
      return current_row;
    }
}
 
std::size_t
pappso::specpeptidoms::SemiGlobalAlignment::perfectShiftPossibleEnd(const QString &sequence,
                                                                    const SpOMSSpectrum &spectrum,
                                                                    std::size_t end_row,
                                                                    std::size_t end_peak) const
{
  std::size_t perfect_shift_end = end_row + 1;
  double missing_mass           = spectrum.getMissingMass(end_peak);
  pappso::MzRange mz_range(missing_mass, m_precision_ptr);
  double aa_mass = 0;
  while(aa_mass < missing_mass && perfect_shift_end < (std::size_t)sequence.size() &&
        !mz_range.contains(aa_mass))
    {
      aa_mass += m_aaCode.getMass(
        sequence.at(perfect_shift_end - 1)
          .toLatin1()); // Aa(sequence.at(perfect_shift_end - 1).unicode()).getMass();
      perfect_shift_end++;
    }
  if(mz_range.contains(aa_mass))
    {
      return perfect_shift_end - 1;
    }
  else
    {
      return end_row;
    }
}
 
pappso::specpeptidoms::SemiGlobalAlignment::~SemiGlobalAlignment()
{
}
 
pappso::specpeptidoms::LocationSaver
pappso::specpeptidoms::SemiGlobalAlignment::getLocationSaver() const
{
  return m_location_saver;
}
 
pappso::specpeptidoms::Scenario
pappso::specpeptidoms::SemiGlobalAlignment::getScenario() const
{
  return m_scenario;
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::saveBestAlignment(const QString &sequence,
                                                              const SpOMSSpectrum &spectrum,
                                                              std::size_t offset)
{
  qDebug();
  m_best_alignment.peaks.clear();
  m_best_alignment.shifts.clear();
  std::vector<QString> temp_interpretation;
  std::size_t previous_row; // FIXME : may be used uninitialised
  std::size_t previous_column = 0;
  std::size_t perfect_shift_end;
  std::pair<std::vector<ScenarioCell>, int> best_alignment = m_scenario.getBestAlignment();
  m_best_alignment.score                                   = best_alignment.second;
  QString skipped_aa;
  double skipped_mass;
  // Retrieving beginning and end
  if(best_alignment.first.front().previous_row > offset)
    {
      throw pappso::PappsoException(
        QString("best_alignment.first.front().previous_row > offset %1 %2")
          .arg(offset)
          .arg(best_alignment.first.front().previous_row));
    }
  if(best_alignment.first.back().previous_row > offset)
    {
      throw pappso::PappsoException(
        QString("best_alignment.first.back().previous_row > offset %1 %2")
          .arg(offset)
          .arg(best_alignment.first.back().previous_row));
    }
  m_best_alignment.beginning = offset - best_alignment.first.front().previous_row;
  m_best_alignment.end       = offset - best_alignment.first.back().previous_row;
 
  qDebug();
  // Filling temp_interpretation and peaks vectors
  for(auto cell : best_alignment.first)
    {
      switch(cell.alignment_type)
        {
          case AlignType::found:
            temp_interpretation.push_back(sequence.at(previous_row - 1));
            if(previous_row > cell.previous_row + 1)
              {
                skipped_mass = m_aaCode.getMass(
                  sequence.at(previous_row - 1)
                    .toLatin1()); // Aa(sequence.at(previous_row - 1).unicode()).getMass();
                skipped_aa =
                  sequence.sliced(cell.previous_row, previous_row - cell.previous_row - 1);
                for(auto aa : skipped_aa)
                  {
                    temp_interpretation.push_back('[' + aa + ']');
                    skipped_mass += m_aaCode.getMass(aa.toLatin1()); // Aa(aa.unicode()).getMass();
                  }
                temp_interpretation.push_back(
                  '[' +
                  QString::number(spectrum.getMZShift(cell.previous_column, previous_column) -
                                  skipped_mass) +
                  ']');
              }
            m_best_alignment.peaks.push_back(cell.previous_column);
            break;
          case AlignType::notFound:
            temp_interpretation.push_back('[' + sequence.at(previous_row - 1) + ']');
            break;
          case AlignType::shift:
            m_best_alignment.peaks.push_back(cell.previous_column);
            temp_interpretation.push_back(sequence.at(previous_row - 1));
            temp_interpretation.push_back(
              '[' +
              QString::number(spectrum.getMZShift(cell.previous_column, previous_column) -
                              m_aaCode.getMass(sequence.at(previous_row - 1).toLatin1())) +
              ']');
            m_best_alignment.shifts.push_back(
              spectrum.getMZShift(cell.previous_column, previous_column) -
              m_aaCode.getMass(sequence.at(previous_row - 1).toLatin1()));
            break;
          case AlignType::perfectShift:
            m_best_alignment.peaks.push_back(cell.previous_column);
            skipped_aa = sequence.sliced(cell.previous_row, previous_row - cell.previous_row);
            std::reverse(skipped_aa.begin(), skipped_aa.end());
            temp_interpretation.push_back(skipped_aa);
            break;
          case AlignType::init:
            previous_row    = cell.previous_row;
            previous_column = cell.previous_column;
            m_best_alignment.peaks.push_back(cell.previous_column);
            break;
        }
      previous_row    = cell.previous_row;
      previous_column = cell.previous_column;
    }
  std::reverse(temp_interpretation.begin(), temp_interpretation.end());
  std::reverse(m_best_alignment.peaks.begin(), m_best_alignment.peaks.end());
 
  qDebug();
  // Compute begin_shift and end_shift
  MzRange zero(0, m_precision_ptr);
  m_best_alignment.begin_shift = spectrum.getMZShift(0, m_best_alignment.peaks.front());
  m_best_alignment.end_shift   = spectrum.getMissingMass(m_best_alignment.peaks.back());
  if(zero.contains(m_best_alignment.end_shift))
    {
      m_best_alignment.end_shift = 0;
    }
 
  qDebug();
  // Computing SPC
  m_best_alignment.SPC = 0;
  for(auto peak : m_best_alignment.peaks)
    {
      switch(spectrum.at(peak).type)
        {
          case specglob::ExperimentalSpectrumDataPointType::native:
            qDebug() << peak << "native";
            m_best_alignment.SPC += 1;
            break;
          case specglob::ExperimentalSpectrumDataPointType::both:
            qDebug() << peak << "both";
            m_best_alignment.SPC += 2;
            break;
          case specglob::ExperimentalSpectrumDataPointType::synthetic:
            qDebug() << peak << "synthetic";
            break;
          case specglob::ExperimentalSpectrumDataPointType::symmetric:
            qDebug() << peak << "symmetric";
            m_best_alignment.SPC += 1;
            break;
        }
    }
 
  qDebug();
  // Final check of the end shift
  if(m_best_alignment.end_shift > 0)
    {
      perfect_shift_end = perfectShiftPossibleEnd(sequence,
                                                  spectrum,
                                                  best_alignment.first.front().previous_row,
                                                  m_best_alignment.peaks.back());
      if(perfect_shift_end != best_alignment.first.front().previous_row)
        {
          skipped_aa =
            sequence.sliced(best_alignment.first.front().previous_row,
                            perfect_shift_end - best_alignment.first.front().previous_row);
          for(auto aa = skipped_aa.begin(); aa != skipped_aa.end(); aa++)
            {
              temp_interpretation.push_back('[' + *aa + ']');
            }
          m_best_alignment.beginning = offset - perfect_shift_end;
          m_best_alignment.end_shift = 0;
        }
      else
        {
          m_best_alignment.score += m_scorevalues.get(ScoreType::foundShift);
        }
    }
 
  qDebug();
  // Writing final interpretation
  m_best_alignment.interpretation = "";
  if(m_best_alignment.end_shift > 0)
    {
      m_best_alignment.interpretation += '[' + QString::number(m_best_alignment.end_shift) + ']';
    }
  for(auto str = temp_interpretation.rbegin(); str != temp_interpretation.rend(); str++)
    {
      m_best_alignment.interpretation += *(str);
    }
  if(m_best_alignment.begin_shift > 0)
    {
      m_best_alignment.interpretation += '[' + QString::number(m_best_alignment.begin_shift) + ']';
    }
  qDebug();
}
 
const pappso::specpeptidoms::Alignment &
pappso::specpeptidoms::SemiGlobalAlignment::getBestAlignment() const
{
  return m_best_alignment;
}
 
std::vector<double>
pappso::specpeptidoms::SemiGlobalAlignment::getPotentialMassErrors(const pappso::AaCode &aa_code,
                                                                   const Alignment &alignment,
                                                                   const QString &protein_seq)
{
  // qWarning() << protein_seq;
  if(alignment.end > (std::size_t)protein_seq.size())
    {
      throw pappso::ExceptionOutOfRange(QString("alignment.end > protein_seq.size() %1 %2")
                                          .arg(alignment.end)
                                          .arg(protein_seq.size()));
    }
  std::vector<double> potential_mass_errors(alignment.shifts);
  double shift = alignment.end_shift;
  std::size_t index;
  if(alignment.beginning > 0)
    { // -1 on unsigned int makes it wrong
      index = alignment.beginning - 1;
      while(shift > 0 && index > 0)
        {
          potential_mass_errors.push_back(shift);
          // qWarning() << " shift=" << shift << " index=" << index
          //           << " letter=" << protein_seq.at(index).toLatin1();
          shift -= aa_code.getMass(
            protein_seq.at(index).toLatin1()); // Aa(protein_seq.at(index).unicode()).getMass();
          index--;
        }
    }
 
  // qWarning() << "second";
  shift = alignment.begin_shift;
  index = alignment.end + 1;
  while(shift > 0 && index < (std::size_t)protein_seq.size())
    {
      potential_mass_errors.push_back(shift);
      qWarning() << " shift=" << shift << " index=" << index
                 << " letter=" << protein_seq.at(index).toLatin1();
      shift -= aa_code.getMass(
        protein_seq.at(index).toLatin1()); // Aa(protein_seq.at(index).unicode()).getMass();
      index++;
    }
  // qWarning();
  return potential_mass_errors;
}