PosteriorErrorProbabilityModel.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: David Wojnar $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/DATASTRUCTURES/DPosition.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/MATH/STATISTICS/GumbelDistributionFitter.h>
#include <OpenMS/MATH/STATISTICS/GumbelMaxLikelihoodFitter.h>
#include <OpenMS/MATH/STATISTICS/GaussFitter.h>
#include <OpenMS/METADATA/PeptideIdentificationList.h>
 
#include <vector>
#include <map>
#include <set>
 
namespace OpenMS
{
  class TextFile;
  class PeptideIdentification;
  class ProteinIdentification;
  class PeptideHit;
  class PeptideIdentificationList;
  namespace Math
  {
 
 
    class OPENMS_DLLAPI PosteriorErrorProbabilityModel :
      public DefaultParamHandler
    {
public:
 
      enum class ModelType
      {
        NONE,                    
        MIXTURE,                 
        TARGET_DECOY_ISOTONIC    
      };
 
      static std::string modelTypeToString(ModelType type);
 
      PosteriorErrorProbabilityModel();
 
      ~PosteriorErrorProbabilityModel() override;
 
      static std::map<std::string, std::vector<std::vector<double>>> extractAndTransformScores(
        const std::vector<ProteinIdentification> & protein_ids,
        const PeptideIdentificationList & peptide_ids,
        const bool split_charge,
        const bool top_hits_only,
        const bool peptide_level = false);
 
      static std::map<std::string, std::vector<std::vector<double>>> extractAndTransformScores(
        const std::vector<ProteinIdentification> & protein_ids,
        const PeptideIdentificationList & peptide_ids,
        const bool split_charge,
        const bool top_hits_only,
        const bool target_decoy_available,
        const double fdr_for_targets_smaller);
 
      static void updateScores(
        const PosteriorErrorProbabilityModel & PEP_model,
        const std::string & search_engine,
        const Int charge,
        const bool prob_correct,
        const bool split_charge,
        std::vector<ProteinIdentification> & protein_ids,
        PeptideIdentificationList & peptide_ids,
        bool & unable_to_fit_data,
        bool & data_might_not_be_well_fit);
 
      static void updateScores(
        const PosteriorErrorProbabilityModel & PEP_model,
        const std::string & search_engine,
        const Int charge,
        const bool prob_correct,
        const bool split_charge,
        std::vector<ProteinIdentification> & protein_ids,
        PeptideIdentificationList & peptide_ids,
        bool & unable_to_fit_data,
        bool & data_might_not_be_well_fit,
        std::set<std::pair<Size, Size>> & updated_hits);
 
      static std::string getSearchEngineKey(const ProteinIdentification & protein_id);
 
      bool fit(std::vector<double>& search_engine_scores, const std::string& outlier_handling);
 
      bool fitGumbelGauss(std::vector<double>& search_engine_scores, const std::string& outlier_handling);
 
      bool fit(std::vector<double>& search_engine_scores, std::vector<double>& probabilities, const std::string& outlier_handling);
 
      bool fitTargetDecoy(const std::vector<double>& target_scores, const std::vector<double>& decoy_scores, const std::string& outlier_handling);
 
      ModelType getModelType() const
      {
        return model_type_;
      }
 
      double getTargetDecoyIncorrectFraction() const
      {
        return td_incorrect_fraction_;
      }
 
      bool hasCalibrationWarning() const
      {
        return calibration_warning_;
      }
 
      double getTargetDecoyAUC() const
      {
        return td_auc_;
      }
 
      std::vector<std::pair<double, double>> getIsotonicPEPCurve() const;
 
      void fillDensities(const std::vector<double> & x_scores, std::vector<double> & incorrect_density, std::vector<double> & correct_density);
      void fillLogDensities(const std::vector<double> & x_scores, std::vector<double> & incorrect_density, std::vector<double> & correct_density);
      void fillLogDensitiesGumbel(const std::vector<double> & x_scores, std::vector<double> & incorrect_density, std::vector<double> & correct_density);
      double computeLogLikelihood(const std::vector<double> & incorrect_density, const std::vector<double> & correct_density) const;
 
      double computeLLAndIncorrectPosteriorsFromLogDensities(
          const std::vector<double>& incorrect_log_density,
          const std::vector<double>& correct_log_density,
          std::vector<double>& incorrect_posterior) const;
 
      std::pair<double, double> pos_neg_mean_weighted_posteriors(const std::vector<double> &x_scores,
                                                                 const std::vector<double> &incorrect_posteriors);
 
      std::pair<double, double> pos_neg_sigma_weighted_posteriors(const std::vector<double> &x_scores,
                                                                 const std::vector<double> &incorrect_posteriors,
                                                                 const std::pair<double, double>& pos_neg_mean);
 
      GaussFitter::GaussFitResult getCorrectlyAssignedFitResult() const
      {
        return correctly_assigned_fit_param_;
      }
 
      GaussFitter::GaussFitResult getIncorrectlyAssignedFitResult() const
      {
        return incorrectly_assigned_fit_param_;
      }
 
      GumbelMaxLikelihoodFitter::GumbelDistributionFitResult getIncorrectlyAssignedGumbelFitResult() const
      {
        return incorrectly_assigned_fit_gumbel_param_;
      }
 
      double getNegativePrior() const
      {
        return negative_prior_;
      }
 
      static double getGumbel_(double x, const GaussFitter::GaussFitResult & params)
      {
        double z = exp((params.x0 - x) / params.sigma);
        return (z * exp(-1 * z)) / params.sigma;
      }
 
      double computeProbability(double score) const;
 
      double computeDecoyProbability(double score) const;
 
      TextFile initPlots(std::vector<double> & x_scores);
 
      const std::string getGumbelGnuplotFormula(const GaussFitter::GaussFitResult & params) const;
 
      const std::string getGaussGnuplotFormula(const GaussFitter::GaussFitResult & params) const;
 
      const std::string getBothGnuplotFormula(const GaussFitter::GaussFitResult & incorrect, const GaussFitter::GaussFitResult & correct) const;
 
      void plotTargetDecoyEstimation(std::vector<double> & target, std::vector<double> & decoy);
 
      inline double getSmallestScore() const
      {
        return smallest_score_;
      }
 
      void tryGnuplot(const std::string& gp_file);
 
private:
      void resetFit_();
 
      void processOutliers_(std::vector<double>& x_scores, const std::string& outlier_handling) const;
 
      bool fitMixture_(std::vector<double> scores, bool use_gumbel, const std::string& outlier_handling);
 
      bool fitTargetDecoyModel_(const std::vector<double>& targets, const std::vector<double>& decoys, const std::string& outlier_handling);
 
      bool fitIsotonic_(const std::vector<double>& target_scores, const std::vector<double>& decoy_scores);
 
      double logIncorrectDensity_(double x) const;
 
      double logCorrectDensity_(double x) const;
 
      double incorrectMode_() const;
 
      static double transformScore_(const std::string& engine, const PeptideHit& hit, const std::string& current_score_type, bool higher_score_better, const std::string& score_source);
 
      static void updateScores_(
        const PosteriorErrorProbabilityModel & PEP_model,
        const std::string & search_engine,
        const Int charge,
        const bool prob_correct,
        const bool split_charge,
        std::vector<ProteinIdentification> & protein_ids,
        PeptideIdentificationList & peptide_ids,
        bool & unable_to_fit_data,
        bool & data_might_not_be_well_fit,
        bool unusable_to_one,
        std::set<std::pair<Size, Size>> & updated_hits);
 
      static std::map<std::string, std::string> chooseScoreSources_(const std::map<std::string, std::string>& run_to_engine, const PeptideIdentificationList& peptide_ids);
 
      static bool getScore_(const std::vector<std::string>& requested_score_types, const PeptideHit & hit, const std::string& actual_score_type, double& value);
 
      PosteriorErrorProbabilityModel & operator=(const PosteriorErrorProbabilityModel & rhs);
      PosteriorErrorProbabilityModel(const PosteriorErrorProbabilityModel & rhs);
      GaussFitter::GaussFitResult incorrectly_assigned_fit_param_;
      GumbelMaxLikelihoodFitter::GumbelDistributionFitResult incorrectly_assigned_fit_gumbel_param_;
      GaussFitter::GaussFitResult correctly_assigned_fit_param_;
      double negative_prior_;
      double smallest_score_;
      const std::string (PosteriorErrorProbabilityModel::* getNegativeGnuplotFormula_)(const GaussFitter::GaussFitResult & params) const;
      const std::string (PosteriorErrorProbabilityModel::* getPositiveGnuplotFormula_)(const GaussFitter::GaussFitResult & params) const;
      ModelType model_type_;
      bool incorrect_is_gumbel_;
      std::vector<double> isotonic_scores_;
      std::vector<double> isotonic_peps_;
      std::vector<double> decoy_scores_;
      std::vector<double> decoy_peps_;
      double td_incorrect_fraction_;
      double td_auc_;
      bool calibration_warning_;
    };
  }
}