xtandemhyperscorebis.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/psm/xtandem/xtandemhyperscorebis.cpp
 * \date 8/9/2016
 * \author Olivier Langella
 * \brief computation of the X!Tandem hyperscore
 */
 
/*******************************************************************************
 * 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 "xtandemhyperscorebis.h"
#include <cmath>
#include "../../peptide/peptiderawfragmentmasses.h"
#include "pappsomspp/core/pappsoexception.h"
 
 
// FIXME pourquoi ici cette directive et pas le namespace habituel ?
using namespace pappso;
 
unsigned int
factorial(unsigned int n)
{
  unsigned int retval = 1;
  for(int i = n; i > 1; --i)
    retval *= i;
  return retval;
}
 
 
XtandemHyperscoreBis::XtandemHyperscoreBis(bool refine_spectrum_synthesis,
                                           PrecisionPtr precision,
                                           const std::vector<Enums::PeptideIon> &ion_list)
{
 
  m_isRefineSpectrumSynthesis = refine_spectrum_synthesis;
  mp_precision                = precision;
  m_ionList                   = ion_list;
  m_totalMatchedIons          = 0;
}
 
XtandemHyperscoreBis::~XtandemHyperscoreBis()
{
}
 
unsigned int
XtandemHyperscoreBis::getMatchedIons(Enums::PeptideIon ion_type) const
{
  return m_ionCount[(std::int8_t)ion_type];
}
 
 
void
XtandemHyperscoreBis::reset()
{
 
  std::fill(m_ionCount, m_ionCount + PEPTIDE_ION_TYPE_COUNT, 0);
 
  m_totalMatchedIons = 0;
  m_protoHyperscore  = 0;
}
 
 
bool
XtandemHyperscoreBis::computeXtandemHyperscore(const MassSpectrum &spectrum,
                                               const Peptide &peptide,
                                               unsigned int parent_charge)
{
  try
    {
 
      unsigned int max_charge = parent_charge;
      if(parent_charge > 1)
        {
          max_charge = parent_charge - 1;
        }
      // PeptideSpectrumMatch psm (spectrum, peptideSp, max_charge, precision,
      // ion_list);
 
      PeptideRawFragmentMasses calc_mass_list(peptide, RawFragmentationMode::full);
 
      std::vector<SimplePeakIonMatch> match_products;
      unsigned int charge_i;
      for(Enums::PeptideIon &ion : m_ionList)
        {
          charge_i = max_charge;
          while(charge_i > 0)
            {
              calc_mass_list.pushBackMatchSpectrum(
                match_products, spectrum, mp_precision, ion, charge_i);
              charge_i--;
            }
        }
      m_totalMatchedIons = match_products.size();
      if(m_totalMatchedIons == 0)
        {
          m_protoHyperscore = 0;
        }
      else
        {
 
          unsigned int charge_ion_count[5][20] = {0};
 
          pappso_double charge_dot_product[5] = {0};
 
          QString sequence = peptide.getSequence();
 
          for(SimplePeakIonMatch &peptide_ion_match : match_products)
            {
              charge_dot_product[peptide_ion_match.ion_charge] +=
                peptide_ion_match.peak.y *
                getXtandemPredictedIonIntensityFactor(
                  sequence, peptide_ion_match.ion_type, peptide_ion_match.ion_size);
              charge_ion_count[peptide_ion_match.ion_charge]
                              [(std::int8_t)peptide_ion_match.ion_type] += 1;
              m_ionCount[(std::int8_t)peptide_ion_match.ion_type] += 1;
            }
          // take the 2 best component
          pappso_double sum_intensity = 0;
          for(unsigned int i = 1; i <= max_charge; i++)
            {
              sum_intensity += charge_dot_product[i];
            }
          for(auto count : m_ionCount)
            {
              sum_intensity *= factorial(count);
            }
 
          m_protoHyperscore = sum_intensity;
        }
      return true;
    }
  catch(PappsoException &exception_pappso)
    {
      QString errorStr = QObject::tr("ERROR computing hyperscore, PAPPSO exception:\n%1")
                           .arg(exception_pappso.qwhat());
      qDebug() << "XtandemHyperscore::XtandemHyperscore PappsoException :\n" << errorStr;
      throw PappsoException(errorStr);
    }
  catch(std::exception &exception_std)
    {
      QString errorStr =
        QObject::tr("ERROR computing hyperscore, std exception:\n%1").arg(exception_std.what());
      qDebug() << "XtandemHyperscore::XtandemHyperscore std::exception :\n" << errorStr;
      throw PappsoException(errorStr);
    }
}
 
 
XtandemHyperscoreBis::AaFactorMap XtandemHyperscoreBis::m_aaIonFactorY = [] {
  AaFactorMap ret;
  // populate ret
  for(long c = 64; c < 126; c++)
    {
      ret.insert(std::pair<char, pappso_double>(c, pappso_double(1.0)));
    }
  ret['P'] = pappso_double(5.0);
  return ret;
}();
 
 
XtandemHyperscoreBis::AaFactorMap XtandemHyperscoreBis::m_aaIonFactorBb = [] {
  AaFactorMap ret;
  // populate ret
  for(long c = 64; c < 126; c++)
    {
      ret.insert(std::pair<char, pappso_double>(c, pappso_double(1.0)));
    }
  ret['D'] = pappso_double(5.0);
  ret['N'] = pappso_double(2.0);
  ret['V'] = pappso_double(3.0);
  ret['E'] = pappso_double(3.0);
  ret['Q'] = pappso_double(2.0);
  ret['I'] = pappso_double(3.0);
  ret['L'] = pappso_double(3.0);
  return ret;
}();
 
 
unsigned int
XtandemHyperscoreBis::getXtandemPredictedIonIntensityFactor(const QString &sequence,
                                                            Enums::PeptideIon ion_type,
                                                            unsigned int ion_size) const
{
  unsigned int Pi(1);
 
  char last_aa_nter('_'), last_aa_cter('_');
 
  if(peptideIonIsNter(ion_type))
    {
      last_aa_nter = sequence[ion_size - 1].toLatin1();
      last_aa_cter = sequence[ion_size].toLatin1();
      if(ion_size == 2)
        {
          if(last_aa_nter == 'P')
            {
              Pi *= 10;
            }
          else
            {
              Pi *= 3;
            }
        }
    }
  else
    {
      unsigned int offset(sequence.size() - ion_size);
      last_aa_nter = sequence[offset - 1].toLatin1();
      last_aa_cter = sequence[offset].toLatin1();
      if((offset) == 2)
        {
          if(last_aa_nter == 'P')
            {
              Pi *= 10;
            }
          else
            {
              Pi *= 3;
            }
        }
    }
  // QDebug << " last_aa_nter=" << QChar(last_aa_nter) << "
  // m_aaIonFactorBb[last_aa_nter]="s ;
  // qDebug() << PeptideFragment::getPeptideIonDirectionName(ion_direction) << "
  // last_aa_nter=" << last_aa_nter << " m_aaIonFactorBb[last_aa_nter]=" <<
  // m_aaIonFactorBb[last_aa_nter] << " last_aa_cter=" << last_aa_cter << "
  // m_aaIonFactorY[last_aa_cter]=" << m_aaIonFactorY[last_aa_cter];
  if(m_isRefineSpectrumSynthesis)
    {
      Pi *= m_aaIonFactorBb[last_aa_nter] * m_aaIonFactorY[last_aa_cter];
    }
 
  return Pi;
}
 
unsigned int
XtandemHyperscoreBis::getTotalMatchedIons() const
{
  return m_totalMatchedIons;
}
 
pappso_double
XtandemHyperscoreBis::getHyperscore() const
{
  try
    {
      if(m_protoHyperscore == 0)
        return 0.0;
      return (log10(m_protoHyperscore) * 4);
    }
  catch(PappsoException &exception_pappso)
    {
      QString errorStr =
        QObject::tr("ERROR in getHyperscore, PAPPSO exception:\n%1").arg(exception_pappso.qwhat());
      qDebug() << "XtandemHyperscore::getHyperscore PappsoException :\n" << errorStr;
      throw PappsoException(errorStr);
    }
  catch(std::exception &exception_std)
    {
      QString errorStr =
        QObject::tr("ERROR in getHyperscore, std exception:\n%1").arg(exception_std.what());
      qDebug() << "XtandemHyperscore::getHyperscore std::exception :\n" << errorStr;
      throw PappsoException(errorStr);
    }
}