locationsaver.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/processing/specpeptidoms/locationsaver.cpp
 * \date 24/03/2025
 * \author Aurélien Berthier
 * \brief save protein subsequences for 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 <algorithm>
#include "locationsaver.h"
#include "pappsomspp/core/processing/specpeptidoms/types.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
 
namespace pappso
{
namespace specpeptidoms
{
 
QString
Location::getPeptideString() const
{
  std::size_t length2;
  if((qsizetype)(beginning + length) <= proteinPtr->size())
    {
      length2 = length;
    }
  else
    {
      length2 = proteinPtr->size() - beginning;
    }
  return proteinPtr->getSequence().sliced(proteinPtr->size() - beginning - length2, length2);
}
 
 
pappso::specpeptidoms::LocationSaver::LocationSaver()
{
  Location location_zero;
  location_zero.beginning  = 0;
  location_zero.length     = 0;
  location_zero.proteinPtr = nullptr;
  location_zero.score      = MIN_ALIGNMENT_SCORE;
  location_zero.tree       = -1;
 
  m_locations_heap.resize(MAX_SAVED_ALIGNMENTS, location_zero);
  // std::make_heap(m_locations_heap.begin(), m_locations_heap.end(),
  // LocationSaver::locationCompare); ?
  // m_tree_scores.reserve ?
}
 
bool
pappso::specpeptidoms::LocationSaver::locationCompare(const Location &loc1, const Location &loc2)
{
  return loc1.score > loc2.score;
}
 
// TODO : check complexity if m_tree_in_heap[tree] == true
void
pappso::specpeptidoms::LocationSaver::addLocation(
  std::size_t beginning, std::size_t length, int tree, int score, const SpOMSProtein *protein_ptr)
{
  try
    {
      m_tree_scores.at(tree) = score;
      if(m_tree_in_heap.at(tree))
        {
          for(std::vector<Location>::iterator iter = m_locations_heap.begin();
              iter != m_locations_heap.end();
              iter++)
            {
              if(iter->tree == tree)
                {
                  iter->score  = score;
                  iter->length = length;
                }
            }
          std::make_heap(
            m_locations_heap.begin(), m_locations_heap.end(), LocationSaver::locationCompare);
        }
      else
        {
          if(m_locations_heap.begin()->tree >= 0)
            {
              m_tree_in_heap.at(m_locations_heap.begin()->tree) = false;
            }
          m_tree_in_heap.at(tree) = true;
          std::pop_heap(
            m_locations_heap.begin(), m_locations_heap.end(), LocationSaver::locationCompare);
          m_locations_heap.pop_back();
          m_locations_heap.push_back({beginning, length, tree, score, protein_ptr});
          std::push_heap(
            m_locations_heap.begin(), m_locations_heap.end(), LocationSaver::locationCompare);
        }
    }
  catch(const std::exception &error)
    {
      throw pappso::PappsoException(
        QObject::tr("addLocation failed std::exception :\n%1").arg(error.what()));
    }
}
 
std::vector<pappso::specpeptidoms::Location>
pappso::specpeptidoms::LocationSaver::getLocations() const
{
  std::vector<Location> locations;
  locations.reserve(m_locations_heap.size());
  for(std::vector<Location>::const_iterator iter = m_locations_heap.begin();
      iter != m_locations_heap.end();
      iter++)
    {
      if(iter->tree >= 0)
        {
          locations.push_back(*iter);
        }
    }
  return locations;
}
 
std::size_t
pappso::specpeptidoms::LocationSaver::getNextTree()
{
  m_tree_scores.push_back(MIN_ALIGNMENT_SCORE);
  m_tree_in_heap.push_back(false);
  return m_tree_scores.size() - 1;
}
 
int
pappso::specpeptidoms::LocationSaver::getMinScore(int tree_id) const
{
  if(m_tree_scores.size() == 0)
    {
      return m_locations_heap.begin()->score;
    }
  else
    {
      if(tree_id > (int)m_tree_scores.size())
        {
 
          throw pappso::ExceptionOutOfRange(
            QObject::tr("LocationSaver::getMinScore failed :\nout of "
                        "range access %1 with m_tree_scores.size() %2")
              .arg(tree_id)
              .arg(m_tree_scores.size()));
        }
 
      /*To be stored, an alignment's score must be higher than: its tree's current score, the
       * current worse stored alignment's score and a score that depends of the current best
       * alignment's score
       */
      return std::max(
        m_tree_scores.at(tree_id),
        std::max(m_locations_heap.begin()->score,
                 (int)std::ceil(ALIGNMENT_FILTER *
                                *(std::max_element(m_tree_scores.begin(), m_tree_scores.end())))));
    }
}
 
pappso::specpeptidoms::LocationSaver::~LocationSaver()
{
}
 
void
pappso::specpeptidoms::LocationSaver::resetLocationSaver()
{
  Location location_zero;
  location_zero.beginning  = 0;
  location_zero.length     = 0;
  location_zero.proteinPtr = nullptr;
  location_zero.score      = MIN_ALIGNMENT_SCORE;
  location_zero.tree       = -1;
 
  std::fill(m_locations_heap.begin(), m_locations_heap.end(), location_zero);
  m_tree_scores.clear();
  m_tree_in_heap.clear();
 
 
  // int m_min_score, m_max_score;
}
} // namespace specpeptidoms
} // namespace pappso