aamodification.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/amino_acid/aamodification.h
 * \date 7/3/2015
 * \author Olivier Langella
 * \brief amino acid modification model
 */
 
/*******************************************************************************
 * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     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/>.
 *
 * Contributors:
 *     Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
 *implementation
 ******************************************************************************/
 
#include <QRegExp>
#include <QDebug>
#include <cmath>
 
#include "aamodification.h"
#include "aa.h"
#include "../pappsoexception.h"
#include "../mzrange.h"
#include "../peptide/peptide.h"
#include "../obo/filterobopsimodsink.h"
#include "../obo/filterobopsimodtermaccession.h"
#include "../exception/exceptionnotfound.h"
 
/*
 
inline void initMyResource() {
    Q_INIT_RESOURCE(resources);
}
*/
 
namespace pappso
{
 
QMutex AaModification::m_mutex;
 
AaModification::AaModification(const QString &accession, pappso_double mass)
  : m_mass(mass), m_accession(accession)
{
  m_atomCount = {{AtomIsotopeSurvey::C, 0},
                 {AtomIsotopeSurvey::H, 0},
                 {AtomIsotopeSurvey::N, 0},
                 {AtomIsotopeSurvey::O, 0},
                 {AtomIsotopeSurvey::S, 0}};
 
  m_mapIsotope = {{Isotope::C13, 0},
                  {Isotope::H2, 0},
                  {Isotope::N15, 0},
                  {Isotope::O17, 0},
                  {Isotope::O18, 0},
                  {Isotope::S33, 0},
                  {Isotope::S34, 0},
                  {Isotope::S36, 0}};
}
 
 
AaModification::AaModification(AaModification &&toCopy) // move constructor
  : m_mass(toCopy.m_mass),
    m_atomCount(std::move(toCopy.m_atomCount)),
    m_mapIsotope(toCopy.m_mapIsotope),
    m_accession(toCopy.m_accession)
{
  m_origin = toCopy.m_origin;
}
 
AaModification::~AaModification()
{
}
 
const QString &
AaModification::getAccession() const
{
 
  qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__;
  return m_accession;
}
const QString &
AaModification::getName() const
{
  return m_name;
}
AaModification::MapAccessionModifications
  AaModification::m_mapAccessionModifications = [] {
    MapAccessionModifications ret;
 
    return ret;
  }();
 
AaModificationP
AaModification::createInstance(const OboPsiModTerm &term)
{
  AaModification *new_mod;
  // qDebug() << " AaModification::createInstance begin";
  new_mod = new AaModification(term.m_accession, term.m_diffMono);
  // xref: DiffFormula: "C 0 H 0 N 0 O 1 S 0"
  new_mod->setDiffFormula(term.m_diffFormula);
  new_mod->setXrefOrigin(term.m_origin);
  new_mod->m_name = term.m_name;
  return new_mod;
}
 
AaModificationP
AaModification::createInstance(const QString &accession)
{
  if(accession == "internal:Nter_hydrolytic_cleavage_H")
    {
      OboPsiModTerm term;
      term.m_accession   = accession;
      term.m_diffFormula = "H 1";
      term.m_diffMono    = MPROTIUM;
      term.m_name        = "Nter hydrolytic cleavage H+";
      return (AaModification::createInstance(term));
    }
  if(accession == "internal:Cter_hydrolytic_cleavage_HO")
    {
      OboPsiModTerm term;
      term.m_accession   = accession;
      term.m_diffFormula = "H 1 O 1";
      term.m_diffMono    = MPROTIUM + MASSOXYGEN;
      term.m_name        = "Cter hydrolytic cleavage HO";
      return (AaModification::createInstance(term));
    }
  if(accession.startsWith("MUTATION:"))
    {
      QRegExp regexp_mutation("^MUTATION:([A-Z])=>([A-Z])$");
      if(regexp_mutation.exactMatch(accession))
        {
          qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__ << " "
                   << regexp_mutation.capturedTexts()[1].at(0) << " "
                   << regexp_mutation.capturedTexts()[2].at(0);
 
          Aa aa_from(
            regexp_mutation.capturedTexts()[1].toStdString().c_str()[0]);
          Aa aa_to(regexp_mutation.capturedTexts()[2].toStdString().c_str()[0]);
          AaModificationP instance_mutation =
            createInstanceMutation(aa_from, aa_to);
          return instance_mutation;
          // m_psiModLabel<<"|";
        }
    }
  // initMyResource();
  FilterOboPsiModSink term_list;
  FilterOboPsiModTermAccession filterm_accession(term_list, accession);
 
  OboPsiMod psimod(filterm_accession);
 
  try
    {
      return (AaModification::createInstance(term_list.getOne()));
    }
  catch(ExceptionNotFound &e)
    {
      throw ExceptionNotFound(QObject::tr("modification not found : [%1]\n%2")
                                .arg(accession)
                                .arg(e.qwhat()));
    }
}
 
void
AaModification::setXrefOrigin(const QString &origin)
{
  // xref: Origin: "N"
  // xref: Origin: "X"
  m_origin = origin;
}
void
AaModification::setDiffFormula(const QString &diff_formula)
{
  QRegExp rx("(^|\\s)([C,H,O,N,H,S])\\s([-]{0,1}\\d+)");
  int pos = 0;
 
  while((pos = rx.indexIn(diff_formula, pos)) != -1)
    {
      qDebug() << rx.cap(2) << " " << rx.cap(3);
      if(rx.cap(2) == "C")
        {
          m_atomCount[AtomIsotopeSurvey::C] = rx.cap(3).toInt();
        }
      else if(rx.cap(2) == "H")
        {
          m_atomCount[AtomIsotopeSurvey::H] = rx.cap(3).toInt();
        }
      else if(rx.cap(2) == "N")
        {
          m_atomCount[AtomIsotopeSurvey::N] = rx.cap(3).toInt();
        }
      else if(rx.cap(2) == "O")
        {
          m_atomCount[AtomIsotopeSurvey::O] = rx.cap(3).toInt();
        }
      else if(rx.cap(2) == "S")
        {
          m_atomCount[AtomIsotopeSurvey::S] = rx.cap(3).toInt();
        }
      pos += rx.matchedLength();
    }
 
  // look for isotopes :
  rx.setPattern("\\(([-]{0,1}\\d+)\\)([C,H,O,N,H,S])\\s([-]{0,1}\\d+)");
  pos = 0;
 
  while((pos = rx.indexIn(diff_formula, pos)) != -1)
    {
      qDebug() << rx.cap(1) << " " << rx.cap(2) << " " << rx.cap(3);
      int number_of_isotopes = rx.cap(3).toInt();
      if(rx.cap(2) == "C")
        {
          if(rx.cap(1) == "13")
            {
              m_mapIsotope.at(Isotope::C13) = number_of_isotopes;
            }
          m_atomCount[AtomIsotopeSurvey::C] += number_of_isotopes;
        }
      else if(rx.cap(2) == "H")
        {
          if(rx.cap(1) == "2")
            {
              m_mapIsotope.at(Isotope::H2) = number_of_isotopes;
            }
          m_atomCount[AtomIsotopeSurvey::H] += number_of_isotopes;
        }
      else if(rx.cap(2) == "N")
        {
          if(rx.cap(1) == "15")
            {
              m_mapIsotope.at(Isotope::N15) = number_of_isotopes;
            }
          m_atomCount[AtomIsotopeSurvey::N] += number_of_isotopes;
        }
      else if(rx.cap(2) == "O")
        {
          if(rx.cap(1) == "17")
            {
              m_mapIsotope.at(Isotope::O17) = number_of_isotopes;
            }
          else if(rx.cap(1) == "18")
            {
              m_mapIsotope.at(Isotope::O18) = number_of_isotopes;
            }
          m_atomCount[AtomIsotopeSurvey::O] += number_of_isotopes;
        }
      else if(rx.cap(2) == "S")
        {
          if(rx.cap(1) == "33")
            {
              m_mapIsotope.at(Isotope::S33) = number_of_isotopes;
            }
          else if(rx.cap(1) == "34")
            {
              m_mapIsotope.at(Isotope::S34) = number_of_isotopes;
            }
          else if(rx.cap(1) == "36")
            {
              m_mapIsotope.at(Isotope::S36) = number_of_isotopes;
            }
          m_atomCount[AtomIsotopeSurvey::S] += number_of_isotopes;
        }
      pos += rx.matchedLength();
    }
 
 
  calculateMassFromChemicalComponents();
}
 
 
void
AaModification::calculateMassFromChemicalComponents()
{
  pappso_double theoreticalm_mass = 0;
  std::map<AtomIsotopeSurvey, int>::const_iterator it_atom =
    m_atomCount.find(AtomIsotopeSurvey::C);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSCARBON * (it_atom->second);
    }
  it_atom = m_atomCount.find(AtomIsotopeSurvey::H);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MPROTIUM * (it_atom->second);
    }
 
  it_atom = m_atomCount.find(AtomIsotopeSurvey::O);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSOXYGEN * (it_atom->second);
    }
 
  it_atom = m_atomCount.find(AtomIsotopeSurvey::N);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSNITROGEN * (it_atom->second);
    }
  it_atom = m_atomCount.find(AtomIsotopeSurvey::S);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSSULFUR * (it_atom->second);
    }
 
  qDebug() << theoreticalm_mass;
 
  theoreticalm_mass += DIFFC12C13 * m_mapIsotope.at(Isotope::C13);
  theoreticalm_mass += DIFFH1H2 * m_mapIsotope.at(Isotope::H2);
  theoreticalm_mass += DIFFN14N15 * m_mapIsotope.at(Isotope::N15);
  theoreticalm_mass += DIFFO16O17 * m_mapIsotope.at(Isotope::O17);
  theoreticalm_mass += DIFFO16O18 * m_mapIsotope.at(Isotope::O18);
  theoreticalm_mass += DIFFS32S33 * m_mapIsotope.at(Isotope::S33);
  theoreticalm_mass += DIFFS32S34 * m_mapIsotope.at(Isotope::S34);
  theoreticalm_mass += DIFFS32S36 * m_mapIsotope.at(Isotope::S36);
 
 
  pappso_double diff = std::fabs((pappso_double)m_mass - theoreticalm_mass);
  if(diff < 0.001)
    {
      m_mass = theoreticalm_mass;
      qDebug() << "AaModification::calculateMassFromChemicalComponents "
               << diff;
    }
  else
    {
      qDebug()
        << "ERROR in AaModification::calculateMassFromChemicalComponents theo="
        << theoreticalm_mass << " m=" << m_mass << " diff=" << diff
        << " accession=" << m_accession;
    }
}
 
AaModificationP
AaModification::getInstanceCustomizedMod(pappso_double modificationMass)
{
  QString accession = QString("%1").arg(modificationMass);
  qDebug() << "AaModification::getInstanceCustomizedMod " << accession;
  QMutexLocker locker(&m_mutex);
  if(m_mapAccessionModifications.find(accession) ==
     m_mapAccessionModifications.end())
    {
      // not found
      m_mapAccessionModifications.insert(std::pair<QString, AaModification *>(
        accession, new AaModification(accession, modificationMass)));
    }
  else
    {
      // found
    }
  return m_mapAccessionModifications.at(accession);
}
 
AaModificationP
AaModification::getInstance(const QString &accession)
{
  try
    {
      QMutexLocker locker(&m_mutex);
      MapAccessionModifications::iterator it =
        m_mapAccessionModifications.find(accession);
      if(it == m_mapAccessionModifications.end())
        {
 
          // not found
          std::pair<MapAccessionModifications::iterator, bool> insert_res =
            m_mapAccessionModifications.insert(
              std::pair<QString, AaModificationP>(
                accession, AaModification::createInstance(accession)));
          it = insert_res.first;
        }
      else
        {
          // found
        }
      return it->second;
    }
  catch(ExceptionNotFound &e)
    {
      throw ExceptionNotFound(
        QObject::tr("ERROR getting instance of : %1 NOT FOUND\n%2")
          .arg(accession)
          .arg(e.qwhat()));
    }
  catch(PappsoException &e)
    {
      throw PappsoException(QObject::tr("ERROR getting instance of %1\n%2")
                              .arg(accession)
                              .arg(e.qwhat()));
    }
  catch(std::exception &e)
    {
      throw PappsoException(QObject::tr("ERROR getting instance of %1\n%2")
                              .arg(accession)
                              .arg(e.what()));
    }
}
 
AaModificationP
AaModification::getInstance(const OboPsiModTerm &oboterm)
{
 
  QMutexLocker locker(&m_mutex);
  MapAccessionModifications::iterator it =
    m_mapAccessionModifications.find(oboterm.m_accession);
  if(it == m_mapAccessionModifications.end())
    {
      // not found
      std::pair<MapAccessionModifications::iterator, bool> insert_res =
        m_mapAccessionModifications.insert(std::pair<QString, AaModificationP>(
          oboterm.m_accession, AaModification::createInstance(oboterm)));
      it = insert_res.first;
    }
  else
    {
      // found
    }
  return it->second;
}
 
 
AaModificationP
AaModification::getInstanceXtandemMod(const QString &type,
                                      pappso_double mass,
                                      const PeptideSp &peptide_sp,
                                      unsigned int position)
{
  PrecisionPtr precision = PrecisionFactory::getDaltonInstance(0.001);
  if(MzRange(mass, precision).contains(getInstance("MOD:00719")->getMass()))
    {
      if(type == "M")
        {
          return getInstance("MOD:00719");
        }
      if(type == "K")
        {
          return getInstance("MOD:01047");
        }
    }
  // accession== "MOD:00057"
  if(MzRange(mass, precision).contains(getInstance("MOD:00408")->getMass()))
    {
      // id: MOD:00394
      // name: acetylated residue
      // potential N-terminus modifications
      if(position == 0)
        {
          return getInstance("MOD:00408");
        }
    }
  if(MzRange(mass, precision).contains(getInstance("MOD:01160")->getMass()))
    {
      //-17.02655
      // loss of ammonia [MOD:01160] -17.026549
      return getInstance("MOD:01160");
    }
 
  if(MzRange(mass, precision).contains(getInstance("MOD:01060")->getMass()))
    {
      //// iodoacetamide [MOD:00397] 57.021464
      if(type == "C")
        {
          return getInstance("MOD:01060");
        }
      else
        {
          return getInstance("MOD:00397");
        }
    }
  if(MzRange(mass, precision).contains(getInstance("MOD:00704")->getMass()))
    {
      // loss of water
      /*
        if (position == 0) {
            if (peptide_sp.get()->getSequence().startsWith("EG")) {
                return getInstance("MOD:00365");
            }
            if (peptide_sp.get()->getSequence().startsWith("ES")) {
                return getInstance("MOD:00953");
            }
            if (type == "E") {
                return getInstance("MOD:00420");
            }
        }
      */
      // dehydrated residue [MOD:00704] -18.010565
      return getInstance("MOD:00704");
    }
  if(MzRange(mass, precision).contains(getInstance("MOD:00696")->getMass()))
    {
      // phosphorylated residue [MOD:00696] 79.966330
      return getInstance("MOD:00696");
    }
  bool isCter = false;
  if(peptide_sp.get()->size() == (position + 1))
    {
      isCter = true;
    }
  if((position == 0) || isCter)
    {
      if(MzRange(mass, precision).contains(getInstance("MOD:00429")->getMass()))
        {
          // dimethyl
          return getInstance("MOD:00429");
        }
      if(MzRange(mass, precision).contains(getInstance("MOD:00552")->getMass()))
        {
          // 4x(2)H labeled dimethyl residue
          return getInstance("MOD:00552");
        }
      if(MzRange(mass, precision).contains(getInstance("MOD:00638")->getMass()))
        {
          // 2x(13)C,6x(2)H-dimethylated arginine
          return getInstance("MOD:00638");
        }
    }
  throw PappsoException(
    QObject::tr("tandem modification not found : %1 %2 %3 %4")
      .arg(type)
      .arg(mass)
      .arg(peptide_sp.get()->getSequence())
      .arg(position));
}
 
pappso_double
AaModification::getMass() const
{
  return m_mass;
}
 
 
int
AaModification::getNumberOfAtom(AtomIsotopeSurvey atom) const
{
  // qDebug() << "AaModification::getNumberOfAtom(AtomIsotopeSurvey atom) NOT
  // IMPLEMENTED";
  return m_atomCount.at(atom);
}
 
 
int
AaModification::getNumberOfIsotope(Isotope isotope) const
{
  try
    {
      return m_mapIsotope.at(isotope);
    }
  catch(std::exception &e)
    {
      throw PappsoException(
        QObject::tr("ERROR in AaModification::getNumberOfIsotope %2")
          .arg(e.what()));
    }
}
 
 
bool
AaModification::isInternal() const
{
  if(m_accession.startsWith("internal:"))
    {
      return true;
    }
  return false;
}
 
AaModificationP
AaModification::createInstanceMutation(const Aa &aa_from, const Aa &aa_to)
{
  QString accession(
    QString("MUTATION:%1=>%2").arg(aa_from.getLetter()).arg(aa_to.getLetter()));
  double diffMono = aa_to.getMass() - aa_from.getMass();
  // not found
  AaModification *instance_mutation;
  // qDebug() << " AaModification::createInstance begin";
  instance_mutation = new AaModification(accession, diffMono);
  // xref: DiffFormula: "C 0 H 0 N 0 O 1 S 0"
 
  for(std::int8_t atomInt = (std::int8_t)AtomIsotopeSurvey::C;
      atomInt != (std::int8_t)AtomIsotopeSurvey::last;
      atomInt++)
    {
      AtomIsotopeSurvey atom = static_cast<AtomIsotopeSurvey>(atomInt);
      instance_mutation->m_atomCount[atom] =
        aa_to.getNumberOfAtom(atom) - aa_from.getNumberOfAtom(atom);
    }
  instance_mutation->m_name = QString("mutation from %1 to %2")
                                .arg(aa_from.getLetter())
                                .arg(aa_to.getLetter());
  return instance_mutation;
}
 
 
AaModificationP
AaModification::getInstanceMutation(const QChar &mut_from, const QChar &mut_to)
{
  QString accession(QString("MUTATION:%1=>%2").arg(mut_from).arg(mut_to));
  try
    {
      QMutexLocker locker(&m_mutex);
      MapAccessionModifications::iterator it =
        m_mapAccessionModifications.find(accession);
      if(it == m_mapAccessionModifications.end())
        {
          Aa aa_from(mut_from.toLatin1());
          Aa aa_to(mut_to.toLatin1());
          AaModificationP instance_mutation =
            createInstanceMutation(aa_from, aa_to);
 
          std::pair<MapAccessionModifications::iterator, bool> insert_res =
            m_mapAccessionModifications.insert(
              std::pair<QString, AaModificationP>(accession,
                                                  instance_mutation));
          it = insert_res.first;
        }
      else
        {
          // found
        }
      return it->second;
    }
  catch(ExceptionNotFound &e)
    {
      throw ExceptionNotFound(
        QObject::tr("ERROR getting instance of : %1 NOT FOUND\n%2")
          .arg(accession)
          .arg(e.qwhat()));
    }
  catch(PappsoException &e)
    {
      throw PappsoException(QObject::tr("ERROR getting instance of %1\n%2")
                              .arg(accession)
                              .arg(e.qwhat()));
    }
  catch(std::exception &e)
    {
      throw PappsoException(QObject::tr("ERROR getting instance of %1\n%2")
                              .arg(accession)
                              .arg(e.what()));
    }
} // namespace pappso
 
} // namespace pappso