aastringcodemassmatching.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/amino_acid/aastringcodemassmatching.cpp
 * \date 10/05/2023
 * \author Olivier Langella
 * \brief convert mass list to amino acid string code list
 */
 
/*******************************************************************************
 * Copyright (c) 2023 Olivier Langella <Olivier.Langella@u-psud.fr>.
 *
 * This file is part of PAPPSOms-tools.
 *
 *     PAPPSOms-tools 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.
 *
 *     PAPPSOms-tools 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 PAPPSOms-tools.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "pappsomspp/core/amino_acid/aastringcodemassmatching.h"
#include "pappsomspp/core/amino_acid/aastringcodec.h"
#include <QDebug>
#include <algorithm>
 
using namespace pappso;
 
AaStringCodeMassMatching::AaStringCodeMassMatching(const AaCode &aa_code,
                                                   std::size_t model_max_size,
                                                   PrecisionPtr precision)
  : m_precision(precision), m_aaCode(aa_code), m_aaCodec(aa_code)
{
  m_base = m_aaCode.getSize() + 1;
 
  std::vector<pappso::CodeToMass> code_to_mass =
    m_aaCodec.generateLlcCodeListUpToMaxPeptideSize(model_max_size);
 
  for(auto &code_mass : code_to_mass)
    {
      aaCodeAndMassRange aaCodeMassRange;
      aaCodeMassRange.code         = code_mass.code;
      double delta                 = precision->delta(code_mass.mass);
      aaCodeMassRange.mz_range_low = code_mass.mass - delta;
      aaCodeMassRange.mz           = code_mass.mass;
      aaCodeMassRange.mz_range_up  = code_mass.mass + delta;
 
      m_codeMassList.push_back(aaCodeMassRange);
    }
 
 
  std::sort(m_codeMassList.begin(),
            m_codeMassList.end(),
            [](const aaCodeAndMassRange &a, const aaCodeAndMassRange &b) {
              return a.mz_range_low < b.mz_range_low;
            });
}
 
AaStringCodeMassMatching::AaStringCodeMassMatching(
  const AaStringCodeMassMatching &other)
  : m_precision(other.m_precision),
    m_aaCode(other.m_aaCode),
    m_aaCodec(other.m_aaCode)
{
  m_base = m_aaCode.getSize() + 1;
}
 
AaStringCodeMassMatching::~AaStringCodeMassMatching()
{
}
 
 
pappso::PrecisionPtr
pappso::AaStringCodeMassMatching::getPrecisionPtr() const
{
  return m_precision;
}
 
 
std::vector<uint32_t>
pappso::AaStringCodeMassMatching::getAaCodeFromMass(double mass) const
{
  std::vector<uint32_t> aa_code_list;
 
  // auto it_aacode = m_codeMassList.begin();
 
  auto it_aacode = std::upper_bound(
    m_codeMassList.begin(),
    m_codeMassList.end(),
    mass,
    [](double mass,
       const pappso::AaStringCodeMassMatching::aaCodeAndMassRange &mass_range) {
      return mass_range.mz_range_up > mass;
    });
 
  bool previous_out_of_range = false;
 
  while(it_aacode != m_codeMassList.end())
    {
      if(mass < it_aacode->mz_range_low)
        {
 
          if(previous_out_of_range)
            break;
          previous_out_of_range = true;
          it_aacode++;
        }
      else
        {
          if(mass <= it_aacode->mz_range_up)
            {
              previous_out_of_range = false;
              aa_code_list.push_back(it_aacode->code);
              it_aacode++;
            }
          else
            {
              it_aacode++;
            }
        }
    }
  return aa_code_list;
}
 
std::vector<uint32_t>
pappso::AaStringCodeMassMatching::getAaCodeFromMassWearingModification(
  double mass, const pappso::Aa &aa, int quantifier) const
{
  std::vector<uint32_t> aa_code_list;
 
  double total_modification_mass = aa.getTotalModificationMass();
  total_modification_mass *= quantifier;
  mass = mass - total_modification_mass;
  if(mass < 0)
    return aa_code_list;
 
  uint8_t aamodCode = m_aaCode.getAaCode(aa.getLetter());
 
  auto it_aacode = m_codeMassList.begin();
 
  while(it_aacode != m_codeMassList.end())
    {
      if(mass < it_aacode->mz_range_low)
        {
          it_aacode++;
        }
      else
        {
          if(mass <= it_aacode->mz_range_up)
            {
              if(m_aaCodec.uniqueCodeContainsAminoAcid(
                   it_aacode->code, aamodCode, quantifier))
                {
                  aa_code_list.push_back(it_aacode->code);
                }
 
              it_aacode++;
            }
          else
            {
              it_aacode++;
            }
        }
    }
  return aa_code_list;
}
 
 
std::vector<uint32_t>
pappso::AaStringCodeMassMatching::getAaCodeFromMassList(
  std::vector<double> &mass_list) const
{
  std::sort(mass_list.begin(), mass_list.end(), [](double a, double b) {
    return a < b;
  });
  std::vector<uint32_t> aa_code_list;
 
  auto it_aacode = m_codeMassList.begin();
  auto it_mass   = mass_list.begin();
 
  while((it_aacode != m_codeMassList.end()) && (it_mass != mass_list.end()))
    {
      if(*it_mass < it_aacode->mz_range_low)
        {
          it_mass++;
        }
      else
        {
          if(*it_mass <= it_aacode->mz_range_up)
            {
              aa_code_list.push_back(it_aacode->code);
              it_aacode++;
            }
          else
            {
              it_aacode++;
            }
        }
    }
  return aa_code_list;
}
 
std::vector<uint32_t>
pappso::AaStringCodeMassMatching::filterCodeList(
  std::vector<uint32_t> &code_list) const
{
  std::sort(code_list.begin(), code_list.end(), [](uint32_t a, uint32_t b) {
    return a < b;
  });
  std::vector<uint32_t> filtered_aa_code_list;
 
  std::vector<uint8_t> aa_ok;
 
  auto it = code_list.begin();
  while(*it < m_base)
    {
      aa_ok.push_back((uint8_t)*it);
      // qDebug() << (uint8_t)*it << " "
      //       << m_aaCode.getAa((uint8_t)*it).getLetter();
      it++;
    }
 
  for(uint32_t code : code_list)
    {
      if(m_aaCodec.codeOnlyContains(code, aa_ok))
        {
          filtered_aa_code_list.push_back(code);
        }
    }
  return filtered_aa_code_list;
}