computeqvalues.cpp

/**
 * \file core/qvalue/computeqvalues.cpp
 * \date 04/09/2019
 * \author Olivier Langella
 * \brief compute q-value for each peptide evidence (PSM) of a project
 */

/*******************************************************************************
 * Copyright (c) 2019 Olivier Langella <olivier.langella@u-psud.fr>.
 *
 * This file is part of XTPcpp.
 *
 *     XTPcpp 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.
 *
 *     XTPcpp 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 XTPcpp.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

#include <set>
#include "computeqvalues.h"
#include "../../utils/utils.h"
#include <pappsomspp/exception/exceptionnotimplemented.h>

ComputeQvalues::ComputeQvalues(ProjectSp project_sp)
{
  std::vector<IdentificationDataSourceSp> identification_source_list =
    project_sp.get()
      ->getIdentificationDataSourceStore()
      .getIdentificationDataSourceList();


  std::vector<PeptideEvidence *> tandem_peptide_evidence_list;
  for(auto &identification_source_sp : identification_source_list)
    {
      if(identification_source_sp.get()->getIdentificationEngine() !=
         IdentificationEngine::XTandem)
        {
          throw pappso::ExceptionNotImplemented(
            QObject::tr("q-value computation for %1 identification engine is "
                        "not yet implemented")
              .arg(
                identification_source_sp.get()->getIdentificationEngineName()));
        }
      else
        {
          tandem_peptide_evidence_list.reserve(
            tandem_peptide_evidence_list.size() +
            identification_source_sp.get()->getPeptideEvidenceStore().size());
          for(auto &pe_sp : identification_source_sp.get()
                              ->getPeptideEvidenceStore()
                              .getPeptideEvidenceList())
            {
              tandem_peptide_evidence_list.push_back(pe_sp.get());
            }
        }
    }

  for(auto &identification_group :
      project_sp.get()->getIdentificationGroupList())
    {
      for(auto &p_protein_match : identification_group->getProteinMatchList())
        {
          if(p_protein_match->getProteinXtpSp().get()->isDecoy())
            {
              p_protein_match->collectPeptideEvidences(
                m_falsePeptideEvidenceList, ValidationState::notValid);
            }
        }
    }
  computeTandemPeptideEvidenceQvalues(tandem_peptide_evidence_list);
}

ComputeQvalues::~ComputeQvalues()
{
}

void
ComputeQvalues::computeTandemPeptideEvidenceQvalues(
  std::vector<PeptideEvidence *> &tandem_peptide_evidence_list) const
{
  std::sort(tandem_peptide_evidence_list.begin(),
            tandem_peptide_evidence_list.end(),
            [](const PeptideEvidence *pepa, const PeptideEvidence *pepb) {
              return pepa->getEvalue() < pepb->getEvalue();
            });
  std::size_t count_decoy  = 0;
  std::size_t count_target = 0;
  std::vector<PeptideEvidence *> buffer_pep;
  for(PeptideEvidence *pep : tandem_peptide_evidence_list)
    {
      if(buffer_pep.size() > 0)
        {
          if(pep->getEvalue() != buffer_pep.back()->getEvalue())
            { // flush buffer and set q-value
              double qvalue = Utils::computeFdr(count_decoy, count_target);
              for(auto &pepb : buffer_pep)
                {
                  pepb->setParam(PeptideEvidenceParam::pappso_qvalue,
                                 QVariant(qvalue));
                }
              buffer_pep.clear();
            }
        }
      if(m_falsePeptideEvidenceList.find(pep) !=
         m_falsePeptideEvidenceList.end())
        {
          count_decoy++;
        }
      else
        {
          count_target++;
        }

      buffer_pep.push_back(pep);
    }

    //don't forget the last one
  if(buffer_pep.size() > 0)
    {
      double qvalue = Utils::computeFdr(count_decoy, count_target);
      for(auto &pepb : buffer_pep)
        {
          pepb->setParam(PeptideEvidenceParam::pappso_qvalue, QVariant(qvalue));
        }
      buffer_pep.clear();
    }
  // check reverse list to clean q-values
  auto rbegin       = tandem_peptide_evidence_list.rbegin();
  double qvalue_max = 999999999;
  while(rbegin != tandem_peptide_evidence_list.rend())
    {
      if((*rbegin)->getParam(PeptideEvidenceParam::pappso_qvalue).toDouble() >
         qvalue_max)
        {
          (*rbegin)->setParam(PeptideEvidenceParam::pappso_qvalue,
                              QVariant(qvalue_max));
        }
      else
        {
          qvalue_max =
            (*rbegin)->getParam(PeptideEvidenceParam::pappso_qvalue).toDouble();
        }
      rbegin++;
    }
}