Percentile calibration for all predictors

Affinity, processing, and presentation share two implementations: HistogramPercentRankTransform preserves the historical histogram lookup; CompactPercentRankTransform uses a continuous, compact approximation to the reference distribution.

New predictor calibration defaults to compact, starting with 64 knots and allowing 128 only when background validation justifies it. Loading an existing predictor preserves its saved calibration, including public models. Calibration changes neither neural-network weights nor raw prediction scores.

What a percentile means

Predictor percentile_ranks(...) methods return values from 0 to 100, with lower always meaning stronger. A percentile estimates how frequently the chosen background produces an equally strong or stronger score. It is not a probability that a peptide binds, is processed, or is presented.

Predictor

Raw score and direction

Compact input coordinate

Low-level transform call

Affinity

IC50 in nM; lower is stronger

log(IC50)

transform(ic50)

Processing

Score in [0, 1]; higher is stronger

logit(score)

transform(scores, survival=True)

Presentation

Score in [0, 1]; higher is stronger

logit(score)

transform(scores, survival=True)

Affinity calibration is per allele, with existing sequence-equivalent allele lookup retained. Processing and presentation each store one score distribution per predictor. A processing percentile is not an affinity percentile, and a presentation percentile is not an average of its component percentiles.

Functional form: what the knots store

For each distinct background score v, the compact fit computes the midrank survival probability, assigning half the tied reference mass to either side:

\[S(v) = \frac{\#\{v_j > v\} + \tfrac12\#\{v_j = v\}}{N}.\]

Each knot stores a pair (x_i, y_i): the transformed raw score and logit(S(v_i)). Logarithms are natural logarithms; logit(p) = log(p / (1-p)). Between adjacent knots the ordinate is linear:

\[y(x) = y_i + \frac{x-x_i}{x_{i+1}-x_i}(y_{i+1}-y_i), \qquad x_i \leq x \leq x_{i+1}.\]

With sigmoid(z) = 1 / (1 + exp(-z)), the two directions are:

\[P_{\mathrm{processing/presentation}}(v) = 100\,\operatorname{sigmoid}(y(x(v))), \qquad P_{\mathrm{affinity}}(v) = 100\,\operatorname{sigmoid}(-y(x(v))).\]

The implementation evaluates each direction directly; it does not obtain a tiny upper-tail percentile by subtracting an already rounded CDF from 100. This is a continuous, piecewise-linear curve in transformed coordinates, not a histogram, a cubic spline, or a new neural network. Its slope can change at a knot, so it is not generally differentiable there.

The fit keeps both endpoints and repeatedly adds the reference point with the largest absolute interpolation error in y. Equal errors choose the earliest source index. It stops at the budget, or earlier if there are too few distinct coordinates or all remaining errors are at most 1e-12. Thus 64 is a knot budget, not a requirement to store 64 redundant points. Storage is two coordinates per knot, two extrapolation slopes, and reference/selection metadata.

Outside the reference range, a line anchored at the corresponding endpoint continues the curve. Its slope is fitted using up to 32 distinct reference coordinates nearest that endpoint, rather than just the last knot interval. Those extrapolated percentiles are model-based estimates: continuity does not establish the accuracy of frequencies in a sparsely observed tail.

Choosing 64 versus 128

There are two distinct APIs:

  • Low-level CompactPercentRankTransform(...).fit(scores, num_knots=64) fits the specified budget directly, without automatic selection.

  • Predictor calibrate_percentile_ranks(..., max_knots=128) and the shared fit_percent_rank_transform helper start with 64 and may select 128. The default value 128 is a ceiling, not a starting size.

The automatic rule uses background scores only, never evaluation labels:

  1. With at least 1,000 groups, hold out one fifth of the groups using a fixed selection seed of 403. Peptide identities define groups in the affinity APIs and calibration commands. Processing/presentation score-only APIs use rows unless the caller supplies groups; pass identities when rows repeat.

  2. Compare 64- and 128-knot fits on the remaining groups. Evaluate percentile cutoffs 0.03, 0.1, 0.3, 1, 3, and 10 percent, retaining only cutoffs with at least ten expected validation observations. At least two cutoffs and two distinct fitting scores are required; otherwise retain 64.

  3. For a cutoff q, let c_q count validation percentiles at or below q and M be the validation row count. The smoothed observed percentage is 100 * (c_q + 0.5) / (M + 1). For each budget, E is the root-mean-square of log10(observed_percentage / q) over the eligible cutoffs.

  4. Select 128 only if E64 > E128 + max(0.01, 0.1 * min(E64, E128)). Fit the selected budget on the complete background.

Setting max_knots=64 disables escalation. Fewer than 1,000 groups also retain 64. This is a predeclared practical tolerance, not a statistical significance test or a guarantee of accuracy at every percentile.

The compact transform’s selection dictionary records the selected budget, actual knot count, reason, seed, reference count/hash, and—when comparison is possible—validation counts, split-mask hash, cutoffs, errors, and tolerance. The fixed selection seed is separate from the CLI --random-seed (default 42), which governs generation of calibration peptides and MHC allele sets.

Python examples

These small synthetic arrays illustrate the API, not a suitable release calibration background. Both classes use instance fit, return self, and support to_dict/from_dict serialization.

>>> import json
>>> import numpy as np
>>> from mhcflurry import HistogramPercentRankTransform, CompactPercentRankTransform
>>> from mhcflurry.percent_rank_transform import PercentRankTransform
>>> PercentRankTransform is HistogramPercentRankTransform
True
>>> background_ic50 = np.array([1., 10., 100., 1000.])
>>> affinity_curve = CompactPercentRankTransform(score_transform="log").fit(background_ic50)
>>> np.round(affinity_curve.transform(background_ic50), 2).tolist()
[12.5, 37.5, 62.5, 87.5]
>>> probability_curve = CompactPercentRankTransform(score_transform="logit").fit([.1, .1, .2, .9])
>>> np.round(probability_curve.transform([.1, .2, .9], survival=True), 2).tolist()
[75.0, 37.5, 12.5]
>>> restored = CompactPercentRankTransform.from_dict(json.loads(json.dumps(probability_curve.to_dict())))
>>> np.array_equal(restored.transform([.2, .8]), probability_curve.transform([.2, .8]))
True
>>> histogram = HistogramPercentRankTransform().fit(background_ic50, bins=[1., 10., 100., 1000.])
>>> np.array_equal(histogram.transform([5., 50.], survival=True), 100 - histogram.transform([5., 50.]))
True

For predictor calibration, affinity takes background peptides and performs inference; processing and presentation take already computed background scores. Their percentile_ranks methods all take raw scores:

# Assume independently prepared backgrounds and loaded predictor instances.
affinity_predictor.calibrate_percentile_ranks(
    peptides=background_peptides, alleles=calibration_alleles,
    method="compact", max_knots=128)
affinity_percentiles = affinity_predictor.percentile_ranks(query_ic50, allele=allele)

processing_predictor.calibrate_percentile_ranks(
    background_processing_scores, groups=background_peptides,
    method="compact", max_knots=128)
processing_percentiles = processing_predictor.percentile_ranks(query_processing_scores)

presentation_predictor.calibrate_percentile_ranks(
    background_presentation_scores, groups=background_peptide_per_score,
    method="compact", max_knots=128)
presentation_percentiles = presentation_predictor.percentile_ranks(query_presentation_scores)

Processing’s raw predict/predict_to_dataframe outputs are unchanged; call percentile_ranks explicitly. An uncalibrated processing predictor raises by default; throw=False warns and returns NaNs. Do not assume historical public processing models include a calibration merely because the API now exists.

Use method="histogram" to explicitly request histogram calibration. Supplying bins also selects it when method is omitted; method="compact" plus bins is rejected. Exact historical reproduction requires the original reference and bin edges, not just the histogram class name. In particular, the current presentation histogram default uses the later tail-adaptive bin policy, not necessarily the grid in an older public artifact.

CompactPresentationPercentiles is only an adapter for the original experiment: it retains its historical classmethod fit and upper-tail default. New code should use the generic class and specify survival=True where appropriate.

Calibration commands and reference policy

The command writes in place. Point it at a separately preserved candidate copy, not a public baseline or the only copy of an earlier calibration. These are command templates; retain the original release’s reference policy, sample count, amino-acid distribution, seed, and MHC allele set sampling for a controlled recalibration. Omitted affinity/presentation reference-generation options use CLI defaults, which need not match a release’s recipe.

mhcflurry calibrate-percentile-ranks \
    --models-dir candidate-affinity-copy \
    --predictor-kind class1_affinity \
    --percentile-method compact --max-percentile-knots 128 \
    --random-seed 42 --num-jobs 0

mhcflurry calibrate-percentile-ranks \
    --models-dir candidate-processing-copy \
    --predictor-kind class1_processing \
    --processing-reference-data independent_background.csv \
    --percentile-method compact --max-percentile-knots 128 \
    --num-jobs 0

mhcflurry calibrate-percentile-ranks \
    --models-dir candidate-presentation-copy \
    --predictor-kind class1_presentation \
    --percentile-method compact --max-percentile-knots 128 \
    --random-seed 42 --num-jobs 0

Processing requires an explicit independent background CSV with peptide and, for flanked ensembles, n_flank and c_flank. Empty flank cells mean missing context. Use the same ensemble and context policy as at prediction time; the command does not fabricate flanks or generate a processing reference from --num-peptides-per-length. It currently supports --num-jobs 0 only and retains inputs and raw scores in percent_rank_reference.npz.

Affinity/presentation commands generate background predictions from their configured reference policy; unlike the processing command, they do not automatically retain a raw reference-score archive. Preserve those predictions separately when needed for calibration plots/reanalysis. Retain model hashes, reference inputs, the command, transform files, and held-out prediction tables with joinable sample/peptide/allele identifiers for all comparisons.

Affinity --only-missing fills gaps, including sequence-equivalent lookup; it does not convert existing histogram calibrations to compact. For an intentional conversion of all requested alleles, omit that flag on a copy.

Saved formats and compatibility

  • Historical percent_ranks.csv loads through the histogram implementation, without refitting. PercentRankTransform remains an import/pickle alias.

  • percent_ranks.json stores versioned compact or mixed collections, including different knots per allele. All-histogram collections use CSV when their edges match; different histogram grids also require JSON.

  • Saving replaces the target calibration file and removes its alternate format only after successful replacement. Saving an empty collection removes saved calibration. Loading prefers JSON if both files are present; it does not fall back to CSV if the JSON is invalid.

  • Older software cannot interpret compact JSON calibration. Do not expect percentile-dependent prediction to work there merely because the network weight files remain readable. Preserve original artifacts for old clients.

Changing ensemble members or weights requires new calibration. Keep raw-score performance comparisons separate from percentile-mapping comparisons.

Numerical limits and interpretation

Compact fitting requires a finite, nonconstant reference. Log coordinates require positive scores; logit coordinates require scores in [0, 1]. Exact zero and one use the nearest interior float64 values for taking logits. Transform calls propagate NaNs but reject infinities and out-of-domain scores. The histogram implementation retains its historical input/boundary behavior.

Compact interpolation is monotone within one distribution, avoiding the deliberate step ties introduced by histograms. It cannot recover ties already present in raw scores or remove floating-point rounding. Extreme percentiles can still underflow to zero or approach 100 without distinguishable digits; log_percentiles(scores, survival=...) retains tiny-tail information before exponentiation, but does not validate extrapolated frequencies.

Monotonic mapping preserves ordering within a shared calibration distribution, subject to those numerical limits. Combining allele-specific calibrations can change ordering across alleles. Better percentile resolution is not an improvement in the neural network’s raw ranking ability, and calibration is not a substitute for held-out affinity/processing/presentation evaluation.

See Presentation percentile calibration for the archived presentation comparison and its distinct 60/20/20 experimental validation design.