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"
 
pappso::specpeptidoms::SemiGlobalAlignment::SemiGlobalAlignment(
  ScoreValues &score_values, const pappso::PrecisionPtr precision_ptr, 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 SpOMSSpectrumCsp &spectrum,
                                                      const QString &protein_seq_in,
                                                      const QString &protein_id)
{
  std::size_t sequence_length = protein_seq_in.size();
  QString protein_seq(protein_seq_in);
  std::reverse(protein_seq.begin(), protein_seq.end());
  pappso::Aa aa('A');
  std::vector<AaPosition> aa_positions;
  int score_found, score_shift, best_score, current_score;
  std::size_t best_column;
  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;
        }
      aa           = Aa(protein_seq[row_number - 1].unicode());
      aa_positions = spectrum->getAaPositions(aa);
      updateAlignmentMatrix(protein_seq, row_number, aa_positions, spectrum, true, protein_id);
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::preciseAlign(const SpOMSSpectrumCsp &spectrum,
                                                         const QString &protein_seq,
                                                         const QString &protein_id,
                                                         const std::size_t beginning,
                                                         const std::size_t length)
{
  std::size_t length2;
  if((qsizetype)(beginning + length) <= protein_seq.size())
    {
      length2 = length;
    }
  else
    {
      length2 = protein_seq.size() - beginning;
    }
 
  QString sequence = protein_seq.sliced(protein_seq.size() - beginning - length2, length2);
  std::reverse(sequence.begin(), sequence.end());
  Aa aa('A');
  std::vector<AaPosition> aa_positions;
  CorrectionTree correction_tree;
 
  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;
    }
  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_scenario.resetScenario();
  for(std::size_t row_number = 1; row_number <= length2; row_number++)
    {
      if(sequence[row_number - 1] == 'U' or sequence[row_number - 1] == 'X')
        {
          continue;
        }
      aa           = Aa(sequence[row_number - 1].unicode());
      aa_positions = spectrum->getAaPositions(aa);
      updateAlignmentMatrix(sequence, row_number, aa_positions, spectrum, false, protein_id);
    }
  saveBestAlignment(sequence, spectrum, protein_seq.size() - beginning);
 
  // 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
    {
      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));
            }
        }
      std::vector<std::vector<std::size_t>> corrections = correction_tree.getPeaks();
      if(corrections.size() > 0)
        {
          std::cin.get();
          m_best_alignment.score = 0; //  Reset the best alignment score (we dont want to keep the
                                      //  original alignment if corrections are needed)
          for(auto peaks_to_remove : corrections)
            {
              correctAlign(
                sequence, protein_id, spectrum, peaks_to_remove, protein_seq.size() - beginning);
            }
          std::cin.get();
          m_best_alignment = m_best_corrected_alignment;
        }
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::correctAlign(const QString &sequence,
                                                         const QString &protein_id,
                                                         const SpOMSSpectrumCsp &spectrum,
                                                         std::vector<std::size_t> peaks_to_remove,
                                                         std::size_t offset)
{
  Aa aa('A');
  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();
  for(qsizetype row_number = 1; row_number <= sequence.size(); row_number++)
    {
      if(sequence[row_number - 1] == 'U' or sequence[row_number - 1] == 'X')
        {
          continue;
        }
      aa           = Aa(sequence[row_number - 1].unicode());
      aa_positions = spectrum->getAaPositions(aa, peaks_to_remove);
      updateAlignmentMatrix(sequence, row_number, aa_positions, spectrum, false, protein_id);
    }
 
  // 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
    {
      saveBestAlignment(sequence, spectrum, offset);
      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(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;
        }
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::postProcessingAlign(const SpOMSSpectrumCsp &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)
    {
      SpOMSSpectrumCsp corrected_spectrum =
        std::make_shared<SpOMSSpectrum>(SpOMSSpectrum(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 SpOMSSpectrumCsp &spectrum,
  const bool fast_align,
  const QString &protein)
{
  int score_found, score_shift, best_score, current_score, alt_score, tree_id;
  uint32_t condition; // FIXME : may be used uninitialised
  std::size_t best_column, alt_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) / Aa('G').getMass()) +
                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 SpOMSSpectrumCsp &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 += 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 SpOMSSpectrumCsp &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 += 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 SpOMSSpectrumCsp &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 += 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 SpOMSSpectrumCsp &spectrum,
                                                              std::size_t offset)
{
  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;
  int score_found, score_shift;
  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;
  AlignType first_align_type;
  // Retrieving beginning and end
  m_best_alignment.beginning = offset - best_alignment.first.front().previous_row;
  m_best_alignment.end       = offset - best_alignment.first.back().previous_row;
 
  // 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 = 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 += 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) -
                              Aa(sequence.at(previous_row - 1).unicode()).getMass()) +
              ']');
            m_best_alignment.shifts.push_back(
              spectrum->getMZShift(cell.previous_column, previous_column) -
              Aa(sequence.at(previous_row - 1).unicode()).getMass());
            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());
 
  // 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;
    }
 
  // 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;
        }
    }
 
  // 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);
        }
    }
 
  // 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) + ']';
    }
}
 
const pappso::specpeptidoms::Alignment &
pappso::specpeptidoms::SemiGlobalAlignment::getBestAlignment(const SpOMSSpectrumCsp &spectrum) const
{
  return m_best_alignment;
}
 
std::vector<double>
pappso::specpeptidoms::SemiGlobalAlignment::getPotentialMassErrors(const Alignment &alignment,
                                                                   const QString &protein_seq)
{
  std::vector<double> potential_mass_errors(alignment.shifts);
  double shift      = alignment.end_shift;
  std::size_t index = alignment.beginning - 1;
  while(shift > 0 && index > 0)
    {
      potential_mass_errors.push_back(shift);
      shift -= Aa(protein_seq.at(index).unicode()).getMass();
      index--;
    }
  shift = alignment.begin_shift;
  index = alignment.end + 1;
  while(shift > 0 && index < (std::size_t)protein_seq.size())
    {
      potential_mass_errors.push_back(shift);
      shift -= Aa(protein_seq.at(index).unicode()).getMass();
      index++;
    }
  return potential_mass_errors;
}