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 <QRegularExpression>
#include <QDebug>
#include <cmath>
 
#include "pappsomspp/core/amino_acid/aamodification.h"
#include "pappsomspp/core/amino_acid/aa.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/mzrange.h"
#include "../peptide/peptide.h"
#include "../obo/filterobopsimodsink.h"
#include "../obo/filterobopsimodtermaccession.h"
#include "pappsomspp/core/exception/exceptionnotfound.h"
#include "../obo/filterobopsimodmap.h"
#include "../obo/obounimod.h"
 
/*
 
inline void initMyResource() {
    Q_INIT_RESOURCE(resources);
}
*/
 
namespace pappso
{
 
QMutex AaModification::m_mutex;
 
AaModification::AaModification(const QString &accession, pappso_double mass)
  : m_accession(accession), m_mass(mass)
{
  m_atomCount = {{Enums::AtomIsotopeSurvey::C, 0},
                 {Enums::AtomIsotopeSurvey::H, 0},
                 {Enums::AtomIsotopeSurvey::N, 0},
                 {Enums::AtomIsotopeSurvey::O, 0},
                 {Enums::AtomIsotopeSurvey::S, 0},
                 {Enums::AtomIsotopeSurvey::P, 0}};
 
  m_mapIsotope = {{Enums::Isotope::C13, 0},
                  {Enums::Isotope::H2, 0},
                  {Enums::Isotope::N15, 0},
                  {Enums::Isotope::O17, 0},
                  {Enums::Isotope::O18, 0},
                  {Enums::Isotope::S33, 0},
                  {Enums::Isotope::S34, 0},
                  {Enums::Isotope::S36, 0}};
}
 
 
AaModification::AaModification(AaModification &&toCopy) // move constructor
  : m_accession(toCopy.m_accession),
    m_name(toCopy.m_name),
    m_mass(toCopy.m_mass),
    m_atomCount(std::move(toCopy.m_atomCount)),
    m_mapIsotope(toCopy.m_mapIsotope)
{
  m_origin = toCopy.m_origin;
}
 
AaModification::~AaModification()
{
}
 
const QString &
AaModification::getAccession() const
{
 
  // qDebug();
  return m_accession;
}
 
 
const QString
AaModification::toProForma() const
{
  QString str_return;
  // qDebug();
  if(m_accession.startsWith("MOD:"))
    return m_accession;
  if(m_accession.startsWith("UNIMOD:"))
    return m_accession;
 
  if(m_accession.startsWith("internal:"))
    return str_return;
 
  if(m_accession.startsWith("MUTATION:"))
    {
      // MUTATION:%1=>%2
      QStringList aa_list = m_accession.mid(9).split("=>");
      // first, find removal accession
      return QString("%1][%2")
        .arg(getInstanceRemovalAccessionByAaLetter(aa_list[0].at(0))->getAccession())
        .arg(getInstanceInsertionAccessionByAaLetter(aa_list[1].at(0))->getAccession());
    }
 
  if(m_mass > 0)
    {
      return QString("+%1").arg(QString::number(m_mass, 'f', 4));
    }
  else
    {
      return QString("%1").arg(QString::number(m_mass, 'f', 4));
    }
  return str_return;
}
 
const QString &
AaModification::getName() const
{
  return m_name;
}
 
const QString &
AaModification::getXrefOrigin() const
{
  return m_origin;
}
 
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.getAccession(), term.m_diffMono);
  // xref: DiffFormula: "C 0 H 0 N 0 O 1 S 0"
  pappso::ChemicalFormula formula;
  formula.setOboPsiModTerm(term);
  new_mod->setDiffFormula(formula);
  new_mod->setXrefOrigin(term.m_origin);
  new_mod->m_name = term.m_name;
 
  qDebug() << new_mod->m_mass;
  return new_mod;
}
 
AaModificationP
AaModification::createInstance(const QString &accession)
{
  if(accession == "internal:Nter_hydrolytic_cleavage_H")
    {
      OboPsiModTerm term;
      term.setAccession(accession);
      term.m_diffFormula = "H 1";
      term.m_diffMono    = MPROTIUM;
      term.m_name        = "Nter hydrolytic cleavage H+";
      return (AaModification::createInstance(term));
    }
  else if(accession == "internal:Cter_hydrolytic_cleavage_HO")
    {
      OboPsiModTerm term;
      term.setAccession(accession);
      term.m_diffFormula = "H 1 O 1";
      term.m_diffMono    = MPROTIUM + MASSOXYGEN;
      term.m_name        = "Cter hydrolytic cleavage HO";
      return (AaModification::createInstance(term));
    }
  else if(accession == "USER:2-Ethynylbenzaldehyde")
    {
      // https://pubchem.ncbi.nlm.nih.gov/compound/640368
      OboPsiModTerm term;
      term.setAccession(accession);
      term.m_diffFormula = "H 4 C 9 O 0"; // C9H6O
      term.m_name        = "2-Ethynylbenzaldehyde";
      ChemicalFormula formula;
      formula.setOboPsiModTerm(term);
      term.m_diffMono = formula.getMass();
 
      return (AaModification::createInstance(term));
    }
  else if(accession.startsWith("MUTATION:"))
    {
      QRegularExpression regexp_mutation("^MUTATION:([A-Z])=>([A-Z])$");
      QRegularExpressionMatch match = regexp_mutation.match(accession);
      if(match.hasMatch())
        {
          qDebug() << match.capturedTexts()[1].at(0) << " " << match.capturedTexts()[2].at(0);
 
          Aa aa_from(match.capturedTexts()[1].toStdString().c_str()[0]);
          Aa aa_to(match.capturedTexts()[2].toStdString().c_str()[0]);
          AaModificationP instance_mutation = createInstanceMutation(aa_from, aa_to);
          return instance_mutation;
          // m_psiModLabel<<"|";
        }
    }
  else if(accession.startsWith("C13N15:"))
    {
      QRegularExpression regexp_label("^C13N15:([A-Z])$");
      QRegularExpressionMatch match = regexp_label.match(accession);
      if(match.hasMatch())
        {
          char labelled_aa(match.capturedTexts()[1].toStdString().c_str()[0]);
 
          AaModificationP instance_labelled_aa_p =
            createInstanceC13N15LabelledAminoAcid(QChar(labelled_aa));
 
 
          return instance_labelled_aa_p;
          // 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)
    {
      pappso::FilterOboPsiModMap obo_unimod_map;
      pappso::OboUnimod unimod_parser(obo_unimod_map);
      // OboUnimod unimod(filterm_accession);
 
      try
        {
          return (AaModification::createInstance(
            obo_unimod_map.getOboPsiModTermWithAccession(accession)));
        }
      catch(ExceptionNotFound &e)
        {
          try
            {
              return (
                AaModification::createInstance(obo_unimod_map.getOboPsiModTermWithName(accession)));
            }
          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::calculateMassFromChemicalComponents()
{
  pappso_double theoreticalm_mass = 0;
  std::map<Enums::AtomIsotopeSurvey, int>::const_iterator it_atom =
    m_atomCount.find(Enums::AtomIsotopeSurvey::C);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSCARBON * (it_atom->second);
    }
  it_atom = m_atomCount.find(Enums::AtomIsotopeSurvey::H);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MPROTIUM * (it_atom->second);
    }
 
  it_atom = m_atomCount.find(Enums::AtomIsotopeSurvey::O);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSOXYGEN * (it_atom->second);
    }
 
  it_atom = m_atomCount.find(Enums::AtomIsotopeSurvey::N);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSNITROGEN * (it_atom->second);
    }
  it_atom = m_atomCount.find(Enums::AtomIsotopeSurvey::S);
  if(it_atom != m_atomCount.end())
    {
      theoreticalm_mass += MASSSULFUR * (it_atom->second);
    }
 
  qDebug() << theoreticalm_mass;
 
  theoreticalm_mass += DIFFC12C13 * m_mapIsotope.at(Enums::Isotope::C13);
  theoreticalm_mass += DIFFH1H2 * m_mapIsotope.at(Enums::Isotope::H2);
  theoreticalm_mass += DIFFN14N15 * m_mapIsotope.at(Enums::Isotope::N15);
  theoreticalm_mass += DIFFO16O17 * m_mapIsotope.at(Enums::Isotope::O17);
  theoreticalm_mass += DIFFO16O18 * m_mapIsotope.at(Enums::Isotope::O18);
  theoreticalm_mass += DIFFS32S33 * m_mapIsotope.at(Enums::Isotope::S33);
  theoreticalm_mass += DIFFS32S34 * m_mapIsotope.at(Enums::Isotope::S34);
  theoreticalm_mass += DIFFS32S36 * m_mapIsotope.at(Enums::Isotope::S36);
 
 
  if(m_mass == 0.0)
    {
      m_mass = theoreticalm_mass;
    }
  else
    {
      pappso_double diff = std::fabs((pappso_double)m_mass - theoreticalm_mass);
      if(diff < 0.001)
        {
          m_mass = theoreticalm_mass;
          qDebug() << 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() << 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::getInstanceRemovalAccessionByAaLetter(const QChar &amino_acid)
{
  switch(amino_acid.toLatin1())
    {
      //    Enums::AminoAcidChar::alanine
      case 'A':
        return getInstance("MOD:01631");
        break;
        // Enums::AminoAcidChar::arginine,
      case 'R':
        return getInstance("MOD:01632");
        break;
        // Enums::AminoAcidChar::cysteine,
      case 'C':
        return getInstance("MOD:01635");
        break;
        // Enums::AminoAcidChar::aspartic_acid,
      case 'D':
        return getInstance("MOD:01634");
        break;
        // Enums::AminoAcidChar::glutamic_acid,
      case 'E':
        return getInstance("MOD:01636");
        break;
 
        // Enums::AminoAcidChar::phenylalanine,
      case 'F':
        return getInstance("MOD:01644");
        break;
        // Enums::AminoAcidChar::glycine,
      case 'G':
        return getInstance("MOD:01638");
        break;
        // Enums::AminoAcidChar::histidine,
      case 'H':
        return getInstance("MOD:01639");
        break;
 
        // Enums::AminoAcidChar::isoleucine,
      case 'I':
        return getInstance("MOD:01640");
        break;
 
        // Enums::AminoAcidChar::lysine,
      case 'K':
        return getInstance("MOD:01642");
        break;
        // Enums::AminoAcidChar::leucine,
      case 'L':
        return getInstance("MOD:01641");
        break;
 
        // Enums::AminoAcidChar::methionine,
      case 'M':
        return getInstance("MOD:01643");
        break;
 
        // Enums::AminoAcidChar::asparagine,
      case 'N':
        return getInstance("MOD:01633");
        break;
 
        // Enums::AminoAcidChar::proline,
      case 'P':
        return getInstance("MOD:01645");
        break;
 
        // Enums::AminoAcidChar::glutamine,
      case 'Q':
        return getInstance("MOD:01637");
        break;
        // Enums::AminoAcidChar::serine,
      case 'S':
        return getInstance("MOD:01646");
        break;
 
        // Enums::AminoAcidChar::threonine,
      case 'T':
        return getInstance("MOD:01647");
        break;
 
 
        // Enums::AminoAcidChar::valine,
      case 'V':
        return getInstance("MOD:01650");
        break;
 
        // Enums::AminoAcidChar::tryptophan,
      case 'W':
        return getInstance("MOD:01648");
        break;
 
        // Enums::AminoAcidChar::tyrosine,
      case 'Y':
        return getInstance("MOD:01649");
        break;
 
      default:
        throw ExceptionNotFound(
          QObject::tr("ERROR getting removal accession instance of amino acid: "
                      "%1 NOT FOUND")
            .arg(amino_acid));
    }
}
 
AaModificationP
AaModification::getInstanceInsertionAccessionByAaLetter(const QChar &amino_acid)
{
  pappso::AaModificationP aa_insert = nullptr;
  switch(amino_acid.toLatin1())
    {
      //    Enums::AminoAcidChar::alanine
      case 'A':
        aa_insert = getInstance("MOD:00010");
        break;
 
        // Enums::AminoAcidChar::arginine,
      case 'R':
        aa_insert = getInstance("MOD:00011");
        break;
 
        // Enums::AminoAcidChar::cysteine,
      case 'C':
        aa_insert = getInstance("MOD:00014");
        break;
 
        // Enums::AminoAcidChar::aspartic_acid,
      case 'D':
        aa_insert = getInstance("MOD:00013");
        break;
 
        // Enums::AminoAcidChar::glutamic_acid,
      case 'E':
        aa_insert = getInstance("MOD:00015");
        break;
 
        // Enums::AminoAcidChar::phenylalanine,
      case 'F':
        aa_insert = getInstance("MOD:00023");
        break;
 
        // Enums::AminoAcidChar::glycine,
      case 'G':
        aa_insert = getInstance("MOD:00017");
        break;
 
        // Enums::AminoAcidChar::histidine,
      case 'H':
        aa_insert = getInstance("MOD:00018");
        break;
 
        // Enums::AminoAcidChar::isoleucine,
      case 'I':
        aa_insert = getInstance("MOD:00019");
        break;
 
        // Enums::AminoAcidChar::lysine,
      case 'K':
        aa_insert = getInstance("MOD:00021");
        break;
        // Enums::AminoAcidChar::leucine,
      case 'L':
        aa_insert = getInstance("MOD:00020");
        break;
 
        // Enums::AminoAcidChar::methionine,
      case 'M':
        aa_insert = getInstance("MOD:00022");
        break;
 
        // Enums::AminoAcidChar::asparagine,
      case 'N':
        aa_insert = getInstance("MOD:00012");
        break;
 
        // Enums::AminoAcidChar::proline,
      case 'P':
        aa_insert = getInstance("MOD:00024");
        break;
 
        // Enums::AminoAcidChar::glutamine,
      case 'Q':
        aa_insert = getInstance("MOD:00016");
        break;
        // Enums::AminoAcidChar::serine,
      case 'S':
        aa_insert = getInstance("MOD:00025");
        break;
 
        // Enums::AminoAcidChar::threonine,
      case 'T':
        aa_insert = getInstance("MOD:00026");
        break;
 
 
        // Enums::AminoAcidChar::valine,
      case 'V':
        aa_insert = getInstance("MOD:00029");
        break;
 
        // Enums::AminoAcidChar::tryptophan,
      case 'W':
        aa_insert = getInstance("MOD:00027");
        break;
 
        // Enums::AminoAcidChar::tyrosine,
      case 'Y':
        aa_insert = getInstance("MOD:00028");
        break;
        // Enums::AminoAcidChar::selenocysteine,
      case 'U':
        aa_insert = getInstance("MOD:00031");
        break;
      default:
        throw ExceptionNotFound(
          QObject::tr("ERROR getting removal accession instance of amino acid: "
                      "%1 NOT FOUND")
            .arg(amino_acid));
    }
 
  return aa_insert;
}
 
AaModificationP
AaModification::getInstance(const OboPsiModTerm &oboterm)
{
 
  QMutexLocker locker(&m_mutex);
 
  std::pair<MapAccessionModifications::iterator, bool> insert_res =
    m_mapAccessionModifications.insert(
      std::pair<QString, AaModificationP>(oboterm.getAccession(), nullptr));
 
  if(!insert_res.second)
    {
      // Failed to insert, was there already by same key.
    }
  else
    {
      // Could insert (albeit with nullptr), was not there already by same key.
      // Just assign to the *correct* pair member the pointer to an
      // extemporaneously allocated AaModification.
      insert_res.first->second = AaModification::createInstance(oboterm);
    }
 
  return insert_res.first->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(Enums::AtomIsotopeSurvey atom) const
{
  // qDebug() << "AaModification::getNumberOfAtom(Enums::AtomIsotopeSurvey atom) NOT
  // IMPLEMENTED";
  return m_atomCount.at(atom);
}
 
 
int
AaModification::getNumberOfIsotope(Enums::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)Enums::AtomIsotopeSurvey::C;
      atomInt != (std::int8_t)Enums::AtomIsotopeSurvey::last;
      atomInt++)
    {
      Enums::AtomIsotopeSurvey atom = static_cast<Enums::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()));
    }
}
 
AaModificationP
AaModification::createInstanceC13N15LabelledAminoAcid(const QChar &from_aa)
{
  switch(from_aa.toLatin1())
    {
      case 'R':
        return createInstance("MOD:00587");
        break;
      case 'K':
        return createInstance("MOD:00582");
        break;
      case 'F':
        return createInstance("MOD:00589");
        break;
    }
 
  // At this point we have no MOD ontology term for all the other amino-acid
  // residue, we need to craft the stuff ourselves.
 
  QString accession(QString("C13N15:%1").arg(from_aa));
 
  // Define the diff formula for each amino acid code not handled above.
 
  QString diff_formula;
 
  switch(from_aa.toLatin1())
    {
      case 'G':
        diff_formula = "(12)C -2 (13)C 2 (14)N -1 (15)N 1";
        break;
      case 'A':
        diff_formula = "(12)C -3 (13)C 3 (14)N -1 (15)N 1";
        break;
      case 'V':
        diff_formula = "(12)C -5 (13)C 5 (14)N -1 (15)N 1";
        break;
      case 'I':
        diff_formula = "(12)C -6 (13)C 6 (14)N -1 (15)N 1";
        break;
      case 'L':
        diff_formula = "(12)C -6 (13)C 6 (14)N -1 (15)N 1";
        break;
      case 'S':
        diff_formula = "(12)C -3 (13)C 3 (14)N -1 (15)N 1";
        break;
      case 'T':
        diff_formula = "(12)C -4 (13)C 4 (14)N -1 (15)N 1";
        break;
      case 'C':
        diff_formula = "(12)C -3 (13)C 3 (14)N -1 (15)N 1";
        break;
      case 'M':
        diff_formula = "(12)C -5 (13)C 5 (14)N -1 (15)N 1";
        break;
      case 'D':
        diff_formula = "(12)C -4 (13)C 4 (14)N -1 (15)N 1";
        break;
      case 'E':
        diff_formula = "(12)C -5 (13)C 5 (14)N -1 (15)N 1";
        break;
      case 'N':
        diff_formula = "(12)C -4 (13)C 4 (14)N -2 (15)N 2";
        break;
      case 'Q':
        diff_formula = "(12)C -5 (13)C 5 (14)N -2 (15)N 2";
        break;
      case 'W':
        diff_formula = "(12)C -11 (13)C 11 (14)N -2 (15)N 2";
        break;
      case 'Y':
        diff_formula = "(12)C -9 (13)C 9 (14)N -1 (15)N 1";
        break;
      case 'H':
        diff_formula = "(12)C -6 (13)C 6 (14)N -3 (15)N 3";
        break;
      case 'P':
        diff_formula = "(12)C -5 (13)C 5 (14)N -1 (15)N 1";
        break;
      default:
        throw PappsoException(
          QObject::tr("ERROR creating C13N15-labelled amino acid residue %1\n").arg(from_aa));
    }
 
  OboPsiModTerm obo_psi_term;
  obo_psi_term.setAccession(accession);
  obo_psi_term.m_diffFormula = diff_formula;
  obo_psi_term.m_name        = QString("Fully C13N15-labelled %1").arg(from_aa);
  obo_psi_term.m_diffMono    = 0;
 
  qDebug() << "obo_psi_term.m_diffFormula: " << obo_psi_term.m_diffFormula;
  return AaModification::createInstance(obo_psi_term);
}
 
} // namespace pappso
 
void
pappso::AaModification::setDiffFormula(const pappso::ChemicalFormula &diff_formula)
{
  m_atomCount[Enums::AtomIsotopeSurvey::C] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::C);
  m_atomCount[Enums::AtomIsotopeSurvey::H] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::H);
  m_atomCount[Enums::AtomIsotopeSurvey::N] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::N);
  m_atomCount[Enums::AtomIsotopeSurvey::O] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::O);
  m_atomCount[Enums::AtomIsotopeSurvey::S] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::S);
  m_atomCount[Enums::AtomIsotopeSurvey::P] =
    diff_formula.getNumberOfAtom(Enums::AtomIsotopeSurvey::P);
 
  for(Enums::Isotope isotope : {Enums::Isotope::C13,
                                Enums::Isotope::H2,
                                Enums::Isotope::N15,
                                Enums::Isotope::O17,
                                Enums::Isotope::O18,
                                Enums::Isotope::S33,
                                Enums::Isotope::S34,
                                Enums::Isotope::S36})
    {
      m_mapIsotope.at(isotope) = diff_formula.getNumberOfIsotope(isotope);
    }
 
  calculateMassFromChemicalComponents();
}