xtandemhyperscore.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/psm/xtandemhyperscore.cpp
 * \date 19/3/2015
 * \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 <QDebug>
#include <cmath>
#include "pappsomspp/core/pappsoexception.h"
#include "../../peptide/peptidefragment.h"
#include "../../peptide/peptidefragmentionlistbase.h"
#include "xtandemhyperscore.h"
#include "../peptidespectrummatch.h"
 
namespace pappso
{
 
unsigned int
factorial(unsigned int n)
{
  unsigned int retval = 1;
  for(int i = n; i > 1; --i)
    retval *= i;
  return retval;
}
 
XtandemHyperscore::XtandemHyperscore(const MassSpectrum &spectrum,
                                     pappso::PeptideSp peptideSp,
                                     unsigned int parent_charge,
                                     PrecisionPtr precision,
                                     std::list<Enums::PeptideIon> ion_list,
                                     bool refine_spectrum_synthesis)
  : _refine_spectrum_synthesis(refine_spectrum_synthesis)
{
  try
    {
      /*
      if ((peptide_ion_sp.get()->getPeptideIonType() == Enums::PeptideIon::c) ||
      (peptide_ion_sp.get()->getPeptideIonType() == Enums::PeptideIon::z)) {
                    if(current_max_charge > 2)  {
                        current_max_charge--;
                    }
                }*/
      unsigned int max_charge = parent_charge;
      if(parent_charge > 1)
        {
          max_charge = parent_charge - 1;
        }
      PeptideSpectrumMatch psm(
        spectrum, peptideSp, max_charge, precision, ion_list);
 
      _ion_count = {{Enums::PeptideIon::a, 0},
                    {Enums::PeptideIon::b, 0},
                    {Enums::PeptideIon::bo, 0},
                    {Enums::PeptideIon::bstar, 0},
                    {Enums::PeptideIon::c, 0},
                    {Enums::PeptideIon::y, 0},
                    {Enums::PeptideIon::yo, 0},
                    {Enums::PeptideIon::ystar, 0},
                    {Enums::PeptideIon::z, 0}};
 
      std::map<Enums::PeptideIon, unsigned int> ion_count;
      for(auto &&ion_type : ion_list)
        {
          ion_count.insert(std::pair<Enums::PeptideIon, unsigned int>(ion_type, 0));
        }
 
      std::map<unsigned int, pappso_double> charge_dot_product;
      std::map<unsigned int, std::map<Enums::PeptideIon, unsigned int>>
        charge_ion_count;
      for(unsigned int i = 1; i <= max_charge; i++)
        {
          charge_dot_product.insert(
            std::pair<unsigned int, pappso_double>(i, 0));
          charge_ion_count.insert(
            std::pair<unsigned int, std::map<Enums::PeptideIon, unsigned int>>(
              i, ion_count));
        }
      QString sequence = peptideSp.get()->getSequence();
      for(auto &&peptide_ion_match : psm)
        {
          Enums::PeptideIon ion_type = peptide_ion_match.getPeptideIonType();
          unsigned int charge = peptide_ion_match.getCharge();
          charge_dot_product[charge] +=
            peptide_ion_match.getPeak().y *
            getXtandemPredictedIonIntensityFactor(
              sequence,
              peptide_ion_match.getPeptideIonDirection(),
              peptide_ion_match.getPeptideFragmentIonSp().get()->size());
          charge_ion_count[charge][ion_type] += 1;
          _ion_count[ion_type] += 1;
        }
 
      qDebug() << __FILE__ << " " << __FUNCTION__ << " " << __LINE__
               << " _ion_count[Enums::PeptideIon::y]=" << _ion_count[Enums::PeptideIon::y];
      // 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 : _ion_count)
        {
          sum_intensity *= factorial(count.second);
        }
 
      _proto_hyperscore = sum_intensity;
    }
  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);
    }
}
 
XtandemHyperscore::XtandemHyperscore(const XtandemHyperscore &other)
{
  _ion_count                 = other._ion_count;
  _proto_hyperscore          = other._proto_hyperscore;
  _refine_spectrum_synthesis = other._refine_spectrum_synthesis;
}
 
unsigned int
XtandemHyperscore::getMatchedIons(Enums::PeptideIon ion_type) const
{
  return _ion_count.at(ion_type);
}
XtandemHyperscore::~XtandemHyperscore()
{
}
 
 
XtandemHyperscore::AaFactorMap XtandemHyperscore::_aa_ion_factor_y = [] {
  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;
}();
 
 
XtandemHyperscore::AaFactorMap XtandemHyperscore::_aa_ion_factor_b = [] {
  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
XtandemHyperscore::getXtandemPredictedIonIntensityFactor(
  const QString &sequence,
  PeptideDirection ion_direction,
  unsigned int ion_size) const
{
  unsigned int Pi(1);
 
  char last_aa_nter('_'), last_aa_cter('_');
 
  if(ion_direction == PeptideDirection::Nter)
    {
      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) << "
  // _aa_ion_factor_b[last_aa_nter]="s ;
  qDebug() << PeptideFragment::getPeptideIonDirectionName(ion_direction)
           << " last_aa_nter=" << last_aa_nter
           << " _aa_ion_factor_b[last_aa_nter]="
           << _aa_ion_factor_b[last_aa_nter] << " last_aa_cter=" << last_aa_cter
           << " _aa_ion_factor_y[last_aa_cter]="
           << _aa_ion_factor_y[last_aa_cter];
  if(_refine_spectrum_synthesis)
    {
      Pi *= _aa_ion_factor_b[last_aa_nter] * _aa_ion_factor_y[last_aa_cter];
    }
 
  return Pi;
}
 
 
pappso_double
XtandemHyperscore::getHyperscore() const
{
  try
    {
      qDebug() << " _proto_hyperscore=" << _proto_hyperscore;
      double hyperscore = (log10(_proto_hyperscore) * 4);
      if(hyperscore < 0)
        return 0;
      return hyperscore;
    }
  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);
    }
}
 
 
} // namespace pappso