Comparison of HLA ligand elution data and binding predictions reveals varying prediction performance for the multiple motifs recognized by HLA-DQ2.5.

Binding prediction tools are commonly used to identify peptides presented on MHC class II molecules. Recently, a wealth of data in the form of naturally eluted ligands has become available and discrepancies between ligand elution data and binding predictions have been reported. Quantitative metrics for such comparisons are currently lacking. In this study, we assessed how efficiently MHC class II binding predictions can identify naturally eluted peptides, and investigated instances with discrepancies between the two methods in detail. We found that, in general, MHC class II eluted ligands are predicted to bind to their reported restriction element with high affinity. But, for several studies reporting an increased number of ligands that were not predicted to bind, we found that the reported MHC restriction was ambiguous. Additional analyses determined that most of the ligands predicted to not bind, are predicted to bind other co-expressed MHC class II molecules. For selected alleles, we addressed discrepancies between elution data and binding predictions by experimental measurements and found that predicted and measured affinities correlate well. For DQA1*05:01/DQB1*02:01 (DQ2.5) however, binding predictions did miss several peptides that were determined experimentally to be binders. For these peptides and several known DQ2.5 binders, we determined key residues for conferring DQ2.5 binding capacity, which revealed that DQ2.5 utilizes two different binding motifs, of which only one is predicted effectively. These findings have important implications for the interpretation of ligand elution data and for the improvement of MHC class II binding predictions.

DQA1*05:01/DQB1*02:01

epstein‐barr virus

human leucocyte antigen

immune epitope database

equilibrium dissociation constant

monoclonal antibody

major histocompatibility complex

mass spectrometry

nanomolar

binding pockets 1‐9 in MHC binding groove

protein database

T cell receptor

INTRODUCTION

Major histocompatibility complex (MHC) class II molecules are expressed on professional antigen‐presenting cells and present antigenic peptides to CD4+ T cells, which play an important role in regulating immune responses. MHC class II, called human leucocyte antigen (HLA) class II in humans, and CD4+ T‐cell responses are of major significance for autoimmunity and antitumor immunity. It is known, for example, that coeliac disease is driven by CD4+ T‐cell responses to gluten proteins presented by certain HLA class II molecules. More recently, it has been shown that neoantigens presented on HLA class II elicit potent antitumor responses in patients and that CD4+ T cells are required for successful anticancer immunotherapy. , , Given these developments, there is an increased interest in developing efficient approaches for identifying epitopes presented by MHC class II. As experimental identification of epitopes remains time consuming and costly, and subject to several technical challenges, MHC class II binding prediction tools are now commonly used.

Most MHC binding prediction tools are based on machine learning algorithms trained on data sets of known peptide–MHC binding affinities. With an increasing amount of data available for training, prediction algorithms have been able to achieve high accuracy. It is however generally accepted that predicting binding to MHC class II is a more challenging problem when compared to predicting binding to MHC class I due to a number of reasons. MHC class II molecules consist of heterodimers of alpha and beta chains. These chains are encoded by genes in the HLA‐DP, HLA‐DQ or HLA‐DR region of chromosome 6 in humans and are highly polymorphic in the general population. Most of the polymorphism is associated with residues forming the MHC’s peptide‐binding region, which accounts in large part for the high degree of variation in binding specificity observed between different alleles. The presence of two potentially polymorphic chains makes the assignment of MHC restriction harder for class II, compared to class I molecules, which consist of a variable alpha chain and the monomorphic beta‐2‐microglobulin molecule.

The main energy of interaction for peptide binding to MHC class I and II is conferred by a peptide core of approximately 9 residues. However, in contrast to MHC class I, the binding groove of MHC class II is open on both ends, allowing peptides of variable length to bind. The presence of amino acids flanking the peptide core is also necessary for efficient binding to MHC class II. As a result, MHC class II has a preferred minimal length of about 13 residues for binding peptides. The presence of such over‐hanging residues, plus the fact that the open ends of the binding groove afford multiple different peptide frames, means that machine learning algorithms for MHC class II binding have to consider multiple possible peptide cores.

The specificity of an MHC molecule can be described as a ‘motif’. Crystallographic analyses of MHC–peptide complexes have shown that different MHC molecules have different pockets within their binding groove. When a peptide is bound, the amino acid side chains of so‐called anchor residues in the peptide specifically insert into these pockets within the MHC binding groove. Both the position and specificity of these anchor residues are similar for different peptides binding the same MHC molecule and together describe the peptide‐binding motif of a given MHC molecule. For MHC class I binding peptides, the anchor residues are located at or near the N‐ and C‐termini of the peptide, usually P2 and P9 in the case of 9‐mer peptides. For MHC class I binding peptides anchor residues can be at various distances from the ends of the peptide and at various positions of the binding core. For DRB1*01:01 and DRB1*15:01, for example, the prototypical anchor positions are at P1, P4, P6, P7 and P9. But other molecules, such as DRB1*11:01 and DQA1*03:01/DQB1*03:01, appear to have slightly different anchor spacing, with the second anchor at P3 instead of P4, or an anchor at P5. While for many MHC class I and class II molecules distinct binding motifs could be defined, there are also instances of MHC molecules that are suggested to carry multiple binding motifs, such as HLA‐DRB1*03:01.

Recent advances in mass spectroscopy (MS)‐based techniques have enabled the development of high‐throughput ligand elution assays, which allow the identification of thousands of natural ligands with a single experiment. , As eluted ligands undergo the natural processing and MHC loading cascade, the resulting elution data will inherently contain valuable biological information that is not available when only peptide binding is considered. Given these advantages, ligand elution experiments have become quite popular, leading to a wealth of data in the form of eluted ligands from both MHC class I and class II molecules.

As eluted ligands were isolated from MHC molecules that presented them on the cell surface, their capacity to bind MHC is inherently supported. With the increasing amount of ligand elution data, discrepancies between MHC class II ligand elution and binding predictions have been reported and the efficiency of binding predictions has been put into doubt. , , , However, quantitative metrics for such comparisons are currently lacking. In this study, we wanted to systematically investigate how well HLA class II binding predictions and ligand elution data agree with each other. To achieve that, we assembled HLA class II ligand elution data reported in literature and performed multiple comparative analyses between elution data and binding predictions: (i) we analysed how many of the eluted ligands were predicted to bind their reported restricting MHC; (ii) for several MHC class II alleles, a significant fraction of reported ligands were not predicted to bind and we investigated these in detail; (iii) we then chose four representative HLA class II alleles and specifically selected peptides for which ligand elution data and binding predictions disagreed and experimentally measured binding affinity to determine which method was correct; and (iv) we found that the correlation between predicted and measured binding is lowest for DQA1*05:01/DQB1*02:01 (DQ2.5). We therefore investigated this allele in detail and performed experiments to the determine key residues for binding. We found that DQ2.5 can bind peptides using alternate binding modes and that, currently, only one of these is being predicted accurately.

MATERIAL AND METHODS

Data assembly

The Immune Epitope Database (IEDB) catalogs experimental data on antibody and T‐cell epitopes as well as MHC ligands. The IEDB was queried to retrieve ligand elution data for MHC class II using the following criteria: ‘Positive Assays Only, Epitope Structure: Linear Sequence, No T cell assays, No B cell assays, MHC ligand assays: MHC ligand elution assay, MHC Restriction Type: Class II’. The collected data included ligand sequence, the MHC class II allele from which the ligand was eluted from, details of the source protein from which the ligands originated from and PubMed identifiers (PMIDs) of the studies that reported the ligand.

We further refined this set to only include peptides that were eluted from HLA molecules typed with 4‐digit resolution, resulting in a set of 24,601 peptides from 27 different HLA class II molecules. The length of the ligands in this set ranged from 3 to 43 residues; based on previous analyses, we only retained peptides of 15–20 amino acids in length. After this step, the data set contained 12,506 peptides from 27 different HLA class II molecules. We then further filtered the set by considering only alleles for which at least 10 ligands were reported. The resulting final data set used for analysis included 12,449 sequences, eluted from 16 different HLA class II molecules. This data set is provided as Table S1.

HLA binding predictions

We performed binding predictions using the standalone version of NetMHCIIpan (version 3.1), as this tool provides predictions for most HLA class II alleles. , HLA class II binding predictions were performed for all 12,449 sequences and their reported HLA class II restriction element. Peptides with predicted percentile ranks of >20% were considered predicted non‐binders. This threshold was selected based on previous studies analysing the binding affinity of HLA class II‐restricted T‐cell epitopes. , , ,

To analyse contamination of ligand elution data, we also performed binding predictions for all peptides to common DR and DQ alleles: HLA‐DRB1*01:01, HLA‐DRB1*03:01, HLA‐DRB1*04:01, HLA‐DRB1*04:05, HLA‐DRB1*07:01, HLA‐DRB1*08:02, HLA‐DRB1*09:01, HLA‐DRB1*11:01, HLA‐DRB1*12:01, HLA‐DRB1*13:02, HLA‐DRB1*15:01, HLA‐DRB3*01:01, HLA‐DRB3*02:02, HLA‐DRB4*01:01, HLA‐DRB5*01:01, HLA‐DQA1*05:01/DQB1*02:01, HLA‐DQA1*05:01/DQB1*03:01, HLA‐DQA1*03:01/DQB1*03:02, HLA‐DQA1*04:01/DQB1*04:02, HLA‐DQA1*01:01/DQB1*05:01, HLA‐DQA1*01:02/DQB1*06:02.

Selection of a peptide set with discrepancies

We selected 4 alleles to further investigate peptides with discrepancies between ligand elution data and binding predictions: HLA‐DQA1*05:01/DQB1*02:01, DRB1*04:01, DRB1*15:01 and DRB1*01:01. These representative alleles were selected considering the number of reported eluted ligands, frequency of the allele in the general human population and whether an in‐house quantitative binding assay was available. In total, 3,443 eluted peptides were available for the 4 selected HLA class II alleles, of which 1,611 had predicted %ranks >20. Accordingly, these were considered as ‘predicted non‐binder’. The composition of the data set after each filtering step is described in detail in Table 1.

Table 1

Composition of eluted ligand dataset after each filtering step and the number of peptides selected to be tested experimentally for the four alleles selected for detailed investigation

Allele	Peptides in ligand elution dataset downloaded from IEDB	Peptides with length 15–20	Peptides with predicted %rank>20	Predicted non‐binder redundant peptides removed	Predicted non‐binders selected for testing	Predicted binders selected for testing
HLA‐DQA1*05:01/DQB1*02:01	2245	2056	1097	811	30	20
HLA‐DRB1*01:01	638	594	114	97	20	15
HLA‐DRB1*04:01	555	463	222	176	30	20
HLA‐DRB1*15:01	388	330	178	132	17	18
4 selected HLA combined	3,826	3,443	1,611	1,216	97	73

Analysis of sequence similarity

Many of the eluted peptides are very similar to each other as they contain the same MHC class II binding core with varying flanking residues of different lengths. When selecting peptides for experimental testing, we wanted to avoid testing the same MHC class II binding core multiple times. Thus, to reduce redundancy, we removed sequences with >80% sequence similarity. We used the IEDB clustering tool, which groups the set of input peptides into clusters based on sequence identity. Here, a cluster is defined as a group of sequences that have a similarity greater than the specified minimum sequence identity threshold (80% in our case). We used all 1611 eluted peptides that were not predicted to bind for the 4 selected alleles as input, and from each cluster, we randomly selected one peptide. As a result, 1216 peptides that had <80% sequence similarity remained in the data set, from which peptides were selected for experimental testing.

HLA binding measurements

Purification of MHC molecules by affinity chromatography was performed as detailed elsewhere. Briefly, Epstein‐Barr virus (EBV) transformed homozygous cell lines or single MHC allele transfected RM3 or fibroblast lines were utilized as sources of HLA class II MHC molecules. MHC class II molecules were purified from cell pellet lysates by repeated passage over Protein A Sepharose beads, conjugated with the monoclonal antibodies (mAbs) LB3.1 (anti‐HLA‐DR), SPV‐L3 (anti‐HLA‐DQ) and B7/21 (anti‐HLA‐DP). Protein purity, concentration and the effectiveness of depletion steps were monitored by SDS‐PAGE and BCA assay.

Classical competition assays to quantitatively measure peptide binding to MHC class II molecules were performed as detailed elsewhere. These assays are based on inhibition of binding of high‐affinity radiolabeled peptides to purified MHC molecules. Briefly, 0·1–1 nM of radiolabeled peptide was co‐incubated at room temperature or 37°C with purified MHC in the presence of a cocktail of protease inhibitors. Following a two‐ to four‐day incubation, MHC bound radioactivity was determined by capturing MHC/peptide complexes on MHC locus specific mAb‐coated Lumitrac 600 plates (Greiner Bio‐one, Frickenhausen, Germany). We then measured bound counts per minute (cpm) using the TopCount (Packard Instrument Co., Meriden, CT, USA) microscintillation counter. The concentration of peptide yielding 50% inhibition of binding of the radiolabeled peptide was calculated. Under the conditions utilized, where the concentration of the radiolabeled peptide is less than the concentration of MHC, and the affinity of the peptide (IC50 nM) is greater than or equal to the concentration of MHC (or more formally, [label] < [MHC] and IC50 ≥ [MHC]) measured IC50 values are reasonable approximations of true equilibrium dissociation constants (Kd). , Each competitor peptide was tested at six different concentrations covering a 100,000‐fold range, and in three or more independent experiments. As a positive control, the unlabelled version of the radiolabeled probe was also tested in each experiment.

Transforming measured IC50 values to percentile ranks

To compare the measured affinities to the predicted percentile ranks, we converted the measured IC50 values to percentile ranks as well. The percentile rank in this context describes the rank of a given IC50 value in a larger reference set of measured IC50 values. We obtained various data sets from the IEDB that contained binding measurements of overlapping 15mers for the four HLA class II alleles of interest: , , , HLA‐DQA1*05:01/DQB1*02:01 (n = 537), DRB1*04:01 (n = 610), DRB1*15:01 (n = 537) and DRB1*01:01 (n = 1,277). Each measured IC50 value was then converted to a percentile rank by calculating its rank in the allele‐specific reference set.

RESULTS

HLA class II reported eluted ligands are, in general, predicted to bind to their restriction element

We first addressed whether ligands eluted from a specific HLA would also be predicted to bind that allele. We utilized a percentile rank threshold of 20% to classify peptides into binders and non‐binders. This threshold is a rather comprehensive; however, as we wanted to assess the binding capacity of reported eluted ligands, a comprehensive threshold was suitable. In a setting where the rationale is to use binding predictions to select the best binders for experimental testing, another strategy would be recommended, such as selecting the top 10% of predicted binders.

To address the extent to which HLA class II eluted peptides are predicted by the 20%‐rank threshold, we used the Immune Epitope Database (IEDB) and downloaded all peptides reported as eluted from HLA class II‐expressing cells. We retained peptides of 15–20 amino acids in length and eluted from HLA molecules typed to the 4‐digit resolution level. We only further considered alleles for which at least 10 ligands were reported. The resulting data set included 12,449 sequences eluted from 16 different HLA class II molecules.

Next, for each of these 16 HLA class II molecules and their associated ligands, we performed predictions using the NetMHCIIpan algorithm. , As was the case with T‐cell epitopes, significant variation in prediction performance across alleles was noted. On average, 57% ±17% SD (Table 2) of eluted peptides was predicted to bind the corresponding allele at the 20% rank threshold.

Table 2

Ligand elution data and associated prediction efficacy

Allele	ligands	%recovered at 20%
HLA‐DRB5*01:01	131	75%
HLA‐DRB3*02:02	10	80%
HLA‐DRB1*15:01	330	46%
HLA‐DRB1*13:02	17	35%
HLA‐DRB1*11:01	50	62%
HLA‐DRB1*04:05	1206	59%
HLA‐DRB1*04:04	35	69%
HLA‐DRB1*04:02	29	31%
HLA‐DRB1*04:01	463	52%
HLA‐DRB1*03:01	79	49%
HLA‐DRB1*01:01	594	81%
HLA‐DQA1*05:05/DQB1*03:01	3908	58%
HLA‐DQA1*05:01/DQB1*02:01	2055	47%
HLA‐DQA1*02:01/DQB1*02:02	3492	36%
HLA‐DQA1*01:02/DQB1*06:02	27	85%
HLA‐DQA1*01:01/DQB1*05:01	23	52%

Thus, while the majority of eluted peptides could be effectively predicted by existing algorithms, the overall prediction efficiency was lower than what has been reported for known T‐cell epitopes. Further, the data also revealed significant inter‐allele variability and identified several alleles for which less than 60% of reported eluted ligands were predicted to bind below 20% rank (Table 2).

Most non‐predicted ligands are predicted to bind other HLA molecules

We addressed whether the 556 DRB1*04:01, DRB1*03:01, DRB1*04:02 and DRB1*13:02 ligands that were reported in these studies could have been eluted from associated and co‐expressed DRB3/4/5 molecules, or other contaminant HLA class II. Out of these 556 reported ligands 277 (50%) were not predicted to bind the reported DRB1 at the 20% rank. Of the 277 peptides not predicted to bind the reported DRB1, 30 (10%) were predicted to bind other representative DRB3/4/5 molecules (i.e. DRB3*01:01, DRB3*02:02, DRB4*01:01 and DRB5*01:01). We additionally found that out of the 279 peptides predicted to bind the reported DRB1, 119 peptides (43%) were also predicted to bind other representative DRB3/4/5 molecules (Table S2).

Next, in light of the observation that several papers had isolated ligands from heterozygous DR samples, we examined whether the 247 ligands that were not predicted to bind, neither the reported DRB1 nor other representative DRB3/4/5 molecules, may be predicted to bind other common DRB1 or DQ alleles. We found that 60% (148/247) of the non‐predicted DRB1‐reported eluted ligands were, in fact, predicted to bind other common DR or DQ alleles. Conducting the same analysis for the 279 peptides that were predicted to the reported DRB1 molecule, we found that 99% (275/279) of them were also predicted to bind other common DRB1 or DQ alleles.

Taken together, of the 277 ligands that were not predicted to bind the reported HLA at the 20% rank threshold, 64% were predicted to bind other HLA molecules, suggesting that a wrong HLA might have been assigned to these ligands. Our results additionally suggest that ligands that were predicted to bind the reported DRB1 molecule could also be promiscuous binders with the capacity to bind multiple different HLA alleles.

Addressing discrepancies between elution data and binding predictions by experimental affinity determinations

Next, we chose representative HLA class II alleles from the set described above, considering the number of reported eluted ligands, frequency of the allele in the general human population and whether an in‐house quantitative binding assay was available. Specifically, we selected HLA‐DQA1*05:01/DQB1*02:01, DRB1*04:01 and DRB1*15:01 for further study, all alleles for which less than 55% of ligands were predicted to bind at the 20% rank threshold. As a control, we selected DRB1*01:01, an allele for which over 80% of ligands were predicted to bind with affinities over the 20% rank.

For these four selected alleles, 3,443 eluted peptides with lengths of 15–20 amino acids were available, of which 1611 (47%) were not predicted to bind at the 20% rank threshold (Table 1). We excluded sequences with >80% sequence similarity to reduce redundancy, resulting in a set of 1,216 peptides. From this set, we selected between 17 and 30 peptides each for HLA (Table 2). This resulted in a total of 97 eluted peptide sequences that were not predicted to bind at 20%‐rank. As a control, we selected between 15 and 20 eluted peptides per allele that were predicted to bind with %‐rank of 10 or better. Binding affinities for the 97 peptides predicted not to bind and 73 predicted binders were measured using standard competitive binding assays, as described in the Methods section.

To compare the measured affinities to the predicted percentile ranks, we converted the measured IC50 values to percentile ranks as well. The percentile rank in this context describes the rank of a given IC50 value in a larger reference set of measured IC50 values. We obtained various data sets from the IEDB that contained binding measurements of overlapping 15mers for the four HLA class II alleles of interest , , , and combined the data to be used as the reference data set of measured IC50 values. Each measured IC50 value was then converted to a percentile rank by calculating its rank in this allele‐specific reference set.

Table S3 lists all tested peptides together with predicted and measured IC50 and %‐rank values.

Overall (Table 3), 67 of the 73 predicted binders (92%) bound the corresponding HLA molecule when tested experimentally (shown as squares in Figure 1), consistent with previous data, and as summarized above. Similarly, 79 of the 97 eluted peptides (81%) predicted not to bind (Table 3) indeed did not bind when tested experimentally (shown as circles in Figure 1). Thus, these results confirmed the general efficacy of the binding predictions, as applied to the experimentally determined binding capacity of reported eluted ligands.

Table 3

Composition of peptide set and associated binding results

HLA allele	Pred. non‐binders selected	Pred. non‐binders that don't bind	Pred. binders selected	Pred. binders that bound
DQA1*05:01/DQB1*02:01	30	22 (73%)	20	18 (90%)
DRB1*01:01	20	17 (85%)	15	14 (93%)
DRB1*04:01	30	28 (93%)	20	17 (85%)
DRB1*15:01	17	12 (71%)	18	18 (100%)
Total	97	79 (81%)	73	67 (92%)

Figure 1

Measured percentile ranks are plotted against predicted ones. Colors indicate different alleles, squares predicted binders, and circles predicted non‐binders. A red outline depicts cases with large discrepancies between measured and predicted affinities, that is measured rank < 10% and predicted rank > 30%.

Importantly, for some alleles a relatively large fraction of eluted peptides did not, in fact, bind the reported corresponding allele when measured experimentally. In light of the complexity and limitations of eluted ligand determinations, it appears likely that at least some of these reported ligands have an incorrectly assigned HLA molecule, or are otherwise contaminants.

Correlation between predicted and measured binding is lowest for DQA1*05:01/DQB1*02:01

Comparison of the NetMHCIIpan predicted %‐ranks to the measured %‐ranks revealed good correlation overall (r = 0·8, Pearson's correlation). However, this analysis also identified some inter‐allele variability and we found that the correlation between predicted and measured %‐ranks is lowest for DQA1*05:01/DQB1*02:01 (r = 0·743, Pearson's correlation). We additionally observed several cases with substantial disagreement between predicted and measured binding. In particular, the analysis identified 8 peptides (outlined in red in Figure 1) that were not predicted to bind and had %‐ranks of >30% but measured %‐ranks were all <10%. This suggests that these peptides are likely real ligands whose binding capacity was missed by the predictions. Further inspection of the affected peptide‐HLA pairs showed that the majority of these peptides (5 of 8) were associated with DQA1*05:01/DQB1*02:01 (Table 4).

Table 4

Peptides showing greatest discrepancy between predicted and measured affinity

Sequence	HLA allele	Measured %‐Rank	Predicted %‐Rank
APDTRFFVPEPGGRGAAP	HLA‐DRB1*01:01	2.19	34.31
WISKQEYDESGPSIVHRK	HLA‐DRB1*15:01	7.26	79.28
SPTEPKNYGSYSTQA	HLA‐DRB1*15:01	2.61	86.56
GEPDYVNGEVAATEA	HLA‐DQA1*05:01/DQB1*02:01	2.98	30.02
GDSDLQLDRISVYYNEA	HLA‐DQA1*05:01/DQB1*02:01	6.7	31.25
RQEEPEYENVVPISRPP	HLA‐DQA1*05:01/DQB1*02:01	3.54	40.81
KPPTADLFTGVLPNGYNPP	HLA‐DQA1*05:01/DQB1*02:01	3.91	43.8
LPGRENYSSVDANGIQ	HLA‐DQA1*05:01/DQB1*02:01	2.23	55.87

Accordingly, we next sought to characterize the respective modes of binding of the 5 ‘discordant’ peptides to DQA1*05:01/DQB1*02:01 (hereafter DQ2.5). As controls, we also included 2 peptides eluted from DQ2.5 that were predicted to bind. We also included one peptide (aTIP 430) previously reported as a DQ2.5‐restricted T‐cell epitope and predicted to bind DQ2.5, and another (DQA 16) previously identified as an endogenously bound DQ2.5 ligand, but not predicted to bind DQ2.5. The DQ2.5 binding capacity of both aTIP 430 and DQA 16 was previously characterized using purified HLA , and confirmed in the present study. Table 4 shows the 9 peptides included in the single amino acid substitution experiments. Table 3 lists the 8 peptides with greatest discrepancy between predicted and measured affinity, where an enrichment for DQ2.5 was noted (5/8). Some peptides are listed in both Tables, as they are relevant in both contexts.

Determining key residues for DQ2.5 binding within the eluted ligands

For each of the ligands in Table 4, we experimentally determined the residues crucial for conferring DQ2.5 binding capacity. We used a previously described single amino acid substitution approach, scanning each peptide sequence by introducing non‐conservative substitutions at each position. Lysine (K) was selected as the basis for the analysis as it is non‐conservative with all amino acids except those with a positive charge (i.e. R, H and K). In cases where the native residue carried a positive charge, a glutamic acid (E) substitution was utilized. Each single substitution analog, and the original unmodified peptide, was tested for its capacity to bind DQ2.5. All peptides together with the measured IC50 values and ranks are detailed in Table S4.

All of the original unmodified peptides bound with IC50 values <1000 nM and %‐ranks <20%. However, several of the single substitution analogs did not bind at all, or bound with considerably lower affinity when compared to the corresponding unmodified peptide. The measured binding affinity of each analog was compared to the unmodified peptide, and a relative binding capacity was calculated as the fold change in IC50 between the peptide and its analog (Figure 2). Positions associated with a fourfold or greater reduction in binding capacity (highlighted by green shading in Figure 2) were defined as potential MHC contacts. These decreases in affinity allow mapping potential MHC contact residues for each ligand. The analysis revealed that on average 5·8 residues/ligand were identified as important for DQ2.5 binding, with a range from 2 to 9 residues. No clear, simple, common pattern was readily apparent across ligands.

Figure 2

The fold change in IC50 between the unmodified peptide and the corresponding analog with substitutions at each position scanning the length of the respective peptide is shown. The dashed line indicated a 10‐fold reduction in binding. Residues for which substitutions resulted in a greater than 10‐fold reduction in binding were defined as anchor residues and are highlighted in green here.

DQA1*05:01/DQB1*02:01 Motif 2 peptides are not effectively predicted

Previous studies indicated that DQ2.5 may be associated with at least two different binding motifs. One motif is characterized by a preference for aromatic or hydrophobic aliphatic residues in position 1 (P1), , , acidic residues in positions 4, 6 and/or 7 (P4, P6 and/or P7), and hydrophobic or aromatic residues at position 9 (P9). An independent study also identified a preference for acidic residues in P4 and P6, and large hydrophobic residues in P9, but with no specificity in P1. Collectively, we hereby refer to these preferences as ‘Motif 1’. Other reported patterns , associate proline or polar residues with P1, acidic or polar residues at various other positions, and a hydrophobic or polar residue in P9 (‘Motif 2’).

We divided all DQ2.5 ligands into three groups according to whether they fit Motif 1, Motif 2 or no motif following the patterns outlined above: Motif1 with aromatic or hydrophobic aliphatic residues (i.e. Phenylalanine, Tryptophan, Tyrosine, Leucine, Isoleucine, Valine or Methionine) in positions P1 and/or P9, and acidic residues (i.e. Aspartate or Glutamate) in positions P4, P6 and/or P7; Motif 2 with Proline or polar residues (i.e. Glutamine, Asparagine, Serine, Threonine, Lysine or Arginine) in P1, hydrophobic or polar residue (i.e. Glutamine, Asparagine, Aspartate, Glutamate, Phenylalanine, Tryptophan, Tyrosine, Leucine, Isoleucine, Valine, Methionine) in P9, acidic residues (i.e. Aspartate or Glutamate) in positions P6 and/or P7 and/or acidic or polar residues (i.e. Aspartate, Glutamate, Glutamine, Asparagine) in position P4. Using the peptides fitting Motif 1, Motif 2 or no motif, we generated three different sequence logos using the Seq2Logo tool (Figure 3A).

Figure 3

(A) All DQ2.5 ligands were divided into three groups according to whether they fit Motif 1, Motif 2, or no motif and sequence logos were generated. (B) The analyzed 9 DQ2.5 binding peptides are grouped according to the motif they conform to. Residues that were identified as anchor residues in single amino acid substitution experiments are highlighted in green. The red boxes highlight residues that explicitly match the position and specificity of Motif 1 or Motif 2. A + in the column prediction indicates that the peptide was predicted to bind DQ2.5 with a percentile rank of < 20%.

Upon closer inspection of the data shown in Figure 2, we were able to classify five ligands as fitting either Motif 1 or Motif 2 (Figure 3B). Specifically, three peptides appear to fit Motif 1, to include PI‐a 70 with Isoleucine at peptide position 8 (I8) as its P1 anchor, and Aspartate at position 11 (D11) and Glutamate at position 13 (E13) as P4 and P7 anchors, respectively. Motif 1 also describes the binding of the TMA 266 ligand, with Glutamate at peptide position 8 (E8) and Leucine at position 13 (L13) as the P4 and P9 anchors, and aTIP 430, with Valine at position 3 (V3), Aspartate at position 9 (D9) and Leucine at position 11 (L11) as P1, P7 and P9 anchors, respectively.

Two other peptides appear to fit Motif 2 (Figure 3B). These are the DQA 16 peptide, for which the single substitution analysis indicated the Proline in position 6 of the peptide (P6) as the P1 anchor and Phenylalanine at position 14 (F14) as the P9 anchor, and the SR 812 ligand where the Proline in position 5 (P5) is the P1 anchor and Isoleucine at position 13 (I13) is the P9 anchor.

The four remaining peptides did not fit either motif. These include the LAM 246 and Cochlin 160 peptides, which have aromatic/hydrophobic residues that fit as P1 and P9 anchors, but do not have acidic residues in the P4, P6 or P7 positions as described by Motif 1; meanwhile, aromatic residues in P1 render the LAM 246 and Cochlin 160 peptides incompatible with Motif 2. For the remaining two peptides (LAT2 21 and Tubulin 38), the single amino acid substitution scan also revealed binding modes not compatible with either Motif 1 or Motif 2.

Interestingly, when we examined the motif groups with respect to predicted binding, it was found that all three peptides fitting Motif 1 were also predicted by the NetMHCIIpan algorithm to bind DQ2.5 with %‐ranks of 20 or better (Table 5). Conversely, none of the peptides associated with Motif 2 or with no motif were predicted by the NetMHCIIpan algorithm. In conclusion, while 3/3 Motif 1‐positive peptides were predicted, 0/6 peptides that did not carry this motif were predicted (P = 0·019, Fisher's exact test).

Table 5

Peptides selected for DQ2.5 single amino acid substitution scan to determine the residues crucial for DQ2.5 binding

Protein	Pos	Sequence	Set	Measured %‐Rank	Predcited %‐Rank
PI‐a	70	NPTVFFDIAVDGEPL	Eluted	0.2	0.8
TMA	266	LAKTAFDEAIAELDT	Eluted	0.9	0.5
aTIP	430	EEVDMTPADALDDFD	Epitope	1.5	1.0
DQA	16	YQSYGPSGQYTHEFD	Eluted	10.6	55.0
SR	812	RQEEPEYENVVPISRPP	Eluted	2.2	40.8
LAM	246	KPPTADLFTGVLPNGYNPP	Eluted	5.6	43.8
Cochlin	160	LPGRENYSSVDANGIQ	Eluted	14.7	55.9
LAT2	21	GEPDYVNGEVAATEA	Eluted	2.1	30.0
Tubulin	38	GDSDLQLDRISVYYNEA	Eluted	6.9	31.3

To validate that the differential predictability of DQ2.5 ligands was associated with the presence of Motif 1, we assembled all described DQ2.5 eluted ligands and evaluated whether those predicted to bind DQ2.5 were enriched for the presence of Motif 1, compared to ligands that were not predicted to bind. Indeed, 20% (196/958) of the ligands predicted to bind did carry Motif 1. By contrast, only 5% (47/958) of ligands predicted to bind carried Motif 2 (P < 0·0001, chi‐square test). The remaining 75% of the ligands predicted to bind DQ2.5 did not conform to any known motif. The neural network architecture utilized by NetMHCIIpan theoretically has the capacity to learn and predict multiple binding motifs for an allele. Our results, however, suggest that the algorithm in practice does not efficiently predict Motif 2.

Taken together, these data highlight that DQ2.5 can bind peptides promiscuously using alternate modes. Further analysis of more peptides is necessary to fully understand what defines the way peptides bind or do not bind this molecule and to improve binding predictions.

DISCUSSION

In this study, we assessed how well MHC class II binding predictions and ligand elution data agree with each other. We found that the majority of experimentally identified MHC class II eluted ligands were predicted to bind to their reported restriction element. For HLA‐DR ligands that showed discrepancies, most of the ligands that were not predicted to bind were actually predicted to bind just below the threshold utilized, or were predicted to bind to other co‐expressed MHC class II molecules.

In contrast to the DR molecules studied, we found that for HLA‐DQA1*05:01/DQB1*02:01 (DQ 2.5) binding predictions did miss several peptides that were determined experimentally to be binders. By determining key residues for DQ2.5 binding in several peptides, we could trace discrepancies back to the existence of two binding motifs, only one of which was effectively predicted.

Our finding that the majority of eluted ligands were predicted to bind their restricting allele supports the use of these prediction tools that are commonly used. , This use had been put into doubt by reports that the predicted binding affinity of peptides to MHC class II correlates poorly with MHC antigen presentation , and peptide immunogenicity. , We have reported previously that in the case of HLA class II‐restricted T‐cell epitopes binding predictions do perform well, as we found that 83·3% of epitopes are predicted to bind their corresponding restricting allele with an affinity of 1000 nM or better.

In the present study, we investigated how efficiently binding predictions can predict eluted ligands, and found that, on average, 57% of eluted peptides are predicted to bind their corresponding allele below the 20%‐rank threshold. These data prove that the majority of eluted ligands can also be effectively predicted, although the percentage is lower than it is in the case for HLA class II‐restricted T‐cell epitopes. It was already shown in the case of HLA class I‐restricted T‐cell epitopes that alleles can vary in terms of the affinity threshold associated with immunogenicity. This may indeed be the case as well for HLA class II ligands and T‐cell epitopes, and needs to be further investigated.

When we scrutinized experimental procedures of the ligand elution studies, we found that the assigned restrictions for several MHC class II ligands are likely incorrect. This reflects experimental challenges of ligand elution assays. The first steps in a typical MHC class II ligand elution assay are lysing the antigen‐presenting cells, purifying the MHC molecules via immunoprecipitation and chromatographically eluting the MHC ligands. Here, the main drawback is that natural ligandomes are multi‐allelic. Different MHC molecules are expressed in the interrogated cell and the eluted ligands can correspond to any one of those MHC molecules. This makes it difficult to unambiguously annotate the ligands to the MHC molecule they were presumed to be eluted from.

Of the eleven DR molecules we studied, 5 (DRB1*04:01, DRB1*03:01, DRB1*04:02, DRB1*13:02 and DRB1*15:01) were associated with lower prediction efficiency. We suspected that the poor‐performing DR alleles may have resulted from an imprecise HLA assignment of the eluted ligands. To further explore this issue, we scrutinized papers that reported the eluted ligands associated with those DR molecules. HLA‐DR loci encompass a monomorphic alpha chain (DRA), which pairs with beta chains expressed by two different loci: DRB1 and DRB3/4/5. These two loci are in strong linkage disequilibrium and are co‐expressed in vivo. The only exceptions to this are the DR1, DR8 and DR10 haplotypes, for which there is no associated DRB3/4/5 locus. Thus, in heterozygous individuals, up to four different DR molecules (two with DRB1‐ and two with DRB3/4/5 beta chains) may be expressed, and even in homozygous humans, two different DR molecules may be co‐expressed. Importantly, all anti‐DR antibodies (L243, LB3.1 and others) commonly used to purify DR molecules and isolate eluted ligands are directed against the DR alpha chain, and thereby do not distinguish between DRB1 from DRB3/4/5 molecules. Thus, unless HLA mono‐allelic lines are used, both DRB1 and the DRB3/4/5 gene products will be co‐purified, and the pool of eluted ligands may contain a mixture of two different HLA class II molecules. All studies pertaining to those DR alleles associated with ‘lower efficiency of binding predictions’ suffered from the aforementioned issues. This included studies using the L243 or L227 monoclonal antibodies, both anti‐DRA, with EBV transformed homozygous , , , , or heterozygous , , cell lines. Other studies , similarly purified DR molecules from heterozygous patient material, yielding a mixture of DRB1 and DRB3 or DRB4 eluted peptides, as well as other undefined DR molecules expressed because of patient heterozygosity. While most authors do not mention issues pertaining to co‐expression, others acknowledge the problem. Koning et al. state that ‘affinity purification… resulted in the isolation of a mixture of the HLA‐DRB1 and the HLA‐DRB4‐locus product DR53’. Similarly, Hill et al. recognize that DRB1 and DRB4 molecules are tightly linked, but conclude that DRB4 is likely to represent a minor fraction of DR molecules, and is therefore ‘unlikely to influence results’.

This issue can be addressed, for example, by using genetically engineered mono‐allelic cell lines. This approach, however, is not feasible when patient samples are used, which will almost always be multi‐allelic. To overcome this, there are currently significant efforts to develop computational methods that can accurately annotate multi‐allelic ligands to their respective MHC. , , However, such methods will still rely on the identification of shared motifs of the eluted peptides.

Once ligands are eluted from MHC molecules, they are sequenced by tandem mass spectrometry (MS/MS). A key step is the identification of peptides by matching acquired MS/MS spectra against a customized protein sequence database. This step has an inherent trade‐off between sensitivity and specificity, which leads to a certain number of false spectrum identifications. , It is thus expected that some ligands identified by these methodologies are not correct. Given that we explicitly focused on peptides that were not predicted to bind, we might have enriched for such peptides that fall into the known false discovery rate of peptide elution experiments.

On the other hand, specifically focusing on these outliers might have enabled us to discover patterns that would have otherwise been missed. When we compared predicted and measured affinities, we identified several peptides that were not predicted to bind their associated allele but were experimentally determined to actually bind. This, therefore, indicates that these peptides are likely real ligands, but were missed by the predictions. As 5 out of 8 of these instances were associated with HLA‐DQA1*05:01/DQB1*02:01 (DQ2.5), we performed a detailed investigation of binding modes to DQ2.5. We performed single amino acid substitutions for these 5 discordant DQ2.5 associated peptides and an additional 4 peptides that are known DQ2.5 binders. Two distinct binding motifs have been described for this allele in the literature, and 5 out of the 9 analysed peptides were determined to fit either of these two motifs, while 4 peptides did not fit any motif.

It might be that these peptides that did not fit any known motif bind DQ2.5 in a mode engaging pockets in the binding groove in a not readily apparent alternate mode. But, as suggested by several elegant thermodynamic studies, , , , it is also possible that the primary energy of binding in these cases is driven by MHC backbone‐peptide side chain interactions. However, without additional data, this remains speculation in the present case.

Crystallography studies of DQ2.5 with bound peptides conforming to Motif 1, Motif 2 and no motif would enable the detailed characterization of the interactions in the different binding modes. Nine crystal structures of DQ2.5 are available in the Protein Database (PDB). However, as DQ2.5 is mainly studied in the context of coeliac disease, in 7 of these studies various deamidated gliadin peptides are reported, which all conform to Motif 2. , , The remaining 2 crystal structures report CLIP1 and CLIP2 peptides, which do not conform to any motif.

Crystal structures of DQ2.5 with bound gliadin peptides or mimics of these peptides in complex with T‐cell receptors (TCRs) showed that the TCR‐peptide contacts are at positions P2, P5, P6 and P8 of the MHC binding core. , Both Motif 1 and Motif 2 do not identify a preference for specific amino acids at positions P2, P5 and P8, and have the same preference for acidic residues at positions P6 and P7. It is hence theoretically possible that peptides conforming to Motif 1 or Motif 2 could be recognized as epitopes. Again, crystal structures of DQ2.5 with bound peptides conforming to Motif 1, ideally in complex with a TCR, would help to clarify the specific interactions of the different binding motifs.

Interestingly, peptides conforming to Motif 2 or no motif were not predicted correctly, indicating that prediction performance might depend on the motif of the peptide in question. We hypothesized that this might be due to a skewed Motif 1 to Motif 2 ratio in the DQ2.5 data set used to train NetMHCIIpan. We investigated the number of peptides in the training data set conforming to each motif and found that the majority of peptides did not fit any motif (53%) and there were less motif 1 peptides than motif 2 peptides (20% and 28%, respectively). The distribution of the different motifs in the training data set does hence not explain prediction efficacy of the different motifs.

While the neural network architecture utilized here hypothetically has the capacity to learn that multiple binding motifs are possible for a single allele, the algorithm in practice does not seem to pick this up. We used the most recent version of NetMHCIIpan (version 4.0) to investigate whether this issue persist in the updated implementation. The NetMHCIIpan 4.0 eluted ligand likelihood prediction (Rank_EL) was trained using eluted ligand data. As that includes all peptides we analysed in this study, the predictions will be overfitted and are not helpful in this context. The NetMHCIIpan 4.0 binding affinity prediction (Rank_BA) in contrast was trained using an extended data set of peptide‐MHC binding affinity data. Considering these updated predictions, all of the three tested Motif 1 peptides were predicted good binders with %‐ranks <1%. The predicted %‐ranks for two Motif 2 peptides in contrast are much weaker: while one peptide could be considered a weak binder with a %‐rank of 17%, the second peptide with a %‐rank of 48% was clearly not predicted to bind. While NetMHCIIpan 3.1 did not predict any of the 4 peptides that did not conform with any motif as binding, NetMHCIIpan 4.0 predicted two of them as binding with %‐ranks <20 (Table S5). While this analysis confirms our results regarding Motif 1, it also indicates that the updated neural network also has issues in predicting half of the peptides that do not conform with Motif 1.

There is already an example in the literature of an allele having two alternative motifs. Interestingly, this is HLA‐DRB1*03:01, which was also found in our study to be associated with lower prediction performance, with only 49% of peptides binding below 20%‐rank. It is possible that alternative binding motifs are a general problem with HLA class II molecules and that currently only peptides conforming to one of the associated motifs are predicted well. Conversely, peptides with other motifs are missed by predictions. This needs to be investigated in more detail by extending the analysis to additional alleles and peptides. The resulting data set from our study will be valuable for training machine learning methods to specifically detect cases with multiple binding motifs and improve prediction performance.

Our study utilizes experimentally eluted MHC ligand data as a set of ‘positive data points’, and examines which of these ligands is not identified with prediction methods, thus focusing on ‘false‐negative predictions’. Examining the flip‐side, namely ‘false‐positive predictions’, requires ligand elution data sets that exhaustively search for all possible MHC ligands in a set of antigens so that a set of ‘negative data points’ that with certainty are not presented by MHC can be identified. The data sets used in this study do not meet these criteria, and generating them in sufficiently large scale is challenging. For MHC class I, a comprehensive data set of eluted ligands derived from vaccinia virus infected murine cells that were further tested for T‐cell recognition has been published, revealing 83 major epitopes. Benchmarking revealed that the best prediction methods were able to capture more than 50% of these epitopes in the top N = 277 predictions within the set of N = 767,788 possible peptides. This means the false‐positive rate, defined as the fraction of negatives that is incorrectly predicted to be a positive is only 0·03%. While this seems close to perfect, when evaluating this performance based on precision, defined as the fraction of positive predictions made that is actually positive, is only 15%. This highlights that metrics of prediction performance have to be considered in the context of the problem considered, and identifying MHC ligands and/or T‐cell epitopes in a complex antigen is a ‘needle‐in‐the‐haystack’ problem. It also has to be pointed out that the prediction methods were asked to identify ligands based on the peptide sequence alone, and other factors that can lead to lack of presentation of a ligand, such as low expression level of the source protein, were not included in this comparison. Incorporating such additional factors should reduce the false‐positive rate. Similar data sets and benchmark comparisons should also to be generated for MHC class II.

CONFLICT OF INTEREST

All authors declare that they have no conflict of interest.

Table S1‐S5. Experimentally measured DQ2.5 peptides with NetMHCIIpan v4 predictions added.

Click here for additional data file.