Recursive confirmation bias in MHC class-I epitope prediction

Code, model weights and data for the recursive-data-corruption experiment: a family of peptide–MHC ranking models trained on progressively corrupted label sets, showing that a model's apparent performance climbs with each round of corruption while its true performance does not.

Reproduce the result

git clone https://github.com/deepflare/IEDB_RCB && cd IEDB_RCB
pip install -r requirements.txt
python reproduce.py

That regenerates results/corruption_figure.png and both metric tables from results/evaluation_data.parquet, which is committed here. No downloads, no credentials, no GPU. On a 2023 MacBook Pro (M2, 16 GB) the clone takes about 10 seconds, pip install about 70 seconds, and reproduce.py about 20 seconds.

The run prints two metric tables. The values that appear in the paper are the iteration-3 row of each series: sensitivity@2% of 0.585 for corrupted against 0.231 for golden, and AUROC of 0.991 for corrupted. The clean_baseline control sits at 0.205 and 0.888. If your numbers match these, the reproduction succeeded.

Red is each model measured against the corrupted labels it was trained to agree with; it rises with every round of corruption. Green is the same model measured against clean labels. Grey is the never-corrupted iteration-0 model on those same clean folds. Green tracking grey is the point: the apparent gain is an artefact of the corrupted evaluation labels, not a real improvement.

System requirements

No non-standard hardware is required. Everything runs on a normal desktop or laptop under any operating system with a working Python; a GPU is optional and is used automatically for scoring if present.

Python	3.11 or newer
Dependencies	requirements.txt — floors are the supported minimums
Disk	~110 MB for the repository, ~1.2 GB more for the full environment
Tested on	macOS 15.6 (Apple silicon) and Linux x86-64

Tested against two dependency sets: the floors in requirements.txt on Python 3.11, and current releases on Python 3.13 (torch 2.13, transformers 5.14, pandas 3.0, numpy 2.5, scipy 1.18, matplotlib 3.11, pyarrow 25). Both reproduce the reported values.

reproduce.py itself needs only pandas, pyarrow, numpy, scipy and matplotlib — it does not import PyTorch. Installing just those five takes about 30 seconds and 410 MB, against 70 seconds and 1.2 GB for the full requirements.txt, and reproduces every reported value identically. Torch and transformers are needed only for predict.py and rank_protein.py.

pip install pandas pyarrow numpy scipy matplotlib && python reproduce.py

The experiment

Five independent dataset versions. Within each, an iteration ladder:

• iteration 0 — trained on the clean seed. Three hyperparameter-sweep replicas.
• iterations 1–3 — each trained on labels rewritten by the previous iteration's own predictions: a positive is kept only if that model ranked it in the top 2% and scored it ≥ 0.9, and is flipped to negative otherwise. Bias compounds.

Six models per version, 30 in total, all listed in models.csv.

Score your own peptides

python predict.py --model MHCRANK-9287 \
    --pairs GILGFVFTL,HLA-A:0002:0001 NLVPMVATV,HLA-A:0002:0001

This needs the checkpoints, which are too large for git. Download checkpoints.zip from the Zenodo deposit at 10.5281/zenodo.21641350 and unpack it at the repository root; see DATA.md. The first run fetches the public ESM-2 config and tokenizer from Hugging Face (a few MB); all model weights come from the checkpoint.

Rank the peptides in a protein

To go from a whole antigen to candidate epitopes, rank_protein.py slides a window over the sequence and scores every candidate:

python rank_protein.py --fasta antigen.fasta \
    --alleles HLA-A:0002:0001 --lengths 9 10 --top 20

Output is ordered best first, with each peptide's position in the protein. Default lengths are 9–12; the models were trained on 8–15mers, so any length in that range can be requested via --lengths.

Scores are comparable within a run — every candidate is padded to the same width, so the ranking does not shift with --batch-size. They are probabilities from a model trained on one dataset, not calibrated affinities, so the ordering is the output rather than the absolute value.

As a sanity check, ranking all 4,452 8–14mers of HSPA8 (P11142) against HLA-A*30:01 puts every peptide observed for that protein and allele in the validation set inside the top 16%, with a median rank of 300. The highest-scoring candidates almost all carry a C-terminal lysine, which is the anchor residue A*30:01 prefers.

Hardware and allele coverage

Both predict.py and rank_protein.py run on CPU and use a GPU automatically when one is available, selecting CUDA or Apple MPS. Force it with --device cpu|cuda|mps. Scores agree across devices to about 1e-6.

Alleles come from mhc_alpha_chains.fasta, which ships 92 class-I alpha chains (31 A, 40 B, 18 C, 3 G). To score an allele that is not listed, add its alpha-chain sequence to that file under the same naming convention.

Architecture

The architecture is a facebook/esm2_t6_8M_UR50D encoder over the joined sequence <cls> MHC [SEP] peptide <eos>, four concatenated pooled representations of the final hidden state (mean, mean/√L, max, cls), and a linear head over two classes. The score is softmax(logits)[:, 1]. MHC input is the 182-residue alpha chain in mhc_alpha_chains.fasta.

Layout

Path	Contents
mhcrank.py	model architecture and checkpoint loading
metrics.py	AUROC and sensitivity@2%, with tie handling
predict.py	score peptide–MHC pairs
rank_protein.py	rank every peptide in a protein for a given allele
reproduce.py	recompute the metrics and redraw the figure
models.csv	the 30 released models, with checkpoint names and SHA-256
results/	committed metrics, the figure, and the labels and scores behind them
datasets/datasets.csv	one row per dataset bundle; bundles themselves listed in DATA.md
predictions/	predictions.csv index; tables listed in DATA.md
mhc_alpha_chains.fasta	MHC alpha-chain sequences used as model input

Dataset names read golden_val_v1_fold0 (clean) and val_corrupted_v1_iter1 (corrupted), where v1–v5 is the dataset version.

Metric definitions

metrics.py computes both metrics from mid-ranks, so tied scores contribute equally regardless of row order. That matters here: the models saturate near 0 and 1 and produce large blocks of identical scores, and resolving those by row order would make the numbers depend on how the data happened to be stored. Sensitivity@2% resolves scores tied at the cutoff proportionally, reporting the expected recall under random tie-breaking.

Uncertainty in the figure is the 95% t-interval across the five dataset versions, which are the independent replicates. The corrupted series has a genuinely narrow interval — the models agree closely about how good they look on their own labels.

Reproducibility notes

predict.py reproduces the released prediction tables to a mean absolute difference of about 2e-3 per score, which moves AUROC by less than 1e-5. The original runs used reduced-precision GPU arithmetic; this code runs full precision on CPU. The prediction tables in the Zenodo deposit are the original outputs, and reproduce.py --from-raw recomputes every reported metric from them.

The pooling feeding the classification head spans the full padded width of a batch, which is how the models were trained. Scores therefore depend slightly on batch padding, and predict.py preserves that training-time behaviour rather than masking padding, so inference matches training.

What is not here

The training harness that produced the checkpoints is not included. Every reported number is reproducible from the released weights, datasets and predictions without it. It is available to editors and reviewers on request.

Citing

See CITATION.cff. Code is MIT licensed; see LICENSE. Underlying observations derive from the Immune Epitope Database and from previously published mass-spectrometry studies, which remain available from their original sources.