proteinintegercode.cpp

Go to the documentation of this file.

/**
 * \file protein/proteinintegercode.cpp
 * \date 22/05/2023
 * \author Olivier Langella
 * \brief transform protein amino acid sequence into vectors of amino acid codes
 */
 
 
/*******************************************************************************
 * Copyright (c) 2023 Olivier Langella
 *<Olivier.Langella@universite-paris-saclay.fr>.
 *
 * This file is part of PAPPSOms++.
 *
 *     PAPPSOms++ 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++ 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++.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "proteinintegercode.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
 
using namespace pappso;
 
ProteinIntegerCode::ProteinIntegerCode(ProteinSp protein,
                                       const AaStringCodec &codec,
                                       std::size_t aa_str_max_size)
{
  msp_protein = protein;
 
  if(aa_str_max_size > 7)
    {
      throw ExceptionOutOfRange(
        QObject::tr("aa_str_max_size exceeds max size"));
    }
 
  QString seq_str = protein.get()->getSequence();
  m_seqAaCode.clear();
 
  for(const QChar &aa_str : seq_str)
    {
      m_seqAaCode.push_back(codec.getAaCode().getAaCode(aa_str.toLatin1()));
    }
 
  for(std::size_t i = 2; i <= aa_str_max_size; i++)
    {
      m_peptideCodedFragments.push_back(computePeptideCodeFragments(codec, i));
    }
}
 
ProteinIntegerCode::ProteinIntegerCode(const ProteinIntegerCode &other)
{
 
  msp_protein = other.msp_protein;
 
  m_seqAaCode             = other.m_seqAaCode;
  m_peptideCodedFragments = other.m_peptideCodedFragments;
}
 
ProteinIntegerCode::~ProteinIntegerCode()
{
}
const std::vector<std::uint8_t> &
pappso::ProteinIntegerCode::getSeqAaCode() const
{
  return m_seqAaCode;
}
 
pappso::ProteinSp
pappso::ProteinIntegerCode::getProteinSp() const
{
  return msp_protein;
}
 
 
std::vector<std::uint32_t>
pappso::ProteinIntegerCode::computePeptideCodeFragments(
  const AaStringCodec &codec, std::size_t fragment_size) const
{
  std::vector<std::uint32_t> fragments;
 
  int max = (m_seqAaCode.size() - fragment_size);
  if(max < 0)
    return fragments;
 
  auto it = m_seqAaCode.begin();
  for(int i = 0; i <= max; i++)
    {
      fragments.push_back(codec.codeLlc(it, fragment_size));
      it++;
    }
 
  return fragments;
}
 
 
const std::vector<std::uint32_t> &
pappso::ProteinIntegerCode::getPeptideCodedFragment(std::size_t size) const
{
  if(size < 2)
    {
 
      throw ExceptionOutOfRange(QObject::tr("size too small"));
    }
  std::size_t indice = size - 2;
  if(indice < m_peptideCodedFragments.size())
    {
      return m_peptideCodedFragments.at(indice);
    }
 
  throw ExceptionOutOfRange(QObject::tr("size too big"));
}
 
 
std::vector<std::pair<std::size_t, std::uint32_t>>
pappso::ProteinIntegerCode::match(
  const std::vector<uint32_t> &code_list_in) const
{
  std::vector<std::pair<std::size_t, std::uint32_t>> return_pos;
  std::vector<uint32_t> code_list = code_list_in;
 
  std::sort(code_list.begin(), code_list.end());
  auto it_end = std::unique(code_list.begin(), code_list.end());
  for(auto it_code = code_list.begin(); it_code != it_end; it_code++)
    {
 
      std::size_t size = 2;
      for(auto &liste_protein_seq_code : m_peptideCodedFragments)
        {
 
          auto it_seq_position = std::find(liste_protein_seq_code.begin(),
                                           liste_protein_seq_code.end(),
                                           *it_code);
          while(it_seq_position != liste_protein_seq_code.end())
            {
              // found
              std::size_t position =
                std::distance(liste_protein_seq_code.begin(), it_seq_position);
              return_pos.push_back({size, position});
 
              it_seq_position = std::find(
                ++it_seq_position, liste_protein_seq_code.end(), *it_code);
              qDebug();
            }
          size++;
          qDebug();
        }
      qDebug();
    }
 
  return return_pos;
}
 
std::vector<double>
pappso::ProteinIntegerCode::convolution(
  const std::vector<uint32_t> &code_list_from_spectrum) const
{
  std::vector<double> convolution_score;
 
 
  std::vector<std::uint8_t>::const_iterator it_aa = m_seqAaCode.begin();
  auto it_couple = m_peptideCodedFragments[0].begin();
  auto it_trio   = m_peptideCodedFragments[1].begin();
  auto it_quatro = m_peptideCodedFragments[2].begin();
  auto it_cinqo  = m_peptideCodedFragments[3].begin();
  for(std::uint8_t aa_code [[maybe_unused]] : m_seqAaCode)
    {
      convolution_score.push_back(convolutionKernel(code_list_from_spectrum,
                                                    it_aa,
                                                    it_couple,
                                                    it_trio,
                                                    it_quatro,
                                                    it_cinqo));
      it_aa++;
      it_couple++;
      it_trio++;
      it_quatro++;
      it_cinqo++;
    }
 
  return convolution_score;
}
 
double
pappso::ProteinIntegerCode::convolutionKernel(
  const std::vector<uint32_t> &spectrum_code_list,
  std::vector<std::uint8_t>::const_iterator it_aa,
  std::vector<std::uint32_t>::const_iterator it_couple,
  std::vector<std::uint32_t>::const_iterator it_trio,
  std::vector<std::uint32_t>::const_iterator it_quatro,
  std::vector<std::uint32_t>::const_iterator it_cinqo) const
{
  double score = 0;
 
  // single :
  auto it_end         = it_aa + 5;
  auto it             = it_aa;
  double single_score = 0;
  while(it != it_end)
    {
      if(std::binary_search(spectrum_code_list.begin(),
                            spectrum_code_list.end(),
                            (std::uint32_t)(*it)))
        {
          single_score += 1;
        }
      else
        {
          single_score += 0.1;
        }
      it++;
    }
 
  // duo
  auto itduo_end   = it_couple + 4;
  auto itduo       = it_couple;
  double duo_score = 0;
  while(itduo != itduo_end)
    {
      if(std::binary_search(
           spectrum_code_list.begin(), spectrum_code_list.end(), *itduo))
        {
          duo_score += 3;
        }
      else
        {
          duo_score += 0.2;
        }
      itduo++;
    }
 
 
  // trio
  auto it3_end      = it_trio + 3;
  auto it3          = it_trio;
  double trio_score = 0;
  while(it3 != it3_end)
    {
      if(std::binary_search(
           spectrum_code_list.begin(), spectrum_code_list.end(), *it3))
        {
          trio_score += 6;
        }
      else
        {
          trio_score += 0.5;
        }
      it3++;
    }
 
 
  // quatro
  auto it4_end        = it_quatro + 2;
  auto it4            = it_quatro;
  double quatro_score = 0;
  while(it4 != it4_end)
    {
      // element < value
      if(std::binary_search(
           spectrum_code_list.begin(), spectrum_code_list.end(), *it4))
        {
          quatro_score += 8;
        }
      quatro_score += 0.8;
      it4++;
    }
 
  // cinqo
  if(std::binary_search(
       spectrum_code_list.begin(), spectrum_code_list.end(), *it_cinqo))
    {
      score += 10;
    }
  else
    {
      score += 1;
    }
  return score * single_score * duo_score * trio_score * quatro_score;
}