# Import data loading functions
from mhcpred.data import get_train_data, get_test_data
# Load training and test data
df_train = get_train_data()
df_test = get_test_data()This post trains a neural network to predict which peptides are presented by a given MHC class I molecule, and benchmarks it against MHCflurry. Section 1 introduces the biology: what MHC molecules are, how classes I and II differ, and why the problem suits machine learning. Section 2 covers the encodings, the model, and the results.
Understanding the major histocompatibility complex (MHC)
The major histocompatibility complex (MHC) is a region of DNA whose genes code for a family of cell-surface proteins, the MHC molecules. These molecules bind fragments of proteins (antigens) and display them on the cell surface, where T cells can inspect them. This display is how the immune system distinguishes the body’s own cells from infected, cancerous, or foreign ones, and it is what triggers an immune response when something looks wrong.
MHC molecules matter for:
- Distinguishing self from non-self, which prevents autoimmune attacks.
- Activating the immune response against infected cells.
- Organ transplantation, where they determine donor compatibility.
- Susceptibility to autoimmune disease, since some MHC variants carry increased risk.
Note
MHC molecules are cellular “display cases”: they present antigens to T cells, triggering an immune response when necessary.
MHC class I and class II
MHC molecules come in two main classes, distinguished by where the antigens they display come from.
MHC class I
Class I molecules are present on most cells, red blood cells excepted, and they present antigens from inside the cell. When a cell is infected or becomes cancerous, proteins within it are broken down into short fragments called epitopes. These are loaded onto class I molecules and displayed at the surface, where killer T cells (cytotoxic T lymphocytes) recognise them and can destroy the cell.
In humans, the main class I molecules are HLA-A, HLA-B and HLA-C, where HLA stands for human leukocyte antigen.
MHC class II
Class II molecules are found mainly on antigen-presenting cells (APCs) such as macrophages, dendritic cells and B cells, and they present antigens from outside the cell. APCs engulf foreign invaders by phagocytosis and break them into epitopes, which are then loaded onto class II molecules. A helper T cell that recognises a foreign epitope becomes activated and releases cytokines, recruiting B cells and killer T cells to fight the infection.
In humans, the main class II molecules are HLA-DP, HLA-DQ and HLA-DR.
Human leukocyte antigens (HLA)
In humans, MHC molecules are called human leukocyte antigens, and the genes coding for them sit in the MHC region on chromosome 6. They matter for:
- Organ transplantation, where HLA matching prevents rejection.
- Autoimmune disease, since certain HLA types carry increased risk.
- Drug response, which HLA variation can influence.
- Population-level resilience, since HLA diversity protects against a wider range of pathogens.
MHC diversity
The MHC is highly polymorphic: each MHC gene exists in many versions (alleles), and different MHC molecules bind different peptides. A diverse population is therefore more likely to contain individuals able to present antigens from a new pathogen. That same polymorphism is what makes the prediction problem hard, since a useful model has to generalise across alleles and not just across peptides.
Public data sources
Two public resources anchor most work in this area. The IPD-IMGT/HLA database curates the sequences of the human MHC/HLA system: allele sequences, nomenclature and associated metadata. Accurate HLA typing depends on it, which makes it central to immunology, transplantation and vaccine development. The Immune Epitope Database (IEDB) catalogues curated experimental data on epitopes recognised by T cells and B cells across diseases and conditions, and supports epitope discovery and vaccine design.
The project below uses neither directly. It works from the NetMHCpan-4.1 training and evaluation data (Reynisson et al. 2020), described in Section 2.
Predicting peptide presentation: a machine learning approach
Relevance to machine learning
Whether an MHC molecule presents a peptide depends on the peptide’s amino acid sequence and on the particular MHC allele. That sequence-to-function relationship, together with experimental data measuring it, is what makes the problem tractable for machine learning. The difficulty is the polymorphism: with many alleles in the population and very uneven amounts of data across them, a model that treats each allele as an independent category has nothing to say about an allele it has never seen.
Key concepts
- Epitope: the part of the peptide recognised by the T cell receptor.
- Polymorphism: the existence of many alleles of the MHC genes within a population, each binding a different set of peptides.
- Binding affinity: the strength of the MHC-peptide interaction, usually reported as an IC50 in nM, where a lower value means tighter binding.
- Presentation: whether the peptide is actually found displayed on the cell surface, typically measured by mass spectrometry of eluted ligands. This is the target used in this project.
Modelling choices
The task can be framed as regression on affinity or as binary classification.
- Input: a peptide sequence and an MHC allele.
- Output: a binding affinity (IC50, Kd), or a binary label — binder/non-binder, or, as here, presented/not presented.
Common approaches:
- Sequence-based features: amino acid composition, n-grams, physicochemical properties.
- Structure-based features: information about the 3D structure of the MHC-peptide complex, when available.
- Allele encoding: one-hot encoding, learned amino acid embeddings, or substitution-matrix encodings of the residues lining the binding groove.
- Algorithms: linear models, support vector machines, random forests, and neural networks including convolutional networks and transformers.
Project overview
Note
The code is available at https://github.com/nbrosse/mhcpred.
The goal is a binary classifier that predicts whether a given peptide is presented by a given MHC class I allele. The data comes from the training and evaluation sets of NetMHCpan-4.1 (Reynisson et al. 2020), a well-established framework for MHC binding prediction. The training data ships as five folds; this project concatenates them and carves out a validation set with a single stratified 90/10 split, so no cross-validation is performed.
The target is the binary hit column, 1 if the peptide is presented by the MHC and 0 otherwise, and there are two features:
peptide: the amino acid sequence of the peptide. These short chains are the potential antigens that could be presented to the immune system.allele: the name of the MHC class I allele. Because MHC molecules are highly polymorphic, each allele has a slightly different binding groove and therefore presents a different set of peptides. The naming convention is documented here.
Note
MHC antigen presentation is a large field, and this project is a deliberately simplified entry point. For a serious treatment, see NetMHCpan (Reynisson et al. 2020) and MHCflurry (O’Donnell et al. 2020) and the references cited within them. The NetMHCpan-4.1 data used here must remain private, so only aggregate statistics and sample rows are shown.
We begin with exploratory data analysis.
Two facts from this analysis drive everything that follows. First, the data is heavily imbalanced — 5.4% positives in training, 4.8% in test — so raw accuracy is close to meaningless and balanced accuracy is the metric to watch. Second, the per-allele counts span five orders of magnitude, from 265,252 samples for HLA-A*02:01 down to 6, and four alleles appear only in the test set. Encoding an allele by its sequence rather than its name is what makes any prediction for those four possible at all.
Using MHCflurry pretrained models for prediction
We use the mhcflurry package as a benchmark; see the associated paper (O’Donnell et al. 2020). MHCflurry predicts how strongly peptides bind MHC class I molecules, using neural networks trained on a large set of experimentally measured binding affinities together with mass spectrometry data on ligand presentation.
The baseline used here is MHCflurry’s binding affinity model, Class1AffinityPredictor. The following assumes mhcflurry is installed and the pretrained models downloaded.
mhcflurry-downloads fetch models_class1_presentation
python scripts/mhcflurry_benchmark.pydef predict_with_mhcflurry() -> pd.DataFrame:
predictor = Class1AffinityPredictor.load()
df_test = get_test_data()
mhcflurry_predictions = predictor.predict_to_dataframe(
peptides=df_test.peptide.values,
alleles=df_test.allele.values,
allele=None,
)
df = pd.merge(df_test, mhcflurry_predictions, on=["allele", "peptide"], how="left")
df.to_csv(str(output_path / "mhcflurry_predictions.csv"), index=False)
return dfThe output is of the form:
| peptide | hit | allele | prediction | prediction_low | prediction_high | prediction_percentile |
|---|---|---|---|---|---|---|
| AAPATRAAL | True | HLA-B*35:03 | 94.297 | 59.902 | 144.624 | 0.205 |
| AAPSAAREL | True | HLA-B*35:03 | 116.19 | 79.847 | 169.241 | 0.262 |
| AEISQIHQSVTD | True | HLA-B*35:03 | 26103.26 | 22695.389 | 28415 | 15.739 |
| ALEEQLQQIRAE | True | HLA-B*35:03 | 24797.131 | 19988.967 | 28062.65 | 13.571 |
| AQDPLLLQM | True | HLA-B*35:03 | 2164.336 | 745.888 | 5390.727 | 1.413 |
| ASAPPGPPA | True | HLA-B*35:03 | 1398.729 | 387.675 | 3293.692 | 1.157 |
| DAHKGVAL | True | HLA-B*35:03 | 84.315 | 54.736 | 133.899 | 0.175 |
| DNPIQTVSL | True | HLA-B*35:03 | 1386.767 | 565.122 | 3667.21 | 1.151 |
| DPEAFLVQI | True | HLA-B*35:03 | 245.485 | 133.986 | 394.752 | 0.484 |
The first three columns come from the test dataset.
- peptide: the amino acid sequence of the peptide being evaluated.
- hit: the ground truth, indicating whether the peptide is known to be presented by the given MHC allele (
True) or not (False). - allele: the name of the MHC class I allele being considered.
The remaining columns are added by the binding affinity model.
- prediction: the predicted IC50 in nM. Lower means tighter binding: in Table 1,
AAPATRAALat 94 nM is a strong binder whileAEISQIHQSVTDat 26,103 nM is not. Values below roughly 500 nM are conventionally treated as binders. - prediction_low / prediction_high: the 5th and 95th percentiles of the predictions made by the individual networks in the ensemble. They measure how much the ensemble disagrees, which is a useful uncertainty proxy but not a calibrated confidence interval.
- prediction_percentile: the rank of the predicted affinity against a background distribution of scores for random peptides on the same allele. Lower is stronger: a percentile of 1.0 means the peptide scores in the top 1% for that allele. This is the column to threshold on, because it is comparable across alleles while raw nM values are not.
A percentile threshold of 2% is the usual cutoff for calling a peptide a likely binder, and is what the evaluation below uses.
Warning
The benchmark is not entirely like for like. The test labels are mass spectrometry presentation hits, but the baseline is MHCflurry’s affinity predictor thresholded at the 2% percentile. MHCflurry 2.0 also ships a Class1PresentationPredictor, which combines affinity with an antigen-processing model and is the appropriate comparison for a presentation task. The numbers reported below therefore understate MHCflurry somewhat.
Prediction: fitting a Class1BinaryNeuralNetwork
We now fit a Class1BinaryNeuralNetwork on the training dataset. The code is available at https://github.com/nbrosse/mhcpred.
Class1BinaryNeuralNetwork subclasses MHCflurry’s Class1NeuralNetwork and overrides only fit_generator and predict. The architecture is therefore MHCflurry’s, but the model is trained from scratch on the binary hit target with a binary cross-entropy loss, and predict returns the sigmoid output — a probability of presentation — rather than an nM affinity.
Here is a glimpse of the training data structure:
| peptide | allele | hit |
|---|---|---|
| YFPLAPFNQL | HLA-C*14:02 | True |
| KESKINQVF | HLA-B*44:02 | True |
| QPHDPLVPLSA | HLA-B*54:01 | True |
| RTIADSLINSF | HLA-B*57:03 | True |
Encoding is the main obstacle: peptides have variable lengths, and alleles arrive as names rather than sequences.
The allele half is solved by allele_sequences.csv, which ships with the MHCflurry model data and maps each allele to the 37 residues that line its binding groove — the pseudo-sequence that actually contacts the peptide.
| Allele | Sequence |
|---|---|
| HLA-A*01:01 | YFAMYQENMAHTDANTLYGIIYDRDYTWVARVYRGYA |
| HLA-A*01:02 | YSAMYQENMAHTDANTLYGIIYDRDYTWVARVYRGYA |
| HLA-A*01:03 | YFAMYQENMAHTDANTLYGIMYDRDYTWVARVYRGYA |
| HLA-A*01:04 | YFAMYQENMAHTDANTLYGIIYDRDYTWVARVYRGYX |
| HLA-A*01:06 | YFAMYQENMAHTDANTLYGIIYDRDYTWVALAYRGYA |
First, we import the necessary libraries, including components from our own mhcpred library, which holds the network and the data loading functions.
import pickle
from pathlib import Path
from typing import Iterator
import numpy as np
import pandas as pd
from mhcflurry.allele_encoding import AlleleEncoding
from mhcflurry.encodable_sequences import EncodableSequences
from sklearn.model_selection import train_test_split
from mhcpred.class1_binary_nn import Class1BinaryNeuralNetwork
from mhcpred.config import settings
from mhcpred.data import get_test_data, get_train_data
from mhcpred.hyperparameters import base_hyperparametersWe load the allele sequences, the training data and the test data.
allele_sequences = pd.read_csv(
str(data_path / "allele_sequences.csv"), index_col=0
).iloc[:, 0]
df_total_train = get_train_data()
df_test = get_test_data()We then restrict the allele sequences to the alleles that actually occur in our data.
alleles_in_use = set(df_total_train.allele).union(set(df_test.allele))
allele_sequences_in_use = allele_sequences[allele_sequences.index.isin(alleles_in_use)]Two MHCflurry classes do the encoding work. The AlleleEncoding class maps allele names to integer indices and to their pseudo-sequences, and caches the encodings so the same mapping is reused everywhere. The EncodableSequences class turns variable-length peptides into fixed-size numerical matrices.
We split the training data with a stratified 90/10 split, which keeps the class balance in both parts, and encode the validation peptides and alleles.
allele_encoding = AlleleEncoding(
alleles=allele_sequences_in_use.index.values,
allele_to_sequence=allele_sequences_in_use.to_dict(),
)
df_train, df_val = train_test_split(
df_total_train, test_size=0.1, shuffle=True, stratify=df_total_train.hit.values
)
val_peptides = EncodableSequences(df_val.peptide.values)
val_alleles = AlleleEncoding(
alleles=df_val.allele.values,
allele_to_sequence=allele_sequences_in_use.to_dict(),
)AlleleEncoding is worth a closer look, because it is what lets the model say anything about an allele it never trained on.
Allele universe versus used alleles. The class separates the complete set of alleles it knows about, defined by the
allele_to_sequencedictionary, from the specific alleles used in a given task, passed as a list at construction time. Theallele_to_indexmapping covers the whole universe, including a special index forNoneused as padding, so a given allele always receives the same index.Padded sequence storage. The pseudo-sequences for the whole universe are stored in a pandas
Seriesand padded to a common length with theXcharacter, which is what makes a fixed-length numerical representation possible.Borrowing. The
borrow_fromparameter creates a newAlleleEncodingthat inherits the universe and the index mapping from an existing instance, so the per-batch encodings built inside the training loop stay consistent with the global one without redeclaring the mapping.Encoding.
allele_representations(encoding_name)encodes the entire universe once, andfixed_length_vector_encoded_sequences(encoding_name)selects the rows for the alleles in use, in order. Theencoding_nameselects the scheme, typicallyBLOSUM62or one-hot.
BLOSUM62 (Blocks Substitution Matrix) is a standard substitution matrix in bioinformatics: for each pair of amino acids it gives a score reflecting how likely one is to be substituted for the other over evolutionary time, with higher scores for more interchangeable pairs and negative scores for unfavourable substitutions. To encode a sequence, each amino acid is replaced by its row of the matrix, a vector of 21 numbers covering the 20 amino acids plus the padding character X. A sequence of length n therefore becomes an n × 21 matrix. Using BLOSUM62 rather than one-hot means chemically similar residues start out with similar representations, which helps when data for a given allele is thin.
The train_data_iterator function yields batches of training data, and drops any allele that appears in the training data but has no sequence available.
def train_data_iterator(
df_train: pd.DataFrame,
train_allele_encoding: AlleleEncoding,
batch_size: int = 1024,
) -> Iterator[tuple[AlleleEncoding, EncodableSequences, np.ndarray]]:
"""
This function creates a data generator for training the neural network.
It iterates over the training data in batches and yields tuples of
(allele_encoding, peptide_sequences, labels). It also handles filtering
of alleles not found in the initial allele encoding.
"""
# Get unique alleles in the training set.
alleles = df_train.allele.unique()
# Filter alleles to keep only those for which sequences are available.
usable_alleles = [
c for c in alleles if c in train_allele_encoding.allele_to_sequence
]
print("Using %d / %d alleles" % (len(usable_alleles), len(alleles)))
print(
"Skipped alleles: ",
[c for c in alleles if c not in train_allele_encoding.allele_to_sequence],
)
df_train = df_train.query("allele in @usable_alleles")
# Calculate the number of batches.
n_splits = np.ceil(len(df_train) / batch_size)
# Infinite loop to allow for multiple epochs.
while True:
# Split the training data into batches.
epoch_dfs = np.array_split(df_train.copy(), n_splits)
for k, df in enumerate(epoch_dfs):
if len(df) == 0:
continue
# Encode peptides and alleles for the current batch.
encodable_peptides = EncodableSequences(df.peptide.values)
allele_encoding = AlleleEncoding(
alleles=df.allele.values,
borrow_from=train_allele_encoding, # Reuse encoding from main allele_encoding
)
# Yield the encoded data and labels (hit column).
yield (allele_encoding, encodable_peptides, df.hit.values)The model is initialised from the base hyperparameters and trained with fit_generator, which takes the generator, the validation data, and the number of epochs and steps per epoch.
batch_size = 1024
train_generator = train_data_iterator(
df_train=df_train,
train_allele_encoding=allele_encoding,
batch_size=batch_size,
)
model = Class1BinaryNeuralNetwork(**base_hyperparameters)
steps_per_epoch = np.ceil(len(df_train) / batch_size)
model.fit_generator(
generator=train_generator,
validation_peptide_encoding=val_peptides,
validation_affinities=df_val.hit.values,
validation_allele_encoding=val_alleles,
validation_inequalities=None,
validation_output_indices=None,
steps_per_epoch=steps_per_epoch,
epochs=2,
)Two things about this call are worth flagging, since they explain the results later. epochs=2 is a deliberate shortcut, far short of convergence. And fit_generator does not support the random negative peptide sampling that fit offers, so the random_negative_* entries in the hyperparameters are inert here and the 5.4% positive rate is fed to the network untouched.
The network takes two inputs.
- Peptide: a 45 × 21 matrix. Peptides are between 5 and 15 amino acids long, and the
left_pad_centered_right_padscheme writes each peptide three times into a 45-slot array — left-aligned, centred, and right-aligned — padding the remaining slots withX. Repeating the peptide in three alignments lets the network pick up motifs anchored at the N-terminus, at the C-terminus, or in the middle without committing to a single convention. Each of the 45 slots is then a 21-dimensional BLOSUM62 vector. - Allele: a single integer index into a lookup table.
These are processed as follows.
Allele representation layer. The allele index goes through a Keras
Embeddinglayer, but the layer is created withtrainable=False: its weights are set to the BLOSUM62 encoding of each allele’s 37-residue pseudo-sequence, 37 × 21 = 777 values per allele, and are never updated by training. It is a frozen lookup table rather than a learned embedding, and that is precisely what allows predictions for alleles absent from the training data — the model has their sequence even without their examples.Flatten layers. The 45 × 21 peptide matrix is flattened into a 945-element vector, and the 37 × 21 allele representation into a 777-element vector.
Concatenate layer. The two are joined into a single 1,722-element vector, which is the step that lets the following layers learn the joint effect of allele and peptide.
Dense layers. Two fully connected layers with
tanhactivations, 1,024 units then 512, learn the non-linear interactions between the peptide and the binding groove.Dropout layers. During training, dropout randomly ignores a fraction of the units, here with probability 0.5, which limits overfitting.
Output layer. A single unit with a sigmoid activation, giving the probability that the peptide is presented by the allele.
Finally, the trained model is applied to the test data, which is encoded exactly as the training data was.
test_peptides = df_test.peptide.values
test_allele_encoding = AlleleEncoding(
alleles=df_test.allele.values,
allele_to_sequence=allele_sequences_in_use.to_dict(),
)
predictions = model.predict(
peptides=test_peptides,
allele_encoding=test_allele_encoding,
)
df_test["predictions"] = predictions
df_test.to_csv(str(output_path / "mhcpred_predictions.csv"), index=False)We evaluate both methods with standard binary classification metrics.
Metrics
This notebook contains training metrics history and classification metrics computed on the predictions by - mhcflurry (benchmark) - mhcpred
from pathlib import Path
import pickle
from mhcpred.config import settings
import pandas as pd
from sklearn.metrics import accuracy_score, confusion_matrix, balanced_accuracy_score
from sklearn.metrics import classification_report
from sklearn.metrics import ConfusionMatrixDisplay
import matplotlib.pyplot as plt
models_path = Path(settings.models_path)
output_path = Path(settings.output_path)Information on the training history
I prefer to use tensorboard, but it is not implemented in the mhcflurry package. The information is quite scarce, but when you execute the code, you have the loss for each step and not only for the whole epoch. Of course, it is a very basic version of logging and should be improved.
with open(str(models_path / "model.pickle"), "rb") as f:
model = pickle.load(f)model.fit_info[{'learning_rate': 0.0010000000474974513,
'loss': [0.09700655937194824, 0.06465369462966919],
'val_loss': [0.06880103051662445, 0.05075661838054657],
'time': 524.7155420780182,
'num_points': 6628048}]Binary classification metrics
We compute the usual binary classification metrics on the imbalanced test dataset: accuracy, balanced accuracy, confusion matrix and the scikit-learn classification report.
We rely on balanced accuracy because the dataset is very imbalanced, so plain accuracy is not informative: a model that always predicts False already scores about 0.95.
mhcflurry metrics
mhcflurry_rank_percentile_threshold = 2 # rank threshold for positive hits
# It comes from the mhcflurry article.df = pd.read_csv(str(output_path / "mhcflurry_predictions.csv"))
y_pred = df.prediction_percentile.values <= mhcflurry_rank_percentile_threshold
y_true = df.hit.values
acc = accuracy_score(y_true=y_true, y_pred=y_pred)
confusion_mat = confusion_matrix(y_true=y_true, y_pred=y_pred)
balanced_acc = balanced_accuracy_score(y_true=y_true, y_pred=y_pred)
class_report = classification_report(y_true=y_true, y_pred=y_pred, output_dict=False)
disp = ConfusionMatrixDisplay(confusion_matrix=confusion_mat)
disp.plot()
plt.show()
print(class_report) precision recall f1-score support
False 0.99 0.98 0.99 900996
True 0.67 0.86 0.76 45423
accuracy 0.97 946419
macro avg 0.83 0.92 0.87 946419
weighted avg 0.98 0.97 0.97 946419
acc, balanced_acc(0.9731936911663861, 0.9217833819652606)The metrics are quite good. Precision on the True class is only 0.67, so the model predicts True too often and produces a lot of false positives: 19,234 of them on the confusion matrix.
mhcpred metrics
mhcpred_proba_threshold = 0.5 # by default, but we try to tune it laterdf = pd.read_csv(str(output_path / "mhcpred_predictions.csv"))
y_true = df.hit.values
y_pred = df.predictions.values >= mhcpred_proba_threshold
acc = accuracy_score(y_true=df.hit.values, y_pred=y_pred)
confusion_mat = confusion_matrix(y_true=df.hit.values, y_pred=y_pred)
balanced_acc = balanced_accuracy_score(y_true=df.hit.values, y_pred=y_pred)
class_report = classification_report(y_true=y_true, y_pred=y_pred, output_dict=False)
disp = ConfusionMatrixDisplay(confusion_matrix=confusion_mat)
disp.plot()
plt.show()
acc, balanced_acc(0.9731657332258088, 0.7775326900487307)print(class_report) precision recall f1-score support
False 0.98 0.99 0.99 900725
True 0.82 0.56 0.67 45416
accuracy 0.97 946141
macro avg 0.90 0.78 0.83 946141
weighted avg 0.97 0.97 0.97 946141
mhcpred performs worse than mhcflurry, as the balanced accuracy shows. Here the problem is the opposite one: recall on the True class is only 0.56, so the model predicts False too often and misses about 20,000 presented peptides (the false negatives on the confusion matrix). Lowering the decision threshold should therefore help.
Threshold tuning
We plot the precision recall curve to try to identify a better threshold.
from sklearn.metrics import precision_recall_curve, PrecisionRecallDisplay
precision, recall, thresholds = precision_recall_curve(y_true=y_true, probas_pred=df.predictions.values)
disp = PrecisionRecallDisplay(precision=precision, recall=recall)
disp.plot()
plt.show()
precision_recall_thresholds = pd.DataFrame({
"precision": precision[:-1],
"recall": recall[:-1],
"thresholds": thresholds,
})precision_recall_thresholds| precision | recall | thresholds | |
|---|---|---|---|
| 0 | 0.048001 | 1.000000 | 0.000114 |
| 1 | 0.048001 | 1.000000 | 0.000116 |
| 2 | 0.048001 | 1.000000 | 0.000117 |
| 3 | 0.048001 | 1.000000 | 0.000125 |
| 4 | 0.048002 | 1.000000 | 0.000125 |
| ... | ... | ... | ... |
| 889313 | 1.000000 | 0.000110 | 0.992152 |
| 889314 | 1.000000 | 0.000088 | 0.992280 |
| 889315 | 1.000000 | 0.000066 | 0.992347 |
| 889316 | 1.000000 | 0.000044 | 0.992431 |
| 889317 | 1.000000 | 0.000022 | 0.992971 |
889318 rows × 3 columns
A threshold of approx. 0.2 seems to be a good compromise for precision/recall.
mhcpred_proba_threshold = 0.2df = pd.read_csv(str(output_path / "mhcpred_predictions.csv"))
y_true = df.hit.values
y_pred = df.predictions.values >= mhcpred_proba_threshold
acc = accuracy_score(y_true=df.hit.values, y_pred=y_pred)
confusion_mat = confusion_matrix(y_true=df.hit.values, y_pred=y_pred)
balanced_acc = balanced_accuracy_score(y_true=df.hit.values, y_pred=y_pred)
class_report = classification_report(y_true=y_true, y_pred=y_pred, output_dict=False)
disp = ConfusionMatrixDisplay(confusion_matrix=confusion_mat)
disp.plot()
plt.show()
acc, balanced_acc(0.9697360118629252, 0.8462451280426622)print(class_report) precision recall f1-score support
False 0.99 0.98 0.98 900725
True 0.68 0.71 0.69 45416
accuracy 0.97 946141
macro avg 0.83 0.85 0.84 946141
weighted avg 0.97 0.97 0.97 946141
Balanced accuracy improves from 0.778 to 0.846. Precision drops but recall rises, which is the right trade-off on a task where missing a presented peptide is the costlier error.
Source: MetricsConclusion
MHC molecules present peptide fragments on the cell surface, which is how T cells detect foreign or abnormal proteins and how immune responses against pathogens and cancer cells are mediated. Predicting which peptides a given allele presents is therefore worth doing computationally, and the sequence-to-function structure of the problem makes it a natural machine learning target.
Two approaches were compared on the same NetMHCpan-4.1 test set of 946,141 peptide-allele pairs, 4.8% of them positive.
| Model | Threshold | Accuracy | Balanced accuracy | Precision (hit) | Recall (hit) |
|---|---|---|---|---|---|
| MHCflurry, affinity model | 2% percentile | 0.973 | 0.922 | 0.67 | 0.86 |
| mhcpred | 0.5 | 0.973 | 0.778 | 0.82 | 0.56 |
| mhcpred | 0.2 | 0.970 | 0.846 | 0.68 | 0.71 |
The accuracy column shows why accuracy is the wrong metric here: all three rows sit at 0.97, which a model that never predicts a hit would nearly match given a 4.8% positive rate. Balanced accuracy separates them.
At its default threshold, mhcpred is precise but timid. A recall of 0.56 means it misses roughly 20,000 of the 45,416 presented peptides in the test set, which is what one expects from a network trained for two epochs on data that is 95% negative with no correction for the imbalance. Lowering the decision threshold to 0.2 trades precision for recall and lifts balanced accuracy from 0.778 to 0.846, still short of MHCflurry’s 0.922.
The gap is unsurprising. MHCflurry is an ensemble trained to convergence on a larger curated dataset, with random negative sampling and per-allele calibrated percentile ranks; mhcpred is a single network trained for two epochs — 524 seconds over 6.6 million training points, validation loss falling from 0.069 to 0.051 — with the class imbalance left untreated. The natural next steps are the ones this shortcut skipped: train to convergence with early stopping, handle the imbalance through class weights or negative sampling, evaluate per allele rather than in aggregate, since the four test-only alleles are the real test of the pseudo-sequence encoding, and benchmark against Class1PresentationPredictor rather than the affinity model.
One caveat on the comparison itself: the merge in mhcflurry_benchmark.py joins on allele and peptide, which duplicates rows wherever a pair occurs more than once, so the MHCflurry results are computed over 946,419 rows against mhcpred’s 946,141. The difference is too small to affect the conclusions, but the two models are not scored on exactly the same rows.
References
Antunes, Dinler A., Jayvee R. Abella, Didier Devaurs, Maurício M. Rigo, and Lydia E. Kavraki. 2018. “Structure-based Methods for Binding Mode and Binding Affinity Prediction for Peptide-MHC Complexes.” Current Topics in Medicinal Chemistry 18 (26): 2239–55. https://doi.org/10.2174/1568026619666181224101744.
O’Donnell, Timothy J, Alex Rubinsteyn, and Uri Laserson. 2020. “MHCflurry 2.0: Improved Pan-Allele Prediction of MHC Class i-Presented Peptides by Incorporating Antigen Processing.” Cell Systems 11 (1): 42–48.
Reynisson, Birkir, Bruno Alvarez, Sinu Paul, Bjoern Peters, and Morten Nielsen. 2020. “NetMHCpan-4.1 and NetMHCIIpan-4.0: improved predictions of MHC antigen presentation by concurrent motif deconvolution and integration of MS MHC eluted ligand data.” Nucleic Acids Research 48 (W1): W449–54. https://doi.org/10.1093/nar/gkaa379.
Citation
BibTeX citation:
@online{brosse2025,
author = {Brosse, Nicolas},
title = {Predicting {MHC-peptide} Presentation with Machine Learning},
date = {2025-02-06},
url = {https://nbrosse.github.io/posts/mhc/mhc.html},
langid = {en}
}
For attribution, please cite this work as:
Brosse, Nicolas. 2025. “Predicting MHC-Peptide Presentation with
Machine Learning.” February 6. https://nbrosse.github.io/posts/mhc/mhc.html.