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 <algorithm>
#include "semiglobalalignment.h"
#include "pappsomspp/core/processing/specpeptidoms/types.h"
#include "pappsomspp/core/amino_acid/aacode.h"
#include "pappsomspp/core/protein/protein.h"
#include "correctiontree.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
 
void
pappso::specpeptidoms::Alignment::reset()
{
  peaks.clear();
  m_peptideModel.reset();
  score       = 0;
  begin_shift = 0.0;
  end_shift   = 0.0;
  shifts.clear();
  SPC       = 0;
  beginning = 0;
  end       = 0;
}
 
 
QString
pappso::specpeptidoms::Alignment::getPeptideString(const QString &protein_sequence) const
{
  return protein_sequence.mid(beginning, end - beginning);
}
 
double
pappso::specpeptidoms::Alignment::getNonAlignedMass() const
{
  double sum_of_elems = std::accumulate(shifts.begin(), shifts.end(), 0);
  return begin_shift + sum_of_elems + end_shift;
}
 
 
std::size_t
pappso::specpeptidoms::Alignment::getPositionStart() const
{
  return beginning;
}
 
 
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;
 
  KeyCell key_cell_init_first;
  key_cell_init_first.n_row     = 0;
  key_cell_init_first.score     = 0;
  key_cell_init_first.beginning = 0;
  key_cell_init_first.tree_id   = 0;
  m_interest_cells.push_back(key_cell_init_first);
}
 
const std::vector<pappso::specpeptidoms::KeyCell> &
pappso::specpeptidoms::SemiGlobalAlignment::getInterestCells() const
{
  return m_interest_cells;
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::fastAlign(const SpOMSSpectrum &spectrum,
                                                      const SpOMSProtein *protein_ptr)
{
  // m_scenario.clear();
  //  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_ptr->size();
 
  KeyCell key_cell_init;
  key_cell_init.n_row     = 0;
  key_cell_init.score     = m_scorevalues.get(ScoreType::init);
  key_cell_init.beginning = 0;
  key_cell_init.tree_id   = 0;
 
  m_interest_cells.resize(spectrum.size());
  std::fill(m_interest_cells.begin(), m_interest_cells.end(), key_cell_init);
 
  m_interest_cells.at(0).score = 0;
 
  // m_location_saver.resetLocationSaver();
  m_updated_cells.clear();
  for(std::size_t row_number = 1; row_number <= sequence_length; row_number++)
    {
      updateAlignmentMatrix(*protein_ptr,
                            row_number,
                            spectrum.getAaPositions(protein_ptr->at(row_number - 1).code),
                            spectrum,
                            true,
                            protein_ptr);
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::preciseAlign(const SpOMSSpectrum &spectrum,
                                                         const SpOMSProtein *protein_ptr,
                                                         const std::size_t beginning,
                                                         const std::size_t length)
{
  try
    {
      qDebug();
      const QString &protein_seq = protein_ptr->getSequence();
      std::size_t length2;
      if((qsizetype)(beginning + length) <= protein_seq.size())
        {
          length2 = length;
        }
      else
        {
          length2 = protein_seq.size() - beginning;
        }
 
      qDebug();
      QString sequence_str = protein_seq.sliced(protein_seq.size() - beginning - length2, length2);
 
      SpOMSProtein sequence("sub_sequence", sequence_str, m_aaCode);
 
      // 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();
          // aa           = Aa(sequence[row_number - 1].unicode());
          updateAlignmentMatrix(sequence,
                                row_number,
                                spectrum.getAaPositions(sequence[row_number - 1].code),
                                spectrum,
                                false,
                                protein_ptr);
        }
      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
        {
          // we only correct alignment if the sequence has a minimum amino acid diversity
          if(checkSequenceDiversity(sequence.getSequence(), 5, 2))
            {
 
              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(sequence,
                                   protein_ptr,
                                   spectrum,
                                   peaks_to_remove,
                                   protein_seq.size() - beginning);
                      qDebug();
                    }
                  qDebug();
                  m_best_alignment = m_best_corrected_alignment;
                }
            }
        }
      else
        {
          // this sequence has too much redundancy
          // we have to lower the score
          m_best_alignment.score = 0;
        }
      qDebug();
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("SemiGlobalAlignment::preciseAlign failed :\n%1").arg(err.qwhat()));
    }
}
 
void
pappso::specpeptidoms::SemiGlobalAlignment::correctAlign(const SpOMSProtein &sequence,
                                                         const SpOMSProtein *protein_ptr,
                                                         const SpOMSSpectrum &spectrum,
                                                         std::vector<std::size_t> &peaks_to_remove,
                                                         std::size_t offset)
{
  std::vector<AaPosition> aa_positions;
  CorrectionTree correction_tree;
  std::vector<std::size_t> final_peaks_to_remove;
 
  KeyCell key_cell_init;
  key_cell_init.beginning = 0;
  key_cell_init.n_row     = 0;
  key_cell_init.score     = m_scorevalues.get(ScoreType::init);
  key_cell_init.tree_id   = 0;
 
  std::fill(m_interest_cells.begin(), m_interest_cells.end(), key_cell_init);
 
  m_interest_cells.at(0).score = 0;
 
  m_scenario.resetScenario();
  qDebug();
  for(qsizetype row_number = 1; row_number <= sequence.size(); row_number++)
    {
      qDebug() << row_number - 1 << " " << sequence.size();
      qDebug() << "sequence[row_number - 1].aa" << (char)sequence[row_number - 1].aa;
      qDebug();
      aa_positions = spectrum.getAaPositions(sequence[row_number - 1].code, peaks_to_remove);
      qDebug();
      updateAlignmentMatrix(sequence, row_number, aa_positions, spectrum, false, protein_ptr);
      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.getSequence();
      qDebug() << offset;
      qDebug() << spectrum.getPrecursorCharge();
      saveBestAlignment(sequence, spectrum, offset);
      qDebug();
      for(std::size_t iter : m_best_alignment.peaks)
        {
          qDebug() << "iter:" << iter << "comp:" << spectrum.getComplementaryPeak(iter);
          if(iter == spectrum.getComplementaryPeak(iter))
            {
              continue;
            }
          else 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 new_peaks_to_remove : corrections)
            {
              final_peaks_to_remove = std::vector<std::size_t>(new_peaks_to_remove);
              final_peaks_to_remove.insert(
                final_peaks_to_remove.end(), peaks_to_remove.begin(), peaks_to_remove.end());
              correctAlign(sequence, protein_ptr, spectrum, final_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 SpOMSProtein *protein_ptr,
                                                                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_ptr, 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 pappso::specpeptidoms::SpOMSProtein &sequence,
  const std::size_t row_number,
  const std::vector<AaPosition> &aa_positions,
  const SpOMSSpectrum &spectrum,
  const bool fast_align,
  const pappso::specpeptidoms::SpOMSProtein *protein_ptr)
{
  int where = 0;
  try
    {
      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;
 
      double smallest_aa_mass = m_aaCode.getMass((std::uint8_t)1);
 
      m_updated_cells.reserve(aa_positions.size());
      where = 1;
      // Computation of the threePeaks condition, see spomsspectrum.h for more details.
      if(fast_align)
        {
          condition = 3;
          if(row_number > 1)
            {
              qDebug() << (char)sequence.at(row_number - 2).aa;
              qDebug() << "condition" << condition;
              condition += 2 << sequence.at(row_number - 2).code;
              qDebug();
              qDebug() << "condition" << condition;
            }
        }
      where = 2;
      for(std::vector<AaPosition>::const_iterator aa_position = aa_positions.begin();
          aa_position != aa_positions.end();
          aa_position++)
        {
 
          where = 3;
          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;
                        }
                    }
                }
 
              where = 4;
              // 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) &&
                     shift <= TOL_PEAKS_MISSING)
                    {
                      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++;
                }
 
              where = 5;
              // 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
              if(aa_position->r_peak <= TOL_PEAKS_MISSING_FIRST_COLUMN)
                {
                  perfect_shift_origin =
                    perfectShiftPossibleFrom0(sequence, spectrum, row_number, aa_position->r_peak);
                }
              else
                {
                  perfect_shift_origin = row_number;
                }
 
              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;
                }
 
              where = 6;
              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) / smallest_aa_mass);
                  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);
                    }
                  where = 7;
                  if(fast_align)
                    {
                      tree_id = m_location_saver.getNextTree();
                    }
                }
 
              where = 8;
              if(best_column != aa_position->r_peak)
                {
                  m_updated_cells.push_back(
                    {aa_position->r_peak, {row_number, best_score, beginning, tree_id}});
                }
 
              where = 9;
              if(best_score > m_location_saver.getMinScore(tree_id) && fast_align)
                {
                  length =
                    row_number - beginning + 1 +
                    std::ceil(spectrum.getMissingMass(aa_position->r_peak) / smallest_aa_mass) +
                    ALIGNMENT_SURPLUS;
                  where = 10;
                  m_location_saver.addLocation(beginning, length, tree_id, best_score, protein_ptr);
                }
              else if(!fast_align)
                {
 
                  where = 11;
                  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);
                    }
                }
            }
        }
 
      where = 30;
      //  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)
        {
          qDebug() << m_interest_cells.size() << " " << m_updated_cells.back().first;
          m_interest_cells.at(m_updated_cells.back().first) = m_updated_cells.back().second;
          m_updated_cells.pop_back();
        }
      where++;
    }
  catch(const std::exception &error)
    {
      throw pappso::PappsoException(
        QObject::tr("updateAlignmentMatrix failed std::exception :\n%1 %2")
          .arg(where)
          .arg(error.what()));
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("updateAlignmentMatrix failed :\n%1").arg(err.qwhat()));
    }
}
 
bool
pappso::specpeptidoms::SemiGlobalAlignment::perfectShiftPossible(
  const pappso::specpeptidoms::SpOMSProtein &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
{
  try
    {
      double missing_mass = 0;
      auto it_end         = sequence.begin() + current_row;
      for(auto iter = sequence.begin() + origin_row; (iter != it_end) && (iter != sequence.end());
          iter++)
        {
          missing_mass += iter->mass; // Aa(iter->unicode()).getMass();
        }
      if(missing_mass > 0)
        {
          pappso::MzRange mz_range(missing_mass, m_precision_ptr);
          return mz_range.contains(spectrum.getMZShift(l_peak, r_peak));
        }
      else
        {
          return false;
        }
    }
  catch(const std::exception &error)
    {
      throw pappso::PappsoException(
        QObject::tr("perfectShiftPossible failed std exception:\n%1").arg(error.what()));
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("perfectShiftPossible failed :\n%1").arg(err.qwhat()));
    }
}
 
std::size_t
pappso::specpeptidoms::SemiGlobalAlignment::perfectShiftPossibleFrom0(
  const pappso::specpeptidoms::SpOMSProtein &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 += sequence.at(perfect_shift_origin - 1)
                   .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 pappso::specpeptidoms::SpOMSProtein &sequence,
  const SpOMSSpectrum &spectrum,
  std::size_t end_row,
  std::size_t end_peak) const
{
  try
    {
      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 += sequence.at(perfect_shift_end - 1)
                       .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;
        }
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("perfectShiftPossibleEnd failed :\n%1").arg(err.qwhat()));
    }
}
 
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 pappso::specpeptidoms::SpOMSProtein &sequence,
  const SpOMSSpectrum &spectrum,
  std::size_t offset)
{
  qDebug();
  m_best_alignment.m_peptideModel.reset();
  m_best_alignment.peaks.clear();
  m_best_alignment.shifts.clear();
  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;
  std::vector<SpOMSAa> 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 - 1;
 
  qDebug();
  AminoAcidModel aa_model;
  aa_model.m_massDifference = 0;
  // Filling temp_interpretation and peaks vectors
  for(auto cell : best_alignment.first)
    {
      switch(cell.alignment_type)
        {
          case AlignType::found:
            aa_model.m_aminoAcid      = sequence.at(previous_row - 1).aa;
            aa_model.m_massDifference = 0;
            aa_model.m_skipped        = false;
            m_best_alignment.m_peptideModel.push_back(aa_model);
            if(previous_row > cell.previous_row + 1)
              {
                skipped_mass = sequence.at(previous_row - 1)
                                 .mass; // Aa(sequence.at(previous_row - 1).unicode()).getMass();
                skipped_aa =
                  sequence.sliced(cell.previous_row, previous_row - cell.previous_row - 1);
                aa_model.m_massDifference = 0;
                aa_model.m_skipped        = true;
                for(auto aa : skipped_aa)
                  {
                    aa_model.m_aminoAcid = aa.aa;
                    m_best_alignment.m_peptideModel.push_back(aa_model);
                    skipped_mass += aa.mass; // Aa(aa.unicode()).getMass();
                  }
                m_best_alignment.m_peptideModel.back().m_massDifference =
                  spectrum.getMZShift(cell.previous_column, previous_column) - skipped_mass;
              }
            m_best_alignment.peaks.push_back(cell.previous_column);
            break;
          case AlignType::notFound:
            aa_model.m_aminoAcid      = sequence.at(previous_row - 1).aa;
            aa_model.m_massDifference = 0;
            aa_model.m_skipped        = true;
            m_best_alignment.m_peptideModel.push_back(aa_model);
            break;
          case AlignType::shift:
 
            aa_model.m_aminoAcid      = sequence.at(previous_row - 1).aa;
            aa_model.m_massDifference = spectrum.getMZShift(cell.previous_column, previous_column) -
                                        aa_model.m_aminoAcid.getMass();
            aa_model.m_skipped = false;
            m_best_alignment.m_peptideModel.push_back(aa_model);
            m_best_alignment.peaks.push_back(cell.previous_column);
            m_best_alignment.shifts.push_back(
              spectrum.getMZShift(cell.previous_column, previous_column) -
              sequence.at(previous_row - 1).mass);
            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());
            aa_model.m_massDifference = 0;
            aa_model.m_skipped        = false;
            for(auto aa : skipped_aa)
              {
                aa_model.m_aminoAcid = aa.aa;
                m_best_alignment.m_peptideModel.push_back(aa_model);
              }
            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(m_best_alignment.peaks.begin(), m_best_alignment.peaks.end());
  std::reverse(m_best_alignment.m_peptideModel.begin(), m_best_alignment.m_peptideModel.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);
          aa_model.m_massDifference = 0;
          aa_model.m_skipped        = true;
          for(auto aa = skipped_aa.begin(); aa != skipped_aa.end(); aa++)
            {
              aa_model.m_aminoAcid = aa->aa;
              m_best_alignment.m_peptideModel.push_back(aa_model);
            }
          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
  if(m_best_alignment.end_shift > 0)
    {
      m_best_alignment.m_peptideModel.setNterShift(m_best_alignment.end_shift);
    }
 
  std::reverse(m_best_alignment.m_peptideModel.begin(), m_best_alignment.m_peptideModel.end());
  if(m_best_alignment.begin_shift > 0)
    {
      m_best_alignment.m_peptideModel.setCterShift(m_best_alignment.begin_shift);
    }
 
  m_best_alignment.m_peptideModel.setPrecursorMass(spectrum.getPrecursorMass());
  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)
{
  // qDebug() << 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);
          // qDebug() << " 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--;
        }
    }
 
  // qDebug() << "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);
      qDebug() << " 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++;
    }
  // qDebug();
  return potential_mass_errors;
}
 
bool
pappso::specpeptidoms::SemiGlobalAlignment::checkSequenceDiversity(const QString &sequence,
                                                                   std::size_t window,
                                                                   std::size_t minimum_aa_diversity)
{
  qDebug() << "sequence=" << sequence << " window=" << window
           << " minimum_aa_diversity=" << minimum_aa_diversity;
  if(sequence.size() < window)
    return false;
  auto it_begin = sequence.begin();
  auto it_end   = sequence.begin() + window;
  QString window_copy(sequence.mid(0, window));
  while(it_end != sequence.end())
    {
      std::partial_sort_copy(it_begin, it_end, window_copy.begin(), window_copy.end());
 
      qDebug() << window_copy;
      std::size_t uniqueCount =
        std::unique(window_copy.begin(), window_copy.end()) - window_copy.begin();
 
      qDebug() << uniqueCount;
      if(uniqueCount < minimum_aa_diversity)
        return false;
      it_begin++;
      it_end++;
    }
  return true;
}