proteinpresenceabsencematrix.cpp

Go to the documentation of this file.

/**
 * \file protein/proteinpresenceabsencematrix.cpp
 * \date 07/02/2024
 * \author Olivier Langella
 * \brief presence/absence matrix of amino acid code along the protein sequence
 */
 
/*******************************************************************************
 * Copyright (c) 2024 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 "proteinpresenceabsencematrix.h"
#include "proteinintegercode.h"
 
pappso::ProteinPresenceAbsenceMatrix::ProteinPresenceAbsenceMatrix()
{
}
 
pappso::ProteinPresenceAbsenceMatrix::~ProteinPresenceAbsenceMatrix()
{
}
 
 
pappso::ProteinPresenceAbsenceMatrix::ProteinPresenceAbsenceMatrix(
  const pappso::ProteinPresenceAbsenceMatrix &other)
{
  m_presenceAbsenceMatrix = other.m_presenceAbsenceMatrix;
}
 
pappso::ProteinPresenceAbsenceMatrix &
pappso::ProteinPresenceAbsenceMatrix::operator=(const pappso::ProteinPresenceAbsenceMatrix &other)
{
  m_presenceAbsenceMatrix = other.m_presenceAbsenceMatrix;
 
  return *this;
}
 
 
const boost::numeric::ublas::matrix<
  pappso::ProteinPresenceAbsenceMatrix::ProteinPresenceAbsenceMatrixElement> &
pappso::ProteinPresenceAbsenceMatrix::getPresenceAbsenceMatrix() const
{
  return m_presenceAbsenceMatrix;
}
 
void
pappso::ProteinPresenceAbsenceMatrix::fillMatrix(
  const pappso::ProteinIntegerCode &coded_protein,
  const std::vector<uint32_t> &code_list_from_spectrum)
{
  const std::vector<std::uint8_t> &seq_aa_code = coded_protein.getSeqAaCode();
 
  m_presenceAbsenceMatrix.resize(seq_aa_code.size(), 5);
  std::vector<std::uint8_t>::const_iterator it_aa = seq_aa_code.begin();
  auto it_couple = coded_protein.getPeptideCodedFragment(2).begin();
  auto it_trio   = coded_protein.getPeptideCodedFragment(3).begin();
  auto it_quatro = coded_protein.getPeptideCodedFragment(4).begin();
  auto it_cinqo  = coded_protein.getPeptideCodedFragment(5).begin();
 
  for(std::size_t i = 0; i < seq_aa_code.size(); i++)
    {
      if(std::binary_search(
           code_list_from_spectrum.begin(), code_list_from_spectrum.end(), (std::uint32_t)(*it_aa)))
        {
          // presence
          m_presenceAbsenceMatrix(i, 0) = ProteinPresenceAbsenceMatrixElement::present;
        }
      else
        {
          m_presenceAbsenceMatrix(i, 0) = ProteinPresenceAbsenceMatrixElement::absent;
        }
      if(std::binary_search(
           code_list_from_spectrum.begin(), code_list_from_spectrum.end(), *it_couple))
        {
          // presence
          m_presenceAbsenceMatrix(i, 1) = ProteinPresenceAbsenceMatrixElement::present;
        }
      else
        {
          m_presenceAbsenceMatrix(i, 1) = ProteinPresenceAbsenceMatrixElement::absent;
        }
      if(std::binary_search(
           code_list_from_spectrum.begin(), code_list_from_spectrum.end(), *it_trio))
        {
          // presence
          m_presenceAbsenceMatrix(i, 2) = ProteinPresenceAbsenceMatrixElement::present;
        }
      else
        {
          m_presenceAbsenceMatrix(i, 2) = ProteinPresenceAbsenceMatrixElement::absent;
        }
      if(std::binary_search(
           code_list_from_spectrum.begin(), code_list_from_spectrum.end(), *it_quatro))
        {
          // presence
          m_presenceAbsenceMatrix(i, 3) = ProteinPresenceAbsenceMatrixElement::present;
        }
      else
        {
          m_presenceAbsenceMatrix(i, 3) = ProteinPresenceAbsenceMatrixElement::absent;
        }
      if(std::binary_search(
           code_list_from_spectrum.begin(), code_list_from_spectrum.end(), *it_cinqo))
        {
          // presence
          m_presenceAbsenceMatrix(i, 4) = ProteinPresenceAbsenceMatrixElement::present;
        }
      else
        {
          m_presenceAbsenceMatrix(i, 4) = ProteinPresenceAbsenceMatrixElement::absent;
        }
      it_aa++;
      it_couple++;
      it_trio++;
      it_quatro++;
      it_cinqo++;
    }
}
 
 
std::vector<double>
pappso::ProteinPresenceAbsenceMatrix::convolution() const
{
 
  std::vector<double> convolution_score;
 
  for(std::size_t ipos = 0; ipos < m_presenceAbsenceMatrix.size1(); ipos++)
    {
      convolution_score.push_back(convolutionKernel(ipos));
    }
 
  return convolution_score;
}
 
double
pappso::ProteinPresenceAbsenceMatrix::getScore(std::size_t seq_position,
                                               std::size_t aa_fragment_size) const
{
  if(m_presenceAbsenceMatrix(seq_position, aa_fragment_size) ==
     ProteinPresenceAbsenceMatrixElement::present)
    {
      switch(aa_fragment_size)
        {
          case 0:
            return 1;
          case 1:
            return 3;
          case 2:
            return 6;
          case 3:
            return 8;
          case 4:
            return 10;
          default:
            break;
        }
    }
  else
    {
      // absent
      switch(aa_fragment_size)
        {
          case 0:
            return 0.1;
          case 1:
            return 0.3;
          case 2:
            return 0.6;
          case 3:
            return 0.8;
          case 4:
            return 1;
          default:
            break;
        }
    }
  return 0.1;
}
 
double
pappso::ProteinPresenceAbsenceMatrix::convolutionKernel(std::size_t position) const
{
  double score = 0;
 
  std::size_t size_seq = m_presenceAbsenceMatrix.size1();
  // single :
  double single_score = 0;
  std::size_t endpos  = std::min(position + 5, size_seq);
  for(std::size_t ipos = position; ipos < endpos; ipos++)
    {
      single_score += getScore(ipos, 0);
    }
 
  // duo
  double duo_score = 0;
  endpos           = std::min(position + 4, size_seq);
  for(std::size_t ipos = position; ipos < endpos; ipos++)
    {
      duo_score += getScore(ipos, 1);
    }
 
 
  // trio
  double trio_score = 0;
  endpos            = std::min(position + 3, size_seq);
  for(std::size_t ipos = position; ipos < endpos; ipos++)
    {
      trio_score += getScore(ipos, 2);
    }
 
 
  // quatro
  double quatro_score = 0;
  endpos              = std::min(position + 2, size_seq);
  for(std::size_t ipos = position; ipos < endpos; ipos++)
    {
      quatro_score += getScore(ipos, 3);
    }
 
  // cinqo
  score += getScore(position, 4);
  return score * single_score * duo_score * trio_score * quatro_score;
}