#include <semiglobalalignment.h>

Public Member Functions
	SemiGlobalAlignment (const ScoreValues &score_values, const pappso::PrecisionPtr precision_ptr, const AaCode &aaCode)

	~SemiGlobalAlignment ()

void	fastAlign (const SpOMSSpectrum &spectrum, const SpOMSProtein *protein_ptr)
	perform the first alignment search between a protein sequence and a spectrum. The member location heap is filled with the candidates locations.

void	preciseAlign (const SpOMSSpectrum &spectrum, const SpOMSProtein *protein_ptr, const std::size_t beginning, const std::size_t length)
	performs the second alignment search between a protein subsequence and a spectrum.

void	postProcessingAlign (const SpOMSSpectrum &spectrum, const SpOMSProtein *protein_ptr, std::size_t beginning, std::size_t length, const std::vector< double > &shifts)
	performs the post-processing : generates corrected spectra and align them

LocationSaver	getLocationSaver () const
	Returns a copy of m_location_saver.

Scenario	getScenario () const
	Returns a copy of m_scenario.

const Alignment &	getBestAlignment () const
	Returns a const ref to m_best_alignment.

const std::vector< KeyCell > &	getInterestCells () const
	convenient function for degub purpose

Static Public Member Functions
static std::vector< double >	getPotentialMassErrors (const pappso::AaCode &aa_code, const Alignment &alignment, const QString &protein_seq)
	Returns a list of the potential mass errors corresponding to the provided alignment in the provided protein sequence.

static bool	checkSequenceDiversity (const QString &sequence, std::size_t window, std::size_t minimum_aa_diversity)
	check that the sequence has a minimum of amino acid checkSequenceDiversity

Private Member Functions
void	saveBestAlignment (const SpOMSProtein &sequence, const SpOMSSpectrum &spectrum, std::size_t offset)
	Stores the best alignment from m_scenario in m_best_alignment.

void	correctAlign (const SpOMSProtein &protein_subseq, const SpOMSProtein *protein_ptr, const SpOMSSpectrum &spectrum, std::vector< std::size_t > &peaks_to_remove, std::size_t offset)
	Recursively performs the correction of the alignment.

void	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)
	updates the scores of the alignment matrix for a given amino acid as well as the location heap/scenario.

bool	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
	indicates if a perfect shift is possible between the provided positions

std::size_t	perfectShiftPossibleFrom0 (const pappso::specpeptidoms::SpOMSProtein &sequence, const SpOMSSpectrum &spectrum, const std::size_t current_row, const std::size_t r_peak) const
	indicates if a perfect shift is possible from the spectrum beginning to the provided peak.

std::size_t	perfectShiftPossibleEnd (const pappso::specpeptidoms::SpOMSProtein &sequence, const SpOMSSpectrum &spectrum, std::size_t end_row, std::size_t end_peak) const
	indicates if a perfect shift is possible between the provided positions

Private Attributes
std::vector< KeyCell >	m_interest_cells

std::vector< std::pair< std::size_t, KeyCell > >	m_updated_cells

const ScoreValues &	m_scorevalues

const int	min_score = 15

pappso::PrecisionPtr	m_precision_ptr

const AaCode &	m_aaCode

LocationSaver	m_location_saver

Scenario	m_scenario

Alignment	m_best_alignment

Alignment	m_best_corrected_alignment

Alignment	m_best_post_processed_alignment

Detailed Description

Definition at line 90 of file semiglobalalignment.h.

Constructor & Destructor Documentation

◆ SemiGlobalAlignment()

pappso::specpeptidoms::SemiGlobalAlignment::SemiGlobalAlignment	(	const ScoreValues &	score_values,
		const pappso::PrecisionPtr	precision_ptr,
		const AaCode &	aaCode
	)

Default constructor

Definition at line 80 of file semiglobalalignment.cpp.

  : 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);
}

References pappso::specpeptidoms::KeyCell::beginning, m_interest_cells, m_precision_ptr, pappso::specpeptidoms::KeyCell::n_row, pappso::specpeptidoms::KeyCell::score, and pappso::specpeptidoms::KeyCell::tree_id.

◆ ~SemiGlobalAlignment()

pappso::specpeptidoms::SemiGlobalAlignment::~SemiGlobalAlignment ( )

Destructor

Definition at line 752 of file semiglobalalignment.cpp.

753{

754}

Member Function Documentation

◆ checkSequenceDiversity()

bool pappso::specpeptidoms::SemiGlobalAlignment::checkSequenceDiversity	(	const QString &	sequence,
		std::size_t	window,
		std::size_t	minimum_aa_diversity
	)

static

check that the sequence has a minimum of amino acid checkSequenceDiversity

Parameters

sequence	protein sequence
window	the size of substring to check
minimum_aa_diversity	minimum number of different amino acid in this window

Definition at line 1010 of file semiglobalalignment.cpp.

{
  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;
}

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment(), and pappso::cbor::psm::PsmSpecPeptidOmsScan::storeAlignment().

◆ correctAlign()

void pappso::specpeptidoms::SemiGlobalAlignment::correctAlign	(	const SpOMSProtein &	protein_subseq,
		const SpOMSProtein *	protein_ptr,
		const SpOMSSpectrum &	spectrum,
		std::vector< std::size_t > &	peaks_to_remove,
		std::size_t	offset
	)

private

Recursively performs the correction of the alignment.

Parameters

protein_seq	Protein reversed sequence to align.
protein_ptr	Protein pointer on the sequence to align.
spectrum	Spectrum to align.
peaks_to_remove	Peaks to remove from the spectrum.
offset	Size of the protein sequence minus beginning of the alignment. Used to compute the position of the alignment in the protein sequence.

Definition at line 269 of file semiglobalalignment.cpp.

{
  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();
}

◆ fastAlign()

void pappso::specpeptidoms::SemiGlobalAlignment::fastAlign	(	const SpOMSSpectrum &	spectrum,
		const SpOMSProtein *	protein_ptr
	)

perform the first alignment search between a protein sequence and a spectrum. The member location heap is filled with the candidates locations.

Parameters

spectrum	Spectrum to align
protein_ptr	Protein pointer on the sequence to align.

Definition at line 103 of file semiglobalalignment.cpp.

{
  // 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);
    }
}

References pappso::specpeptidoms::KeyCell::beginning, pappso::specpeptidoms::SpOMSSpectrum::getAaPositions(), pappso::specpeptidoms::init, pappso::specpeptidoms::KeyCell::n_row, pappso::specpeptidoms::KeyCell::score, and pappso::specpeptidoms::KeyCell::tree_id.

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ getBestAlignment()

const pappso::specpeptidoms::Alignment & pappso::specpeptidoms::SemiGlobalAlignment::getBestAlignment ( ) const

Returns a const ref to m_best_alignment.

Definition at line 959 of file semiglobalalignment.cpp.

{
  return m_best_alignment;
}

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ getInterestCells()

const std::vector< pappso::specpeptidoms::KeyCell > & pappso::specpeptidoms::SemiGlobalAlignment::getInterestCells ( ) const

convenient function for degub purpose

Definition at line 97 of file semiglobalalignment.cpp.

{
  return m_interest_cells;
}

◆ getLocationSaver()

pappso::specpeptidoms::LocationSaver pappso::specpeptidoms::SemiGlobalAlignment::getLocationSaver ( ) const

Returns a copy of m_location_saver.

Definition at line 757 of file semiglobalalignment.cpp.

{
  return m_location_saver;
}

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ getPotentialMassErrors()

std::vector< double > pappso::specpeptidoms::SemiGlobalAlignment::getPotentialMassErrors	(	const pappso::AaCode &	aa_code,
		const Alignment &	alignment,
		const QString &	protein_seq
	)

static

Returns a list of the potential mass errors corresponding to the provided alignment in the provided protein sequence.

Parameters

aa_code	the amino acid code of reference to get aminon acid masses
alignment	Alignment for which to get the potential mass errors.
protein_seq	Protein sequence corresponding to the provided alignment.

Definition at line 965 of file semiglobalalignment.cpp.

{
  // 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;
}

References pappso::specpeptidoms::Alignment::begin_shift, pappso::specpeptidoms::Alignment::beginning, pappso::specpeptidoms::Alignment::end, pappso::specpeptidoms::Alignment::end_shift, pappso::AaCode::getMass(), pappso::specpeptidoms::shift, and pappso::specpeptidoms::Alignment::shifts.

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ getScenario()

pappso::specpeptidoms::Scenario pappso::specpeptidoms::SemiGlobalAlignment::getScenario ( ) const

Returns a copy of m_scenario.

Definition at line 763 of file semiglobalalignment.cpp.

{
  return m_scenario;
}

◆ perfectShiftPossible()

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

private

indicates if a perfect shift is possible between the provided positions

Parameters

sequence	Reversed sequence of the protein being aligned
spectrum	Spectrum being aligned
origin_row	beginning row of the aa gap to verify (== index of the first missing aa in sequence)
current_row	row being processed (== index of the current AaPosition in sequence)
l_peak	left peak index of the mz gap to verify
r_peak	right peak index of the mz gap to verify

Definition at line 650 of file semiglobalalignment.cpp.

{
  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 &stderr)
    {
      throw pappso::PappsoException(
        QObject::tr("perfectShiftPossible failed std exception:\n%1").arg(stderr.what()));
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("perfectShiftPossible failed :\n%1").arg(err.qwhat()));
    }
}

References pappso::MzRange::contains(), pappso::specpeptidoms::SpOMSSpectrum::getMZShift(), and pappso::PappsoException::qwhat().

◆ perfectShiftPossibleEnd()

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

private

indicates if a perfect shift is possible between the provided positions

Parameters

sequence	Reversed sequence of the protein being aligned
spectrum	Spectrum being aligned
end_row	Index of the last aligned row.
end_peak	Index of the last aligned peak.

Definition at line 717 of file semiglobalalignment.cpp.

{
  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()));
    }
}

References pappso::MzRange::contains(), pappso::specpeptidoms::SpOMSSpectrum::getMissingMass(), and pappso::PappsoException::qwhat().

◆ perfectShiftPossibleFrom0()

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

private

indicates if a perfect shift is possible from the spectrum beginning to the provided peak.

Parameters

sequence	Reversed sequence of the protein being aligned
spectrum	Spectrum being aligned
current_row	row being processed (== index of the current AaPosition in sequence)
r_peak	right peak index of the mz gap to verify

Definition at line 690 of file semiglobalalignment.cpp.

{
  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;
    }
}

References pappso::MzRange::contains(), and pappso::specpeptidoms::SpOMSSpectrum::getMZShift().

◆ postProcessingAlign()

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
	)

performs the post-processing : generates corrected spectra and align them

Parameters

spectrum	Spectrum to align
protein_ptr	Protein pointer on the sequence to align.
beginning	Index of the beginning of the subsequence to align.
length	Length of the subsequence to align.
shifts	List of potential precursor mass errors to test.

Definition at line 353 of file semiglobalalignment.cpp.

{
  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;
    }
}

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ preciseAlign()

void pappso::specpeptidoms::SemiGlobalAlignment::preciseAlign	(	const SpOMSSpectrum &	spectrum,
		const SpOMSProtein *	protein_ptr,
		const std::size_t	beginning,
		const std::size_t	length
	)

performs the second alignment search between a protein subsequence and a spectrum.

Parameters

spectrum	Spectrum to align
protein_ptr	Protein pointer on the sequence to align.
beginning	Index of the beginning of the subsequence to align.
length	Length of the subsequence to align.

Definition at line 139 of file semiglobalalignment.cpp.

{
  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()));
    }
}

References pappso::specpeptidoms::CorrectionTree::addPeaks(), pappso::specpeptidoms::SpOMSSpectrum::getAaPositions(), pappso::specpeptidoms::SpOMSSpectrum::getComplementaryPeak(), pappso::specpeptidoms::CorrectionTree::getPeaks(), pappso::specpeptidoms::SpOMSProtein::getSequence(), pappso::specpeptidoms::init, pappso::specpeptidoms::MIN_ALIGNMENT_SCORE(), and pappso::PappsoException::qwhat().

Referenced by pappso::cbor::psm::PsmSpecPeptidOmsScan::sequenceAlignment().

◆ saveBestAlignment()

void pappso::specpeptidoms::SemiGlobalAlignment::saveBestAlignment	(	const SpOMSProtein &	sequence,
		const SpOMSSpectrum &	spectrum,
		std::size_t	offset
	)

private

Stores the best alignment from m_scenario in m_best_alignment.

Parameters

sequence	reversed sequence of the current alignment.
spectrum	Spectrum currently being aligned.
offset	Size of the protein sequence minus beginning of the alignment. Used to compute the position of the alignment in the protein sequence.

Definition at line 769 of file semiglobalalignment.cpp.

{
  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;
 
  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();
}

◆ updateAlignmentMatrix()

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
	)

private

updates the scores of the alignment matrix for a given amino acid as well as the location heap/scenario.

Parameters

sequence	Reversed sequence of the protein being aligned
row_number	number of the row to update (== index in sequence of the amino acid being aligned)
aa_positions	list of the AaPositions of the current amino acid
spectrum	Spectrum being aligned
fast_align	Whether to use the fast version of the algorithm (for 1st alignemnt step)
protein_ptr	Protein pointer on the sequence to align.

Definition at line 379 of file semiglobalalignment.cpp.

{
  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))
                    {
                      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
              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;
                }
 
              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 &stderr)
    {
      throw pappso::PappsoException(
        QObject::tr("updateAlignmentMatrix failed std::exception :\n%1 %2")
          .arg(where)
          .arg(stderr.what()));
    }
  catch(const pappso::PappsoException &err)
    {
      throw pappso::PappsoException(
        QObject::tr("updateAlignmentMatrix failed :\n%1").arg(err.qwhat()));
    }
}

Member Data Documentation

◆ m_aaCode

const AaCode& pappso::specpeptidoms::SemiGlobalAlignment::m_aaCode

private

Definition at line 275 of file semiglobalalignment.h.

◆ m_best_alignment

Alignment pappso::specpeptidoms::SemiGlobalAlignment::m_best_alignment

private

Definition at line 278 of file semiglobalalignment.h.

◆ m_best_corrected_alignment

Alignment pappso::specpeptidoms::SemiGlobalAlignment::m_best_corrected_alignment

private

Definition at line 279 of file semiglobalalignment.h.

◆ m_best_post_processed_alignment

Alignment pappso::specpeptidoms::SemiGlobalAlignment::m_best_post_processed_alignment

private

Definition at line 280 of file semiglobalalignment.h.

◆ m_interest_cells

std::vector<KeyCell> pappso::specpeptidoms::SemiGlobalAlignment::m_interest_cells

private

Definition at line 270 of file semiglobalalignment.h.

Referenced by SemiGlobalAlignment().

◆ m_location_saver

LocationSaver pappso::specpeptidoms::SemiGlobalAlignment::m_location_saver

private

Definition at line 276 of file semiglobalalignment.h.

◆ m_precision_ptr

pappso::PrecisionPtr pappso::specpeptidoms::SemiGlobalAlignment::m_precision_ptr

private

Definition at line 274 of file semiglobalalignment.h.

Referenced by SemiGlobalAlignment().

◆ m_scenario

Scenario pappso::specpeptidoms::SemiGlobalAlignment::m_scenario

private

Definition at line 277 of file semiglobalalignment.h.

◆ m_scorevalues

const ScoreValues& pappso::specpeptidoms::SemiGlobalAlignment::m_scorevalues

private

Definition at line 272 of file semiglobalalignment.h.

◆ m_updated_cells

std::vector<std::pair<std::size_t, KeyCell> > pappso::specpeptidoms::SemiGlobalAlignment::m_updated_cells

private

Definition at line 271 of file semiglobalalignment.h.

◆ min_score

const int pappso::specpeptidoms::SemiGlobalAlignment::min_score = 15

private

Definition at line 273 of file semiglobalalignment.h.

The documentation for this class was generated from the following files:

pappsomspp/core/processing/specpeptidoms/semiglobalalignment.h
pappsomspp/core/processing/specpeptidoms/semiglobalalignment.cpp

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ SemiGlobalAlignment()

◆ ~SemiGlobalAlignment()

Member Function Documentation

◆ checkSequenceDiversity()

◆ correctAlign()

◆ fastAlign()

◆ getBestAlignment()

◆ getInterestCells()

◆ getLocationSaver()

◆ getPotentialMassErrors()

◆ getScenario()

◆ perfectShiftPossible()

◆ perfectShiftPossibleEnd()

◆ perfectShiftPossibleFrom0()

◆ postProcessingAlign()

◆ preciseAlign()

◆ saveBestAlignment()

◆ updateAlignmentMatrix()

Member Data Documentation

◆ m_aaCode

◆ m_best_alignment

◆ m_best_corrected_alignment

◆ m_best_post_processed_alignment

◆ m_interest_cells

◆ m_location_saver

◆ m_precision_ptr

◆ m_scenario

◆ m_scorevalues

◆ m_updated_cells

◆ min_score