peptide.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/masschroq/peptide.cpp
 * \date 24/09/2024
 * \author Olivier Langella
 * \brief peptide model in masschroqlite
 */
 
/*******************************************************************************
 * Copyright (c) 2024 Olivier Langella
 *<Olivier.Langella@universite-paris-saclay.fr>.
 *
 * This file is part of MassChroQ.
 *
 *     MassChroQ 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.
 *
 *     MassChroQ 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 MassChroQ.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#include "peptide.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionnotfound.h"
#include "../peptide/peptideproformaparser.h"
#include "msrungroup.h"
#include "peptidelabel.h"
 
pappso::masschroq::Peptide::Peptide(const QString &id,
                                    const pappso::PeptideSp &peptide_sp,
                                    const std::vector<ProteinSp> &protein_list)
  : PeptideBase(peptide_sp), m_id(id), m_proteinSpList(protein_list)
{
}
 
pappso::masschroq::Peptide::Peptide(const Peptide &other)
  : PeptideBase(other.msp_peptide), m_id(other.m_id), m_proteinSpList(other.m_proteinSpList)
{
  m_mods = other.m_mods;
}
 
pappso::masschroq::Peptide::~Peptide()
{
  if(mpa_peptideNaturalIsotopeList != nullptr)
    delete mpa_peptideNaturalIsotopeList;
}
 
void
pappso::masschroq::Peptide::computeIsotopologues(double ni_min_abundance)
{
  qDebug();
  if(ni_min_abundance > 0)
    {
      if(m_peptideLabelMap.size() > 0)
        {
          qDebug();
          for(auto &pair_label : m_peptideLabelMap)
            {
              qDebug() << pair_label.second.get();
              pair_label.second->computeIsotopologues(ni_min_abundance);
              qDebug();
            }
          qDebug();
        }
      else
        {
          qDebug();
          PeptideBase::computeIsotopologues(ni_min_abundance);
        }
    }
  qDebug();
}
 
const QString &
pappso::masschroq::Peptide::getId() const
{
  return m_id;
}
 
void
pappso::masschroq::Peptide::setMods(const QString &mods)
{
  m_mods = mods;
}
 
const QString &
pappso::masschroq::Peptide::getMods() const
{
  return m_mods;
}
 
const std::map<QString, pappso::masschroq::PeptideLabelSp> &
pappso::masschroq::Peptide::getPeptideLabelMap() const
{
  return m_peptideLabelMap;
}
 
void
pappso::masschroq::Peptide::addObservedChargeState(std::uint8_t charge)
{
  if(m_allObservedChargeStateList.end() ==
     std::find(m_allObservedChargeStateList.begin(), m_allObservedChargeStateList.end(), charge))
    {
      m_allObservedChargeStateList.push_back(charge);
    }
}
 
const std::vector<std::uint8_t> &
pappso::masschroq::Peptide::getAllObservedChargeStateList() const
{
  return m_allObservedChargeStateList;
}
 
 
void
pappso::masschroq::Peptide::addObservedInMsRunSp(const pappso::masschroq::MsRunSp &msrun_sp)
{
  auto it = std::find(m_observedInMsRunSpList.begin(), m_observedInMsRunSpList.end(), msrun_sp);
  if(it == m_observedInMsRunSpList.end())
    {
      m_observedInMsRunSpList.push_back(msrun_sp);
    }
}
 
 
const std::vector<pappso::masschroq::ProteinSp> &
pappso::masschroq::Peptide::getProteinSpList() const
{
  return m_proteinSpList;
}
 
void
pappso::masschroq::Peptide::addAlignedPeakMeasurement(
  const PeptideMeasurements::Measurement &one_xic_measure,
  const pappso::MsRunRetentionTime<QString> &msrun_retention_time)
{
  if(one_xic_measure.m_tracePeakSp.get() != nullptr)
    {
      AlignedPeakPositionElement position;
      position.intensity = one_xic_measure.m_tracePeakSp.get()->getMaxXicElement().y;
      if(msrun_retention_time.isAligned())
        {
          position.alignedRetentionTime =
            msrun_retention_time.translateOriginal2AlignedRetentionTime(
              one_xic_measure.m_tracePeakSp.get()->getMaxXicElement().x);
          position.alignedRetentionTimeCenter =
            msrun_retention_time.translateOriginal2AlignedRetentionTime(
              (one_xic_measure.m_tracePeakSp.get()->getLeftBoundary().x +
               one_xic_measure.m_tracePeakSp.get()->getRightBoundary().x) /
              2);
        }
      else
        {
          position.alignedRetentionTime = one_xic_measure.m_tracePeakSp.get()->getMaxXicElement().x;
          position.alignedRetentionTimeCenter =
            ((one_xic_measure.m_tracePeakSp.get()->getLeftBoundary().x +
              one_xic_measure.m_tracePeakSp.get()->getRightBoundary().x) /
             2);
        }
      /*
    position.alignedRetentionTime =
      position.alignedRetentionTime +
      ((position.alignedRetentionTimeCenter - position.alignedRetentionTime) /
       4);*/
      position.wide = one_xic_measure.m_tracePeakSp.get()->getArea();
      //  one_xic_measure.m_tracePeakSp.get()->getRightBoundary().x -
      //  one_xic_measure.m_tracePeakSp.get()->getLeftBoundary().x;
      m_alignedPeakPositionElementList.push_back(position);
    }
}
 
void
pappso::masschroq::Peptide::addAlignedPeptideObservation(
  const pappso::masschroq::PeptideObservation &peptide_observation,
  const pappso::MsRunRetentionTime<QString> &msrun_retention_time)
{
  double best_rt = peptide_observation.getBestXicCoord().get()->rtTarget;
  if(msrun_retention_time.isAligned())
    {
      best_rt = msrun_retention_time.translateOriginal2AlignedRetentionTime(best_rt);
    }
 
  m_referenceMs2ObservationList.push_back(best_rt);
}
 
void
pappso::masschroq::Peptide::computeConsensusRetentionTime()
{
  m_consensusAlignedPeakRetentionTime = 0;
 
  if(m_alignedPeakPositionElementList.size() > 0)
    {
      std::sort(m_alignedPeakPositionElementList.begin(),
                m_alignedPeakPositionElementList.end(),
                [](AlignedPeakPositionElement const &a, AlignedPeakPositionElement const &b) {
                  return a.wide > b.wide;
                });
 
      // take only the 90% best measures
      std::size_t limit = m_alignedPeakPositionElementList.size() * 0.9;
      if(limit < 5)
        {
          limit = m_alignedPeakPositionElementList.size();
          //  throw pappso::PappsoException("limit < 3");
        }
      double total_intensity = std::accumulate(m_alignedPeakPositionElementList.begin(),
                                               m_alignedPeakPositionElementList.begin() + limit,
                                               (double)0.0,
                                               [](double sum, AlignedPeakPositionElement &element) {
                                                 sum += element.intensity;
                                                 return (sum);
                                               });
      double sum_rt          = std::accumulate(m_alignedPeakPositionElementList.begin(),
                                      m_alignedPeakPositionElementList.begin() + limit,
                                      (double)0.0,
                                      [](double sum, AlignedPeakPositionElement &element) {
                                        sum += element.intensity * element.alignedRetentionTime;
                                        return (sum);
                                      });
 
      m_consensusAlignedPeakRetentionTime = sum_rt / total_intensity;
    }
  m_alignedPeakPositionElementList.clear();
 
  double sum_rt = std::accumulate(
    m_referenceMs2ObservationList.begin(), m_referenceMs2ObservationList.end(), (double)0.0);
  m_consensusAlignedMs2RetentionTime = sum_rt / (double)m_referenceMs2ObservationList.size();
 
  m_referenceMs2ObservationList.clear();
}
 
bool
pappso::masschroq::Peptide::isObservedInMsRunSp(const pappso::masschroq::MsRun *msrun_p)
{
  auto it =
    std::find_if(m_observedInMsRunSpList.begin(),
                 m_observedInMsRunSpList.end(),
                 [msrun_p](const pappso::masschroq::MsRunSp &a) { return (a.get() == msrun_p); });
  if(it == m_observedInMsRunSpList.end())
    return false;
  return true;
}
 
void
pappso::masschroq::Peptide::setReferenceMsRunRetentionTimePtr(
  const pappso::MsRunRetentionTime<QString> *msrun_retention_time_reference_p)
{
  mp_referenceMsRunRetentionTime = msrun_retention_time_reference_p;
}
 
double
pappso::masschroq::Peptide::getConsensusPeakRetentionTime() const
{
  return m_consensusAlignedPeakRetentionTime;
}
 
double
pappso::masschroq::Peptide::getConsensusMs2RetentionTime() const
{
  return m_consensusAlignedMs2RetentionTime;
}
 
 
void
pappso::masschroq::Peptide::addMsRunXicCoordCharge(
  const pappso::masschroq::Peptide::MsRunXicCoordCharge &msrun_xic_coord)
{
  m_msRunXicCoordChargeList.push_back(msrun_xic_coord);
}
 
 
pappso::XicCoordSPtr
pappso::masschroq::Peptide::getBestXicCoordSPtrForChargeInMsRunGroup(
  const pappso::masschroq::MsRunGroup *msrun_group_p, std::uint8_t charge) const
{
  // get XIC coordinates of the most intense precursor
  pappso::pappso_double intensity = -1;
  pappso::XicCoordSPtr best_xic_coord;
 
  for(auto &xic_coord_charge : m_msRunXicCoordChargeList)
    {
      if(msrun_group_p->contains(xic_coord_charge.msrun_p))
        {
          if(xic_coord_charge.charge == charge)
            {
              if(xic_coord_charge.intensity > intensity)
                {
                  intensity      = xic_coord_charge.intensity;
                  best_xic_coord = xic_coord_charge.xic_coord_sp;
                  qDebug();
                }
            }
        }
    }
 
  return best_xic_coord;
}
 
 
void
pappso::masschroq::Peptide::populateIonMobilityGrid(
  pappso::IonMobilityGrid *ion_mobility_grid_p) const
{
 
  for(auto &xic_coord_charge_a : m_msRunXicCoordChargeList)
    {
      // qInfo() << " observed_in_a";
      const pappso::MsRunId &msrun_id_a =
        *(xic_coord_charge_a.msrun_p->getMsRunReaderSPtr().get()->getMsRunId().get());
      pappso::XicCoordSPtr xic_coord_a = xic_coord_charge_a.xic_coord_sp;
 
      for(auto &xic_coord_charge_b : m_msRunXicCoordChargeList)
        {
          if(&xic_coord_charge_a == &xic_coord_charge_b)
            {
              break;
            }
          else
            {
 
              // qInfo() << " observed_in_b";
              const pappso::MsRunId &msrun_id_b =
                *(xic_coord_charge_b.msrun_p->getMsRunReaderSPtr().get()->getMsRunId().get());
              if(msrun_id_b == msrun_id_a)
                {
                }
              else
                {
                  pappso::XicCoordSPtr xic_coord_b = xic_coord_charge_b.xic_coord_sp;
 
                  if(xic_coord_charge_a.charge == xic_coord_charge_b.charge)
                    {
                      ion_mobility_grid_p->storeObservedIdentityBetween(
                        msrun_id_a, xic_coord_a.get(), msrun_id_b, xic_coord_b.get());
                    }
                }
            }
        }
    }
}
 
 
pappso::XicCoordSPtr
pappso::masschroq::Peptide::getBestIonMobilityXicCoordToExtractOverallMsRunGroup(
  const pappso::masschroq::MsRunGroup &msrun_group,
  const pappso::MsRunId &targeted_msrun,
  std::uint8_t charge) const
{
  pappso::XicCoordSPtr xic_coord_mean;
  std::size_t count                            = 0;
  pappso::IonMobilityGrid *ion_mobility_grid_p = msrun_group.getIonMobilityGridSp().get();
 
  for(auto &msrun_xic_coord_charge : m_msRunXicCoordChargeList)
    {
      if(msrun_group.contains(msrun_xic_coord_charge.msrun_p))
        {
 
          if(msrun_xic_coord_charge.charge == charge)
            {
              // compute the mean coordinate :
              pappso::XicCoordSPtr xic_coord_to_add = msrun_xic_coord_charge.xic_coord_sp;
 
              if(ion_mobility_grid_p != nullptr)
                {
                  // ion mobility coordinate translation for the targeted msrun
                  xic_coord_to_add = ion_mobility_grid_p->translateXicCoordFromTo(
                    *xic_coord_to_add.get(),
                    *(msrun_xic_coord_charge.msrun_p->getMsRunReaderSPtr()
                        .get()
                        ->getMsRunId()
                        .get()),
                    targeted_msrun);
                }
 
              if(xic_coord_mean.get() == nullptr)
                {
                  xic_coord_mean = xic_coord_to_add.get()->initializeAndClone();
                }
              else
                {
                  xic_coord_mean = xic_coord_mean.get()->addition(xic_coord_to_add);
                }
              count++;
              qDebug() << " xic_coord_mean.get()->toString=" << xic_coord_mean.get()->toString();
            }
        }
    }
  if(xic_coord_mean.get() != nullptr)
    {
      xic_coord_mean = xic_coord_mean.get()->divideBy(count);
    }
 
 
  return xic_coord_mean;
}
 
void
pappso::masschroq::Peptide::setJsonLabelList(const QJsonObject &json_label_list)
{
  auto it = json_label_list.begin();
 
  while(it != json_label_list.end())
    {
      QString label = it.key();
 
      QString proforma = it.value().toObject().value("proforma").toString();
 
      pappso::masschroq::PeptideLabelSp peptide_label_sp;
 
      if(proforma.isEmpty())
        {
          throw pappso::ExceptionNotFound(
            QObject::tr("ERROR: label %1 %2 must contain proforma sequence").arg(label).arg(m_id));
        }
      else
        {
          peptide_label_sp = std::make_shared<pappso::masschroq::PeptideLabel>(
            pappso::PeptideProFormaParser::parseString(proforma), this, label);
        }
      m_peptideLabelMap.insert({label, peptide_label_sp});
      it++;
    }
}
 
pappso::masschroq::PeptideLabel *
pappso::masschroq::Peptide::getPeptideLabelPtr(const QString &label) const
{
  auto it = m_peptideLabelMap.find(label);
  if(it == m_peptideLabelMap.end())
    {
      throw pappso::ExceptionNotFound(
        QObject::tr("ERROR: label %1 not found in peptide %2").arg(label).arg(m_id));
    }
  return it->second.get();
}