PosteriorErrorProbabilityModel Class Reference

Estimates posterior error probabilities (PEPs) of peptide-spectrum matches from search engine scores. More...

#include <OpenMS/MATH/STATISTICS/PosteriorErrorProbabilityModel.h>

Inheritance diagram for PosteriorErrorProbabilityModel:

Public Types
enum class	ModelType { NONE , MIXTURE , TARGET_DECOY_ISOTONIC }
	The kind of model that was estimated by the last call to one of the fit functions. More...

Public Member Functions
	PosteriorErrorProbabilityModel ()
	default constructor
	~PosteriorErrorProbabilityModel () override
	Destructor.
bool	fit (std::vector< double > &search_engine_scores, const std::string &outlier_handling)
	fits the distributions to the data points(search_engine_scores). Estimated parameters for the distributions are saved in member variables. computeProbability can be used afterwards.
bool	fitGumbelGauss (std::vector< double > &search_engine_scores, const std::string &outlier_handling)
	fits the distributions to the data points(search_engine_scores). Estimated parameters for the distributions are saved in member variables. computeProbability can be used afterwards. Uses Gumbel+Gauss (independent of the parameter incorrectly_assigned). Fits the Gumbel distribution by maximizing the (weighted) likelihood.
bool	fit (std::vector< double > &search_engine_scores, std::vector< double > &probabilities, const std::string &outlier_handling)
	fits the distributions to the data points(search_engine_scores) and writes the computed probabilities into the given vector (the second one).
bool	fitTargetDecoy (const std::vector< double > &target_scores, const std::vector< double > &decoy_scores, const std::string &outlier_handling)
	estimates the PEP model from scores of target and decoy PSMs. computeProbability() (targets) and computeDecoyProbability() (decoys) can be used afterwards.
ModelType	getModelType () const
	returns the kind of model estimated by the last fit
double	getTargetDecoyIncorrectFraction () const
	returns the target-decoy estimate of the fraction of incorrect target PSMs, i.e. (number of decoys) / (number of targets) / decoy_ratio (capped at 1), as computed by the last call to fitTargetDecoy(). Returns -1 if it is not available (also if the last fit failed).
bool	hasCalibrationWarning () const
	returns true if the last call to fitTargetDecoy() found the target-decoy assumptions to be violated at low scores
double	getTargetDecoyAUC () const
	returns the area under the ROC curve for separating target from decoy PSMs by score (ties count half), as computed by the last call to fitTargetDecoy(), or -1 if it is not available (no decoys, decoy_handling = ignore, or the last fit failed).
std::vector< std::pair< double, double > >	getIsotonicPEPCurve () const
	returns the PEP curve of the non-parametric target-decoy estimate as pairs (score, PEP), sorted by score.
void	fillDensities (const std::vector< double > &x_scores, std::vector< double > &incorrect_density, std::vector< double > &correct_density)
	Writes the distributions densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.
void	fillLogDensities (const std::vector< double > &x_scores, std::vector< double > &incorrect_density, std::vector< double > &correct_density)
	Writes the log distributions densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.
void	fillLogDensitiesGumbel (const std::vector< double > &x_scores, std::vector< double > &incorrect_density, std::vector< double > &correct_density)
	Writes the log distributions of gumbel and gauss densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.
double	computeLogLikelihood (const std::vector< double > &incorrect_density, const std::vector< double > &correct_density) const
	computes the Likelihood with a log-likelihood function.
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
	returns estimated parameters for correctly assigned sequences. Fit should be used before.
GaussFitter::GaussFitResult	getIncorrectlyAssignedFitResult () const
	returns estimated parameters for incorrectly assigned sequences (for the Gumbel distribution, x0 and sigma are the location and scale parameter). Fit should be used before.
GumbelMaxLikelihoodFitter::GumbelDistributionFitResult	getIncorrectlyAssignedGumbelFitResult () const
	returns estimated Gumbel parameters for incorrectly assigned sequences. Fit should be used before.
double	getNegativePrior () const
	returns the estimated negative prior probability.
double	computeProbability (double score) const
double	computeDecoyProbability (double score) const
	returns the posterior error probability of a decoy PSM with the given score
TextFile	initPlots (std::vector< double > &x_scores)
	initializes the plots
const std::string	getGumbelGnuplotFormula (const GaussFitter::GaussFitResult &params) const
	returns the gnuplot formula of the fitted gumbel distribution. Only x0 and sigma are used as local parameter alpha and scale parameter beta, respectively.
const std::string	getGaussGnuplotFormula (const GaussFitter::GaussFitResult &params) const
	returns the gnuplot formula of the fitted gauss distribution.
const std::string	getBothGnuplotFormula (const GaussFitter::GaussFitResult &incorrect, const GaussFitter::GaussFitResult &correct) const
	returns the gnuplot formula of the fitted mixture distribution.
void	plotTargetDecoyEstimation (std::vector< double > &target, std::vector< double > &decoy)
	plots the score distributions of target and decoy PSMs together with the estimated model
double	getSmallestScore () const
	returns the smallest score used in the last fit
void	tryGnuplot (const std::string &gp_file)
	try to invoke 'gnuplot' on the file to create PDF automatically
Public Member Functions inherited from DefaultParamHandler
	DefaultParamHandler (const std::string &name)
	Constructor with name that is displayed in error messages.
	DefaultParamHandler (const DefaultParamHandler &rhs)
	Copy constructor.
virtual	~DefaultParamHandler ()
	Destructor.
DefaultParamHandler &	operator= (const DefaultParamHandler &rhs)
	Assignment operator.
virtual bool	operator== (const DefaultParamHandler &rhs) const
	Equality operator.
void	setParameters (const Param &param)
	Sets the parameters.
const Param &	getParameters () const
	Non-mutable access to the parameters.
const Param &	getDefaults () const
	Non-mutable access to the default parameters.
const std::string &	getName () const
	Non-mutable access to the name.
void	setName (const std::string &name)
	Mutable access to the name.
const std::vector< std::string > &	getSubsections () const
	Non-mutable access to the registered subsections.

Static Public Member Functions
static std::string	modelTypeToString (ModelType type)
	human-readable name of a ModelType
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)
	extract and transform score types to a range and score orientation that the PEP model can handle
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)
	extract and transform score types to a range and score orientation that the PEP model can handle
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)
	update score entries with PEP (or 1-PEP) estimates
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)
	update score entries with PEP (or 1-PEP) estimates and report which hits were updated
static std::string	getSearchEngineKey (const ProteinIdentification &protein_id)
	the key under which the scores of a search run are reported by extractAndTransformScores()
static double	getGumbel_ (double x, const GaussFitter::GaussFitResult &params)
	computes the gumbel density at position x with parameters params.
Static Public Member Functions inherited from DefaultParamHandler
static void	writeParametersToMetaValues (const Param &write_this, MetaInfoInterface &write_here, const std::string &key_prefix="")
	Writes all parameters to meta values.

Private Member Functions
void	resetFit_ ()
	resets the state of the last fit (model type, target-decoy diagnostics, non-parametric curve)
void	processOutliers_ (std::vector< double > &x_scores, const std::string &outlier_handling) const
	remove (or clamp) outliers for fitting
bool	fitMixture_ (std::vector< double > scores, bool use_gumbel, const std::string &outlier_handling)
	fits a two-component mixture model with the EM algorithm
bool	fitTargetDecoyModel_ (const std::vector< double > &targets, const std::vector< double > &decoys, const std::string &outlier_handling)
	fits the target-decoy model (or, without decoys or with decoy_handling 'ignore', the mixture model) to finite target and decoy scores
bool	fitIsotonic_ (const std::vector< double > &target_scores, const std::vector< double > &decoy_scores)
	estimates the non-parametric target-decoy model; returns false if it cannot be estimated
double	logIncorrectDensity_ (double x) const
	log density of the incorrect component at x (in the shifted score space)
double	logCorrectDensity_ (double x) const
	log density of the correct component at x (in the shifted score space)
double	incorrectMode_ () const
	mode of the incorrect component (in the shifted score space)
PosteriorErrorProbabilityModel &	operator= (const PosteriorErrorProbabilityModel &rhs)
	assignment operator (not implemented)
	PosteriorErrorProbabilityModel (const PosteriorErrorProbabilityModel &rhs)
	Copy constructor (not implemented)

Static Private Member Functions
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)
	implementation of updateScores(); hits without a usable score get a PEP of 1 if unusable_to_one is set, otherwise they are left unchanged
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)

Private Attributes
GaussFitter::GaussFitResult	incorrectly_assigned_fit_param_
	stores parameters for incorrectly assigned sequences. If gumbel fit was used, A can be ignored. Furthermore, in this case, x0 and sigma are the local parameter alpha and scale parameter beta, respectively.
GumbelMaxLikelihoodFitter::GumbelDistributionFitResult	incorrectly_assigned_fit_gumbel_param_
GaussFitter::GaussFitResult	correctly_assigned_fit_param_
	stores gauss parameters
double	negative_prior_
	stores final prior probability for negative peptides
double	smallest_score_
	smallest score which was used for fitting the model
const std::string(PosteriorErrorProbabilityModel::*	getNegativeGnuplotFormula_ )(const GaussFitter::GaussFitResult &params) const
	points either to getGumbelGnuplotFormula or getGaussGnuplotFormula depending on whether one uses the gumbel or the gaussian distribution for incorrectly assigned sequences.
const std::string(PosteriorErrorProbabilityModel::*	getPositiveGnuplotFormula_ )(const GaussFitter::GaussFitResult &params) const
	points to getGumbelGnuplotFormula
ModelType	model_type_
	kind of the last estimated model
bool	incorrect_is_gumbel_
	whether the incorrect component of the (last) mixture model is a Gumbel distribution (otherwise Gauss)
std::vector< double >	isotonic_scores_
	knots of the non-parametric PEP curve: scores (ascending) and PEPs (non-increasing)
std::vector< double >	isotonic_peps_
std::vector< double >	decoy_scores_
	knots of the PEP curve of decoys (label-swapped estimates at the scores of the fitted decoys): scores (ascending) and PEPs (non-increasing)
std::vector< double >	decoy_peps_
double	td_incorrect_fraction_
	target-decoy estimate of the fraction of incorrect targets (-1 if not available)
double	td_auc_
	area under the ROC curve of targets vs. decoys of the last target-decoy fit (-1 if not available)
bool	calibration_warning_
	whether the last target-decoy fit found too few decoys among the lowest-scoring PSMs (see hasCalibrationWarning())

Additional Inherited Members
Protected Member Functions inherited from DefaultParamHandler
virtual void	updateMembers_ ()
	This method is used to update extra member variables at the end of the setParameters() method.
void	defaultsToParam_ ()
	Updates the parameters after the defaults have been set in the constructor.
Protected Attributes inherited from DefaultParamHandler
Param	param_
	Container for current parameters.
Param	defaults_
	Container for default parameters. This member should be filled in the constructor of derived classes!
std::vector< std::string >	subsections_
	Container for registered subsections. This member should be filled in the constructor of derived classes!
std::string	error_name_
	Name that is displayed in error messages during the parameter checking.
bool	check_defaults_
	If this member is set to false no checking if parameters in done;.
bool	warn_empty_defaults_
	If this member is set to false no warning is emitted when defaults are empty;.

Detailed Description

Estimates posterior error probabilities (PEPs) of peptide-spectrum matches from search engine scores.

Two estimators are available:

• Unsupervised mixture model (fit(), and fitTargetDecoy() without decoys or with decoy_handling = ignore, where targets and decoys are pooled): a two-component mixture of a Gumbel (or Gauss) distribution for incorrect and a Gauss distribution for correct assignments is fitted to all scores with the EM algorithm. The very same densities are used for fitting and for computing probabilities. The Gumbel component is fitted by weighted maximum likelihood.
• Non-parametric target-decoy estimate (fitTargetDecoy() with at least one decoy and decoy_handling = auto, the default): Under target-decoy competition (concatenated target+decoy search, top hit per spectrum), each decoy winner stands in for (about) one incorrect target winner. The probability p(s) that a PSM with score s is a decoy is estimated by weighted isotonic regression (pool-adjacent-violators, tied scores pooled) and the PEP of a target PSM is PEP(s) = p(s) / (1 - p(s)) / decoy_ratio (cf. Percolator/qvality, Käll et al. 2008). The blocks and their PEPs only depend on the ordering of the scores, so they are invariant to monotone score transformations (e-values, p-values, raw scores, discrete or censored scores). Block decoy fractions are shrunk with a Jeffreys prior, and the PEPs of the targets scoring above all decoys sum up to 1 / decoy_ratio (the expected number of incorrect targets above the best decoy), so that no PEP is exactly 0. Even a single decoy (or a single target) is informative. Within a block, the PEP decreases linearly with the score around the block's PEP (keeping the mean PEP of the block), such that PSMs are ordered by score; between blocks, the PEP is interpolated linearly. These ramps use the score values, i.e. they (but not the block PEPs) depend on the score scale. As in qvality, Percolator and mokapot, the PEP is not forced to 1 at the lowest scores, where it follows the local ratio of decoys and targets (see hasCalibrationWarning()). Decoys get their PEPs from a symmetric (label-swapped) estimate, see computeDecoyProbability(): the curve is fitted with the decoy labels, which makes it more favorable for targets than for decoys at the same score. Without the label swap, an incorrect target can score above the best decoy and get a low PEP, but a decoy never can. Downstream target-decoy FDR estimates on the PEPs (e.g. protein FDR after protein inference across many small runs) would then miss these incorrect targets.

Note

Without decoys, the (unsupervised) mixture model assumes that the data contain correct assignments. For data with few or no correct assignments, it might assign low PEPs to a few PSMs in the extreme upper tail of the scores. Its parametric densities are also only an approximation of real score distributions (on Comet e-values of a yeast/UPS1 dataset, its PEPs were about 1.5 times too high in total). Use target/decoy information whenever possible. With fewer than min_mixture_scores scores (default 35, as PeptideProphet; also counted after outlier handling), the mixture model is not fitted at all.

One can output the fit as a gnuplot formula using getGumbelGnuplotFormula() and getGaussGnuplotFormula() after fitting.

Note

All parameters are stored in GaussFitResult. In the case of the Gumbel distribution x0 and sigma represent the local parameter alpha and the scale parameter beta, respectively.

Todo:

allow non-parametric via kernel density estimation

Parameters of this class are:

Name	Type	Default	Restrictions	Description
out_plot	string			If given, the some output files will be saved in the following manner: _scores.txt for the scores and which contains the fitted values for each step of the EM-algorithm, e.g., out_plot = /usr/home/OMSSA123 leads to /usr/home/OMSSA123_scores.txt, /usr/home/OMSSA123 will be written. If no directory is specified, e.g. instead of '/usr/home/OMSSA123' just OMSSA123, the files will be written into the working directory.
number_of_bins	int	100		Number of bins used for visualization. Only needed if each iteration step of the EM-Algorithm will be visualized
incorrectly_assigned	string	Gumbel	Gumbel, Gauss	Distribution of the scores of incorrectly assigned sequences in the mixture model. For 'Gumbel', the Gumbel distribution is used (fitted by maximum likelihood). For 'Gauss', the Gauss distribution is used.
max_nr_iterations	int	1000		Bounds the number of iterations for the EM algorithm when convergence is slow.
neg_log_delta	int	6		The negative logarithm of the convergence threshold for the likelihood increase.
outlier_handling	string	ignore_iqr_outliers	ignore_iqr_outliers, set_iqr_to_closest_valid, ignore_extreme_percentiles, none	What to do with outliers (mixture model only): - ignore_iqr_outliers: ignore outliers outside of 3*IQR from Q1/Q3 for fitting - set_iqr_to_closest_valid: set IQR-based outliers to the last valid value for fitting - ignore_extreme_percentiles: ignore everything outside 99.9th and 1st percentile (also removes equal values like potential censored max values in XTandem) - none: do nothing
decoy_handling	string	auto	auto, ignore	How target/decoy annotations are used: - auto: if decoy PSMs are present (even a single one), non-parametric target-decoy estimate from the local ratio of decoys and targets (isotonic regression); its blocks depend only on the order of the scores, i.e. it is suitable for any score type. Without decoys, mixture model of the target scores. - ignore: ignore the annotations, i.e. targets and decoys are pooled and modeled without supervision (behavior of OpenMS 3.5 and earlier) The target-decoy estimate assumes a search against a concatenated target+decoy database (target-decoy competition). Only top hits are used for fitting.
decoy_ratio	float	1.0	min: 1.0e-06	Number of decoy sequences per target sequence in the searched database (cf. DecoyDatabase's 'shuffle_decoy_ratio'). Under target-decoy competition, the number of incorrect targets is estimated as the number of decoys divided by this ratio.
min_decoys	int	10	min: 1	Only used by IDPosteriorErrorProbability with 'split_charge': charge states with fewer decoy PSMs use the target-decoy data of all charge states (1: separate models for all charge states with decoys).
min_mixture_scores	int	35	min: 2	Minimum number of (finite) scores to fit the mixture model (used without decoys, or with decoy_handling 'ignore'), before and after outlier handling. With fewer scores, no model is fitted: the scores of so few PSMs without decoys do not determine the two components (PeptideProphet also requires 35 PSMs per model). IDPosteriorErrorProbability with 'split_charge' then uses the mixture model of all charge states.

Note:

• If a section name is documented, the documentation is displayed as tooltip.
• Advanced parameter names are italic.

Member Enumeration Documentation

ModelType

enum class ModelType

strong

The kind of model that was estimated by the last call to one of the fit functions.

Enumerator
NONE	nothing was fitted (yet), or the last fit failed
MIXTURE	unsupervised two-component mixture model
TARGET_DECOY_ISOTONIC	non-parametric target-decoy estimate (isotonic regression of the decoy probability)

Constructor & Destructor Documentation

PosteriorErrorProbabilityModel() [1/2]

PosteriorErrorProbabilityModel

(

)

default constructor

~PosteriorErrorProbabilityModel()

~PosteriorErrorProbabilityModel ( )

override

Destructor.

PosteriorErrorProbabilityModel() [2/2]

PosteriorErrorProbabilityModel ( const PosteriorErrorProbabilityModel &  rhs )

private

Copy constructor (not implemented)

Member Function Documentation

chooseScoreSources_()

static std::map< std::string, std::string > chooseScoreSources_ ( const std::map< std::string, std::string > &  run_to_engine, const PeptideIdentificationList &  peptide_ids  )

staticprivate

For each engine (key) of run_to_engine (run identifier -> engine key): the score type to use for all of its hits. For engines with dedicated support, the most preferred of their score types that any hit has; otherwise (or if no hit has any of them) an empty string, i.e. the main score of the identifications. One score type per engine keeps the scores of a model on the same scale.

computeDecoyProbability()

double computeDecoyProbability

(

double

score

)

const

returns the posterior error probability of a decoy PSM with the given score

For the non-parametric target-decoy estimate, the PEP of a decoy that was used for fitting is the PEP it would get as a target, i.e. from the estimate in which this decoy is labeled as a target (label swap). This makes decoy and incorrect target PSMs exchangeable (the same distribution of PEPs), which target-decoy FDR estimates based on the PEPs require. Decoy PEPs are therefore lower than or equal to the PEPs of targets with the same score. Exact for the best scoring decoys (all decoys of smaller datasets); for decoys in large blocks of large datasets, where the effect is small, the change of the block PEP is used. Scores in between the scores of the fitted decoys are interpolated linearly. For the mixture model, the PEP equals computeProbability().

computeLLAndIncorrectPosteriorsFromLogDensities()

double computeLLAndIncorrectPosteriorsFromLogDensities	(	const std::vector< double > &	incorrect_log_density,
		const std::vector< double > &	correct_log_density,
		std::vector< double > &	incorrect_posterior
	)		const

computes the posteriors for the datapoints to belong to the incorrect distribution

Parameters

[in]	incorrect_log_density	...
[in]	correct_log_density	...
[out]	incorrect_posterior	resulting posteriors; same size as incorrect_log_density

Returns

the log-likelihood of the model

computeLogLikelihood()

double computeLogLikelihood	(	const std::vector< double > &	incorrect_density,
		const std::vector< double > &	correct_density
	)		const

computes the Likelihood with a log-likelihood function.

computeProbability()

double computeProbability

(

double

score

)

const

Returns the computed posterior error probability for a given score.

The probability is non-increasing with the score.

Note

: fit has to be used before using this function. Otherwise this function will compute nonsense.

extractAndTransformScores() [1/2]

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

extract and transform score types to a range and score orientation that the PEP model can handle

Scores are transformed such that higher values are better (e.g. -log10 of e-values).

Parameters

[in]	protein_ids	the protein identifications
[in]	peptide_ids	the peptide identifications
[in]	split_charge	whether different charge states should be treated separately
[in]	top_hits_only	only consider rank 1 (affects the first vector of the result only)
[in]	peptide_level	report one score per peptide (modified sequence; with split_charge: per peptide and charge state) instead of one per hit: in [0] the best score of the considered hits of the peptide, in [1] and [2] the best top hit of the peptide (after the competition per spectrum). Hits without a sequence are counted individually.

Returns

engine (and optional charge state) id -> vector of three score vectors: [0] scores of all considered hits (decoys included, no labels used), [1] scores of the top hits of all spectra that are annotated as target ("target" or "target+decoy"), [2] scores of the top hits of all spectra that are annotated as decoy. [1] and [2] are empty if no hit is annotated as decoy. Groups with very few scores are included as well (unlike in the deprecated overload), although they might be too small to fit a model.

Note

Supported engines (i.e. with a dedicated score selection) are: XTandem, OMSSA, Mascot, SpectraST, MyriMatch, SimTandem, MSGFPlus/MS-GF+, Comet, MSFragger, Sage, SimpleSearchEngine and ProSE. For other engines, the main score of the identifications is used, if it is not on a probability scale (e.g. q-values or PEPs).

extractAndTransformScores() [2/2]

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

extract and transform score types to a range and score orientation that the PEP model can handle

Deprecated:

Use the overload without target_decoy_available and fdr_for_targets_smaller. Target and decoy scores are now determined from the target/decoy annotation of the hits, fdr_for_targets_smaller is ignored.

Parameters

[in]	protein_ids	the protein identifications
[in]	peptide_ids	the peptide identifications
[in]	split_charge	whether different charge states should be treated separately
[in]	top_hits_only	only consider rank 1
[in]	target_decoy_available	if false, the target and decoy score vectors are left empty
[in]	fdr_for_targets_smaller	ignored

Returns

engine (and optional charge state) id -> vector of triplets (score, target, decoy); only groups with more than two scores

fillDensities()

void fillDensities	(	const std::vector< double > &	x_scores,
		std::vector< double > &	incorrect_density,
		std::vector< double > &	correct_density
	)

Writes the distributions densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.

fillLogDensities()

void fillLogDensities	(	const std::vector< double > &	x_scores,
		std::vector< double > &	incorrect_density,
		std::vector< double > &	correct_density
	)

Writes the log distributions densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.

fillLogDensitiesGumbel()

void fillLogDensitiesGumbel	(	const std::vector< double > &	x_scores,
		std::vector< double > &	incorrect_density,
		std::vector< double > &	correct_density
	)

Writes the log distributions of gumbel and gauss densities into the two vectors for a set of scores. Incorrect_densities represent the incorrectly assigned sequences.

fit() [1/2]

bool fit	(	std::vector< double > &	search_engine_scores,
		const std::string &	outlier_handling
	)

fits the distributions to the data points(search_engine_scores). Estimated parameters for the distributions are saved in member variables. computeProbability can be used afterwards.

Fits a Gumbel (parameter incorrectly_assigned = Gumbel) or a Gauss distribution for the incorrectly assigned and a Gauss distribution for the correctly assigned sequences with the EM algorithm. The same distributions are used to calculate the final probabilities.

Parameters

[in,out]	search_engine_scores	a vector which holds the data points
[in]	outlier_handling	Valid values are in the Param 'outlier_handling'

Returns

true if algorithm has run through. Else false will be returned. In that case no plot and no probabilities are calculated.

Note

the vector is sorted from smallest to biggest value!

non-finite scores are ignored for fitting

with fewer finite scores than the parameter min_mixture_scores (default 35), before or after outlier handling, no model is fitted (returns false)

fit() [2/2]

bool fit	(	std::vector< double > &	search_engine_scores,
		std::vector< double > &	probabilities,
		const std::string &	outlier_handling
	)

fits the distributions to the data points(search_engine_scores) and writes the computed probabilities into the given vector (the second one).

Parameters

[in,out]	search_engine_scores	a vector which holds the data points
[out]	probabilities	Probability for each data point after running this function. If it has some content it will be overwritten.
[in]	outlier_handling	Valid values are in the Param 'outlier_handling'

Returns

true if algorithm has run through. Else false will be returned. In that case no plot and no probabilities are calculated.

Note

the vectors are sorted from smallest to biggest value!

fitGumbelGauss()

bool fitGumbelGauss	(	std::vector< double > &	search_engine_scores,
		const std::string &	outlier_handling
	)

fits the distributions to the data points(search_engine_scores). Estimated parameters for the distributions are saved in member variables. computeProbability can be used afterwards. Uses Gumbel+Gauss (independent of the parameter incorrectly_assigned). Fits the Gumbel distribution by maximizing the (weighted) likelihood.

Parameters

[in,out]	search_engine_scores	a vector which holds the data points
[in]	outlier_handling	Valid values are in the Param 'outlier_handling'

Returns

true if algorithm has run through. Else false will be returned. In that case no plot and no probabilities are calculated.

Note

the vector is sorted from smallest to biggest value!

fitIsotonic_()

bool fitIsotonic_ ( const std::vector< double > &  target_scores, const std::vector< double > &  decoy_scores  )

private

estimates the non-parametric target-decoy model; returns false if it cannot be estimated

fitMixture_()

bool fitMixture_ ( std::vector< double >  scores, bool  use_gumbel, const std::string &  outlier_handling  )

private

fits a two-component mixture model with the EM algorithm

Parameters

[in]	scores	(finite) scores of the PSMs
[in]	use_gumbel	use a Gumbel (instead of a Gauss) distribution for the incorrect component
[in]	outlier_handling	see parameter 'outlier_handling'

Returns

true on success

fitTargetDecoy()

bool fitTargetDecoy	(	const std::vector< double > &	target_scores,
		const std::vector< double > &	decoy_scores,
		const std::string &	outlier_handling
	)

estimates the PEP model from scores of target and decoy PSMs. computeProbability() (targets) and computeDecoyProbability() (decoys) can be used afterwards.

The scores must be oriented such that higher scores are better (as returned by extractAndTransformScores()). They should stem from a target-decoy competition, i.e. a search against a concatenated target+decoy database, with one (the top) hit per spectrum. With at least one decoy, the non-parametric target-decoy estimate is used; without decoys (or with decoy_handling = ignore), the unsupervised mixture model (see class documentation). Use getModelType() to find out which estimator was used.

Parameters

[in]	target_scores	scores of target PSMs
[in]	decoy_scores	scores of decoy PSMs
[in]	outlier_handling	Valid values are in the Param 'outlier_handling' (only used for the mixture model)

Returns

true if a model could be estimated

fitTargetDecoyModel_()

bool fitTargetDecoyModel_ ( const std::vector< double > &  targets, const std::vector< double > &  decoys, const std::string &  outlier_handling  )

private

fits the target-decoy model (or, without decoys or with decoy_handling 'ignore', the mixture model) to finite target and decoy scores

getBothGnuplotFormula()

const std::string getBothGnuplotFormula	(	const GaussFitter::GaussFitResult &	incorrect,
		const GaussFitter::GaussFitResult &	correct
	)		const

returns the gnuplot formula of the fitted mixture distribution.

getCorrectlyAssignedFitResult()

GaussFitter::GaussFitResult getCorrectlyAssignedFitResult ( ) const

inline

returns estimated parameters for correctly assigned sequences. Fit should be used before.

getGaussGnuplotFormula()

const std::string getGaussGnuplotFormula

(

const GaussFitter::GaussFitResult &

params

)

const

returns the gnuplot formula of the fitted gauss distribution.

getGumbel_()

static double getGumbel_ ( double  x, const GaussFitter::GaussFitResult &  params  )

inlinestatic

computes the gumbel density at position x with parameters params.

References GaussFitter::GaussFitResult::sigma, and GaussFitter::GaussFitResult::x0.

getGumbelGnuplotFormula()

const std::string getGumbelGnuplotFormula

(

const GaussFitter::GaussFitResult &

params

)

const

returns the gnuplot formula of the fitted gumbel distribution. Only x0 and sigma are used as local parameter alpha and scale parameter beta, respectively.

getIncorrectlyAssignedFitResult()

GaussFitter::GaussFitResult getIncorrectlyAssignedFitResult ( ) const

inline

returns estimated parameters for incorrectly assigned sequences (for the Gumbel distribution, x0 and sigma are the location and scale parameter). Fit should be used before.

getIncorrectlyAssignedGumbelFitResult()

GumbelMaxLikelihoodFitter::GumbelDistributionFitResult getIncorrectlyAssignedGumbelFitResult ( ) const

inline

returns estimated Gumbel parameters for incorrectly assigned sequences. Fit should be used before.

getIsotonicPEPCurve()

std::vector< std::pair< double, double > > getIsotonicPEPCurve

(

)

const

returns the PEP curve of the non-parametric target-decoy estimate as pairs (score, PEP), sorted by score.

Probabilities of scores in between are interpolated linearly, probabilities of scores outside are constant. Empty unless the last fit was of type ModelType::TARGET_DECOY_ISOTONIC.

getModelType()

ModelType getModelType ( ) const

inline

returns the kind of model estimated by the last fit

getNegativePrior()

double getNegativePrior ( ) const

inline

returns the estimated negative prior probability.

getScore_()

static bool getScore_ ( const std::vector< std::string > &  requested_score_types, const PeptideHit &  hit, const std::string &  actual_score_type, double &  value  )

staticprivate

gets a specific score (either main score [preferred] or metavalue)

Parameters

[in]	requested_score_types	the requested score_types in order of preference (will be tested with a "_score" suffix as well)
[in]	hit	the PeptideHit to extract from
[in]	actual_score_type	the current score type to take preference if matching
[out]	value	the score

Returns

true if one of the requested scores was found

getSearchEngineKey()

static std::string getSearchEngineKey ( const ProteinIdentification &  protein_id )

static

the key under which the scores of a search run are reported by extractAndTransformScores()

This is the upper-case name of the search engine without non-alphanumeric characters, with aliases resolved (e.g. "X! Tandem" -> "XTANDEM", "MS-GF+" -> "MSGFPLUS").

getSmallestScore()

double getSmallestScore ( ) const

inline

returns the smallest score used in the last fit

getTargetDecoyAUC()

double getTargetDecoyAUC ( ) const

inline

returns the area under the ROC curve for separating target from decoy PSMs by score (ties count half), as computed by the last call to fitTargetDecoy(), or -1 if it is not available (no decoys, decoy_handling = ignore, or the last fit failed).

Values of 0.5 or below mean that targets do not score better than decoys (e.g. a wrong score orientation, or a group of PSMs without correct identifications whose targets and decoys did not compete fairly).

getTargetDecoyIncorrectFraction()

double getTargetDecoyIncorrectFraction ( ) const

inline

returns the target-decoy estimate of the fraction of incorrect target PSMs, i.e. (number of decoys) / (number of targets) / decoy_ratio (capped at 1), as computed by the last call to fitTargetDecoy(). Returns -1 if it is not available (also if the last fit failed).

hasCalibrationWarning()

bool hasCalibrationWarning ( ) const

inline

returns true if the last call to fitTargetDecoy() found the target-decoy assumptions to be violated at low scores

A warning is raised if targets do not score better than decoys (see getTargetDecoyAUC()), or if even the lowest-scoring PSMs contain significantly (binomial test) and clearly too few decoys for incorrect PSMs (decoy fraction below decoy_ratio / (1 + decoy_ratio), PEP below 0.8), so that the PEPs stay well below 1 down to the lowest scores. Then either correct PSMs extend to the lowest scores, or decoys are underrepresented among the incorrect PSMs (e.g. a wrong decoy_ratio), which makes all target-decoy PEPs (and FDRs) too low. The set of lowest-scoring PSMs with the highest decoy fraction among those with at least a quarter of the estimated number of incorrect PSMs (number of decoys * (1 + 1 / decoy_ratio) / 4) is tested. The check is skipped if this minimum set size is below 100 PSMs (i.e. with fewer than 200 decoys at decoy_ratio 1), and with decoy_handling = ignore.

incorrectMode_()

double incorrectMode_ ( ) const

private

mode of the incorrect component (in the shifted score space)

initPlots()

TextFile initPlots

(

std::vector< double > &

x_scores

)

initializes the plots

logCorrectDensity_()

double logCorrectDensity_ ( double  x ) const

private

log density of the correct component at x (in the shifted score space)

logIncorrectDensity_()

double logIncorrectDensity_ ( double  x ) const

private

log density of the incorrect component at x (in the shifted score space)

modelTypeToString()

static std::string modelTypeToString ( ModelType  type )

static

human-readable name of a ModelType

operator=()

PosteriorErrorProbabilityModel & operator= ( const PosteriorErrorProbabilityModel &  rhs )

private

assignment operator (not implemented)

plotTargetDecoyEstimation()

void plotTargetDecoyEstimation	(	std::vector< double > &	target,
		std::vector< double > &	decoy
	)

plots the score distributions of target and decoy PSMs together with the estimated model

The decoy histogram is scaled such that, under target-decoy competition, it estimates the score distribution of the incorrect target PSMs.

Parameters

[in,out]	target	scores of target PSMs (will be sorted)
[in,out]	decoy	scores of decoy PSMs (will be sorted)

pos_neg_mean_weighted_posteriors()

std::pair< double, double > pos_neg_mean_weighted_posteriors	(	const std::vector< double > &	x_scores,
		const std::vector< double > &	incorrect_posteriors
	)

Parameters

[in]	x_scores	Scores observed "on the x-axis"
[in]	incorrect_posteriors	Posteriors/responsibilities of belonging to the incorrect component

Returns

New estimate for the mean of the correct (pair.first) and incorrect (pair.second) component

Note

only for Gaussian estimates

pos_neg_sigma_weighted_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
	)

Parameters

[in]	x_scores	Scores observed "on the x-axis"
[in]	incorrect_posteriors	Posteriors/responsibilities of belonging to the incorrect component
[in]	pos_neg_mean	Positive(correct) and negative(incorrect) means, respectively

Returns

New estimate for the std. deviation of the correct (pair.first) and incorrect (pair.second) component

Note

only for Gaussian estimates

processOutliers_()

void processOutliers_ ( std::vector< double > &  x_scores, const std::string &  outlier_handling  ) const

private

remove (or clamp) outliers for fitting

resetFit_()

void resetFit_ ( )

private

resets the state of the last fit (model type, target-decoy diagnostics, non-parametric curve)

transformScore_()

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  )

staticprivate

transform different score types to a range and score orientation that the model can handle (engine string is assumed in upper-case)

Parameters

[in]	engine	the search engine key as returned by getSearchEngineKey()
[in]	hit	the PeptideHit to extract transformed scores from
[in]	current_score_type	the current score type of the PeptideIdentification to take precedence
[in]	higher_score_better	the score orientation of the PeptideIdentification (used for the main score of unknown engines)
[in]	score_source	the score type to use (see chooseScoreSources_()); empty: the main score of the identification

Returns

transformed score (higher is better), or NaN if no usable score is available

tryGnuplot()

void tryGnuplot

(

const std::string &

gp_file

)

try to invoke 'gnuplot' on the file to create PDF automatically

updateScores() [1/2]

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

update score entries with PEP (or 1-PEP) estimates

Parameters

[in]	PEP_model	the PEP model used to update the scores
[in]	search_engine	the score of search_engine will be updated (key as returned by extractAndTransformScores())
[in]	charge	identifications with the given charge will be updated
[in]	prob_correct	report 1-PEP
[in]	split_charge	if charge states have been treated separately
[in,out]	protein_ids	the protein identifications
[in,out]	peptide_ids	the peptide identifications
[out]	unable_to_fit_data	there was a problem fitting the data (probabilities are all smaller 0 or larger 1)
[out]	data_might_not_be_well_fit	fit was successful but of bad quality (probabilities are all smaller 0.8 and larger 0.2)

updateScores() [2/2]

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

update score entries with PEP (or 1-PEP) estimates and report which hits were updated

Same as the overload without updated_hits, except for hits without a usable score: these are left unchanged (and not added to updated_hits), so that the caller can treat them consistently with the other hits of their run (the overload without updated_hits assigns them a PEP of 1).

Parameters

[in]	PEP_model	the PEP model used to update the scores
[in]	search_engine	the score of search_engine will be updated (key as returned by extractAndTransformScores())
[in]	charge	identifications with the given charge will be updated
[in]	prob_correct	report 1-PEP
[in]	split_charge	if charge states have been treated separately
[in,out]	protein_ids	the protein identifications
[in,out]	peptide_ids	the peptide identifications
[out]	unable_to_fit_data	there was a problem fitting the data (probabilities are all smaller 0 or larger 1)
[out]	data_might_not_be_well_fit	fit was successful but of bad quality (probabilities are all smaller 0.8 and larger 0.2)
[in,out]	updated_hits	pairs (index into peptide_ids, index of the hit) of the updated hits are added

updateScores_()

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  )

staticprivate

implementation of updateScores(); hits without a usable score get a PEP of 1 if unusable_to_one is set, otherwise they are left unchanged

Member Data Documentation

calibration_warning_

bool calibration_warning_

private

whether the last target-decoy fit found too few decoys among the lowest-scoring PSMs (see hasCalibrationWarning())

correctly_assigned_fit_param_

GaussFitter::GaussFitResult correctly_assigned_fit_param_

private

stores gauss parameters

decoy_peps_

std::vector<double> decoy_peps_

private

decoy_scores_

std::vector<double> decoy_scores_

private

knots of the PEP curve of decoys (label-swapped estimates at the scores of the fitted decoys): scores (ascending) and PEPs (non-increasing)

getNegativeGnuplotFormula_

const std::string(PosteriorErrorProbabilityModel::* getNegativeGnuplotFormula_) (const GaussFitter::GaussFitResult &params) const

private

points either to getGumbelGnuplotFormula or getGaussGnuplotFormula depending on whether one uses the gumbel or the gaussian distribution for incorrectly assigned sequences.

getPositiveGnuplotFormula_

const std::string(PosteriorErrorProbabilityModel::* getPositiveGnuplotFormula_) (const GaussFitter::GaussFitResult &params) const

private

points to getGumbelGnuplotFormula

incorrect_is_gumbel_

bool incorrect_is_gumbel_

private

whether the incorrect component of the (last) mixture model is a Gumbel distribution (otherwise Gauss)

incorrectly_assigned_fit_gumbel_param_

GumbelMaxLikelihoodFitter::GumbelDistributionFitResult incorrectly_assigned_fit_gumbel_param_

private

incorrectly_assigned_fit_param_

GaussFitter::GaussFitResult incorrectly_assigned_fit_param_

private

stores parameters for incorrectly assigned sequences. If gumbel fit was used, A can be ignored. Furthermore, in this case, x0 and sigma are the local parameter alpha and scale parameter beta, respectively.

isotonic_peps_

std::vector<double> isotonic_peps_

private

isotonic_scores_

std::vector<double> isotonic_scores_

private

knots of the non-parametric PEP curve: scores (ascending) and PEPs (non-increasing)

model_type_

ModelType model_type_

private

kind of the last estimated model

negative_prior_

double negative_prior_

private

stores final prior probability for negative peptides

smallest_score_

double smallest_score_

private

smallest score which was used for fitting the model

td_auc_

double td_auc_

private

area under the ROC curve of targets vs. decoys of the last target-decoy fit (-1 if not available)

td_incorrect_fraction_

double td_incorrect_fraction_

private

target-decoy estimate of the fraction of incorrect targets (-1 if not available)