Improved methods for predicting peptide binding affinity to MHC class II molecules

Data sets

The data set used to generate the new versions of NetMHCII and NetMHCIIpan contains peptide–MHC II binding affinities retrieved from the IEDB (http://www.iedb.org) in 2016. All data points are experimental IC₅₀ binding values, which were log‐transformed to fall in the range between 0 and 1 using the relation 1−log(IC₅₀ nm)/log(50 000) as explained by Nielsen et al.27. The 2016 data set contains 134 281 data points, covering 36 HLA‐DR, 27 HLA‐DQ, 9 HLA‐DP and 8 H‐2 molecules. The data set was split into five groups by clustering the common motif of peptides as described by Nielsen et al.28 and these five groups were used for a five‐fold cross‐validation. This 2016 data set is publicly available at http://www.cbs.dtu.dk/suppl/immunology/NetMHCIIpan-3.2. The data set used to develop the previous versions of NetMHCII and NetMHCIIpan is available at http://www.cbs.dtu.dk/suppl/immunology/NetMHCIIpan-3.0.

A summary of the data included in the 2013 and 2016 data sets is shown in Table 1 and a description of the full 2016 data set is available in the Supplementary material (Table S1).

Network training

The NetMHCII method was implemented as described by Nielsen and Lund15 and the NetMHCIIpan method was implemented as described by Andreatta et al.16 NetMHCII is an allele‐specific method that contains a specific predictor for each MHC molecule in the data set and it can therefore only predict binding affinities for MHC molecules found in the training data, whereas NetMHCIIpan is a pan‐specific method that can make predictions for any MHC molecule with a known protein sequence. To achieve its pan‐specificity, NetMHCIIpan incorporates information about the MHC‐II molecule, using a pseudo sequence consisting of residues that are considered important for peptide binding. This pseudo sequence is constructed using the method described by Karosiene et al.11 and is composed of 34 residues: 15 from the α‐chain and 19 from the β‐chain. Both methods were trained using a five‐fold cross‐validation set‐up. For each fold, we generate a network ensemble of individual networks trained without early stopping for 500 cycles with 10, 15, 40 and 60 hidden neurons using 10 different initial configurations, generating a total of 40 networks. This was done for each of the five training/test set combinations leading to a total of 200 networks. The peptide and the MHC pseudo sequence were encoded using the BLOSUM50 matrix and the PFR was encoded using the average BLOSUM scores on a maximum window of three amino acids at either end of the binding core.29 For each peptide core, the input to the neural network therefore consisted of the peptide core (9 × 20 = 180 inputs), the PFRs (2 × 20 = 40 inputs), the peptide length (2 inputs), the length of the C‐terminal and N‐terminal PFRs (2 × 2 = 4 inputs), resulting in a total of 226 input values for NetMHCII and 906 for NetMHCIIpan (an additional 34 × 20 = 680 input values from the pseudo sequence).

Binding core predictions

To improve the binding core predictions, we include the offset correction step to both NetMHCII and NetMHCIIpan. We followed the procedure described by Andreatta et al.16 and we evaluated the performance of this offset correction using the benchmark data set of 51 crystal structures of peptide–MHC‐II complexes.

Performance measures

The predictive performance of the different methods was measured using the area under the receiver operating characteristics curve (AUC). To classify peptides into binders and non‐binders, a binding threshold of 500 nm was used, classifying all peptides with an IC₅₀ binding value < 500 nm as binders. All performance values shown in this paper are averages of the AUC performance per MHC molecule using only molecules with more than 20 peptides and at least four binders.

Leave‐one‐molecule‐out network training

To assess the predictive performance of NetMHCIIpan in the situation where a molecule is not part of the training data, a leave‐one‐molecule‐out (LOMO) approach was applied.

To estimate LOMO performance for MHC molecule X, the NetMHCIIpan networks were trained using the fivefold cross‐validation set‐up from above. In the LOMO cross‐validation set‐up all binding data from molecule X were removed from the training sets and all test sets only include binding data from molecule X. This set‐up ensures that the method is trained without peptides binding to molecule X and it can therefore be used to evaluate the ability of the method to predict peptide binding of uncharacterized MHC‐II molecules.

Nearest neighbour distance calculation

The nearest neighbour distance is estimated from the alignment score of the HLA pseudo sequences using the relation $d=(s(\text{A,B}))/\left(\sqrt{s(\text{A,A})·s(\text{B,B})}\right)$. In this equation s(A,B) is the BLOSUM50 alignment score between the pseudo sequences for MHC molecules A and B, respectively.29 Nearest neighbours are found from the subset of molecules characterized with at least 50 data points and at least 10 binders.

Sequence logos

Sequence logos were constructed from the predicted binding cores of the top 1% strongest predicted binders using 200 000 natural random 15‐mer peptides and was visualized using seq2logo 30 with default settings.

Generation of HLA‐II distance trees

The HLA‐II distance tree was generated for each of the HLA‐DR, ‐DQ and ‐DP molecules in our data set using MHCcluster.31 To make the tree we first predicted the binding affinity for 200 000 natural random 15‐mer peptides using the new version of NetMHCIIpan. We then used MHCcluster to find the functional similarity between any two MHC molecules. MHCcluster calculates the similarity between two MHC molecules by correlating the union of the predicted top 10% strongest binding peptides. Using the bootstrap method in MHCcluster we generated 100 distance matrices and converted these to distance trees using the unweighted pair group method with arithmetic mean clustering. These trees were then combined into a consensus tree and visualized in splitstree.32 Sequence logos were constructed as explained above.

T‐cell epitope benchmark

A set of MHC‐II restricted T‐cell epitopes identified by multimer/tetramer staining assays was downloaded from IEDB. Only fully typed restrictions were included; that is, fully typed α‐ and β‐chains for HLA‐DQ and HLA‐DP, and a fully typed β‐chain for HLA‐DR (where the α‐chain is invariant). Epitopes with non‐natural amino acids were excluded. Also, epitopes with identical match to the peptides in the training data were excluded. The source protein sequence for each epitope was identified by mapping the annotated IEDB protein ID to the NCBI protein database. The final validation data set consisted of 1698 epitopes, restricted to 33 distinct MHC‐II molecules. For performance evaluation, the epitope source protein was split into overlapping peptides of the length of the epitope, and AUC and Frank values were calculated for each epitope–MHC pair annotating the epitopes as positive and all others as negatives. Here, Frank is the ratio of the number of peptides with a prediction score higher than the positive peptide to the number of peptides contained within the source protein. Hence, the Frank value is 0 if the positive peptide has the highest prediction value of all peptides within the source protein and a value of 0·5 in cases in which an equal number of peptides has a higher and lower prediction value compared with the positive peptide.

Comparing NetMHCII and NetMHCIIpan on a shared evaluation set

Using the data set from 2016, we retrained NetMHCII15 and NetMHCIIpan11 using a five‐fold cross‐validation setup to generate two new versions of these methods, named NetMHCII‐2.3 and NetMHCIIpan‐3.2. We then investigated how these new versions performed compared with the previous versions, which are NetMHCII‐2.2 and NetMHCIIpan‐3.1, trained on the 2013 data set. To make the comparison, we used the same fivefold cross‐validation set‐up and compared peptide data points in common between the 2013 and 2016 data sets. The result from this analysis in shown in Table 2.

The new versions of NetMHCII and NetMHCIIpan improved performance compared with the older versions (Table 2); but the performance gain was not statistically significant (P > 0·1 in both cases). Another interesting point is that the allele‐specific NetMHCII‐2.3 obtained a higher average performance than the pan‐specific NetMHCIIpan‐3.2, but this effect will be discussed later.

Performance of NetMHCIIpan on new data points for common MHC molecules

Using the five‐fold cross‐validation setup, we then evaluated the performance of the two versions of NetMHCII and NetMHCIIpan using only the subset of new peptides for the MHC molecules common between the old and new data sets. The result of this analysis is shown in Table 3 and it demonstrates a significant gain in predictive performance of the new versions (NetMHCII, P < 0.001 and NetMHCIIpan, P < 0.0003, using paired t‐test). This result underlines the importance of expanding the size of the training data even for previously characterized MHC molecules. [Correction added on 02 April 2018, after first online publication: In the preceding sentence, P < 0·005 and P < 0·001 was corrected to P < 0·001 and P < 0·0003 respectively.]

Binding core predictions

We evaluated the accuracy for binding core identification of the two updated MHC‐II binding prediction methods on the data set of peptide–MHC crystal structures described by Andreatta et al.16 Overall we find that (i) the inclusion of the offset correction has a substantial impact on the accuracy of binding core identification for both methods, and (ii) the overall accuracy of both methods is improved compared with the earlier version. For details see the Supplementary material (Table S2).

Performance of a consensus method

For predicting binding affinities to MHC‐I, it has been shown that a simple combination of the predictions from NetMHC27 and NetMHCpan10 gives a higher performance than using each method individually.33 We therefore made a similar combination of the predictions from NetMHCII‐2.3 and NetMHCIIpan‐3.2 to investigate if the performance could be improved for MHC‐II using this consensus approach. In the consensus method, we use an average of the prediction scores (values between 0 and 1) from NetMHCII‐2.3 and NetMHCIIpan‐3.2 to define the consensus method. The result of this analysis is shown in Fig. 1 and detailed performance values are found in the Supplementary material (Table S3). Figure 1(a) shows that the combination of NetMHCII‐2.3 and NetMCHIIpan‐3.2 has a significantly improved performance compared with each individual method and Fig. 1(b) shows that NetMHCIIpan‐3.2 outperforms NetMHCII‐2.3, especially for MHC molecules where only a few peptides are found in the data set.

Performance of NetMHCIIpan for previously uncharacterized MHC molecules

For NetMHCIIpan, we also tested the performance on MHC molecules that were not part of the 2013 data set (see Table 4). As expected, we observed that the new version of NetMHCIIpan had a significant increase in the predictive performance when compared with the previous version of NetMHCIIpan (P = 3·6 × 10⁻⁵, using a paired t‐test); this result therefore demonstrates the importance of expanding the allotypic coverage of the training data.

Leave‐one‐molecule‐out performance

The pan‐specific method is capable of making predictions for uncharacterized MHC molecules, so to assess the predictive performance of the NetMHCIIpan method in these situations we conducted a LOMO experiment. In the LOMO, the binding data for the MHC molecule in question were excluded from training and the resulting model was then evaluated using only binding data for the MHC molecule in question (for details see the Materials and methods). The LOMO experiment was made for all MHC molecules shared between the 2013 and the 2016 data sets, and the performance was evaluated on peptides shared between the two data sets. The result of this LOMO benchmark is shown in Table 5, together with the pseudo distances of the MHC molecule to each of the two training data sets estimated from the nearest neighbour sequence similarity as described in Materials and methods.

Table 5 shows an increased performance for NetMHCIIpan‐3.2‐LOMO compared with netMHCIIpan‐3.1‐LOMO. This gain is in general most pronounced for the MHC molecules that share a decrease in the pseudo sequence distance.

To further investigate this last observation, the LOMO performance evaluation was extended to include all MHC molecules in the 2016 data set. The result from this analysis is shown in Fig. 2 with a scatterplot of the relationship between the distance to the nearest neighbour in the training data set and the LOMO performance. The complete data used to create Fig. 2 can be found in Table S4. The figure shows that the HLA‐DQ and the HLA‐DP molecules have close nearest neighbours whereas the HLA‐DR and H‐2 molecules tend to have more distant neighbours. This figure also demonstrates a weak but statistically significant (P = 0·04 with exact permutation test) correlation between the LOMO performance and the distance to the nearest neighbour in the training data. This is in agreement with earlier findings for both MHC‐I and MHC‐II molecules10, 11 and shows how the predictive performance of the pan‐specific method depends on the distance to the nearest neighbour.

Distance tree for HLA molecules

Having arrived at the final retrained versions of NetMHCIIpan, we next use the MHccluster method31 to evaluate the similarities of binding motifs between the HLA molecules included in the 2016 training data. In short, the MHCcluster method estimates the similarity between two MHC molecules using the correlation between predicted binding values for a large set of random natural peptides. The similarity is 1 if the two molecules have a perfect binding specificity overlap and −1 if the two molecules share no specificity overlap (for details see Materials and methods). Comparing the binding pattern similarity between any two HLA class II molecules in the 2016 training data, we constructed the distance tree shown in Fig. 3. This figure confirms the earlier findings by Karosiene et al.:11 (i) the different loci show limited overlap in binding preference, (ii) HLA‐DP is less diverse compared with HLA‐DQ and HLA‐DR, and (iii) the diversity of HLA‐DQ can largely be split into three groups; one with preference for negatively charged amino acids towards the C‐terminus, one with a preference for positively charged amino acids towards the C‐terminus, and one with a preference for small amino acids at the anchor positions.

T‐cell epitope benchmark

We next evaluated the predictive performance of the two NetMHCIIpan methods on an IEDB T‐cell epitope data set. We queried the IEDB for MHC‐II‐restricted epitopes identified by tetramer/multi‐mer staining, which is the reference standard for epitope identification with known MHC restriction. For each epitope–MHC‐II pair, we calculated AUC and Frank values for the two NetMHCIIpan methods by predicting binding affinities to the MHC‐II restriction element of the epitope for all overlapping peptides with the same length as the epitope in the source protein sequence, annotating the epitope as positive and the remaining peptides as negative. This annotation is very stringent because peptides that share the same ligand binding‐core are counted as negatives even though they could be presented by the human MHC molecule; the set‐up will therefore most likely underestimate the predictive performance. The details from this analysis are found in Table S5 and the results are summarized in Fig. 4.

The Frank value is 0 if the positive peptide has the highest prediction value of all peptides within the source protein, and a value of 0·5 in cases where an equal number of peptides has a higher and lower prediction value compared with the positive peptide. Figure 4(a) shows that the Frank score for NetMHCIIpan‐3.1 is significantly lower than NetMHCIIpan‐3.1. It further shows that NetMHCIIpan‐3.2 has a median < 0·2 indicating that the positive peptide was found among the top 20% of the peptides from the source protein if sorted on their predicted peptide binding affinity. Figure 4(b) demonstrates a significant improvement in the AUC performance of NetMHCIIpan‐3.2 compared with NetMHCIIpan‐3.1. We speculate that the gain in predictive performance of NetMHCIIpan‐3.2 could be attributed to at least two factors, the inclusion of binding data for additional MHC‐II molecules in the training data, and the expansion of the number of data points for MHC‐II molecules already included in the old training data. Figure 4(c,d) quantifies that both of these factors indeed contribute to the performance gain. Figure 4(c) shows the performance gain as a function of the change in distance of the query molecule to the nearest neighbour of the training data. From this plot, we see that the gain in predictive performance is related to a decrease in the nearest neighbour distance, and hence directly related to the inclusion of binding data for additional MHC‐II molecules in the new data set. Figure 4(d) shows the performance gain as a function of the change in the number of data points between the two data sets used for training. We here only include molecules shared between the two data sets used for training NetMHCIIpan‐3.1 and NetMHCIIpan‐3.2, as we in the previous analysis demonstrated how the distance to the nearest neighbour influences the performance. Figure 4(d) shows that the gain in performance is correlated to change in the number of data points for the given MHC molecules. This indicates that the performance gain of the new NetMHCIIpan version is also driven by the increase in the number of data points for molecules already included in the 2013 data set. The one data point in Figure 4(c,d) with increased nearest neighbour distance and decreased number of data points corresponds to the HLA‐DPA10103‐DPB10201 molecule for which faulty data were removed in the 2016 data set.

Comparing NetMHCII and NetMHCIIpan

Using the data points shared by the old and updated data sets, we first compared the different versions of NetMHCII and NetMHCIIpan. We showed how the new versions of the methods outperformed the previous versions for both NetMHCII and NetMHCIIpan. We then evaluated the performance of the two versions of the methods using only ‘new’ peptides, for the MHC molecules covered both by the old and the updated data sets. The result of this analysis showed that both methods on this data set gained a significant improvement in the predictive performance, supporting the importance of expanding the size of the training data even for MHC molecules already characterized by binding data. When evaluating new peptides one has to keep in mind that MHC binding predictors are often used to select peptides for experimental validation and new data sets may be less diverse than historic data sets generated sampling the entire space of a given set of protein sequences.34

The main difference between NetMHCII and NetMHCIIpan is that NetMHCII is an allele‐specific method trained separately for each MHC molecule, whereas NetMHCIIpan is a pan‐specific method that contains a single ensemble of networks using information from all MHC molecules in the data set. We would therefore expect that the allele‐specific method outperforms the pan‐specific method for MHC molecules where sufficient data are available to accurately characterize the binding motif, and we would expect the pan‐specific method to outperform the allele‐specific method when data are scarcer. This is exactly what we observed when we compared the predictive performances of NetMHCII‐2.3 and NetMHCpan‐3.2. Earlier work has shown a similar result, namely that when allele‐specific neural network prediction algorithms rely on a sufficient number of peptide binders to achieve high predictive performances.33, 35 This illustrates how the allele‐specific method is preferable only if a large amount of data is available for the MHC molecule in question, but highlights the strength of the pan‐specific methods, which can benefit from the data of related MHC molecules to make reliable predictions for MHC molecules with limited data. Because of this difference between the allele‐specific and pan‐specific methods, we implemented a simple combination of two methods as this has been shown to improve the predictive performance for MHC‐I molecules.33 This analysis showed that NetMHCIIpan‐3.2 outperforms NetMHCII‐2.3 for MHC molecules, which has been trained with very few peptides, but that a combination of the predictions from the two MHC‐II methods still outperformed each individual method.

Leave‐one‐molecule‐out performance for NetMHCIIpan

One of the main powers of the NetMHCIIpan method is that it can predict binding affinities for uncharacterized MHC molecules. To assess the performance of the method in such a task, we constructed a LOMO experiment where we tested the performance of the NetMHCIIpan method for predicting binding affinity for MHC molecules not included in the training data for the method. From this analysis, we could show that the pan‐specific method is capable of prediction binding affinity for MHC molecules where no binding affinity data are available and further demonstrate that the predictive performance is dependent on the distance to the nearest neighbour. This last observation indicated that the predictive performance of the NetMHCIIpan method could be further improved by including more uncharacterized MHC molecules into the training data and it is therefore important to generate experimental peptide binding affinity data points in a targeted fashion for MHC molecules not yet characterized.

Distance tree for HLA class II molecules

To understand the different groups of HLA class II molecules, we generated a fictional distance tree using NetMHCIIpan‐3.2. The groups shown in this distance tree can be used to understand how peptides interact with different MHC molecules and can be used to discriminate between binders and non‐binders. The distance tree can also be used to identify T‐cell epitopes with similar properties important for the design of epitope‐based vaccines. Another aspect that can be observed for the tree is that most MHC molecules have strong anchor positions at P1, P4, P6 and P9, which have also been observed in previous studies.8

T‐cell epitope benchmark

Accurate predictions of peptide binding affinities to MHC molecules are important for understanding the cell‐mediated immune response and for generating better screening methods for cost‐effective identification of immunogenic peptides. We therefore wanted to test the predictive performance of the two versions of NetMHCIIpan on a T‐cell epitope data set, and doing this we demonstrated how the new version of NetMHCIIpan obtained a significantly improved predictive performance compared with the earlier version. Two main factors explain this performance gain: (i) including data for new MHC‐II molecules decreases the distance to the nearest neighbour, (ii) including an increased number of data points allows the method better characterizing the specificity of a given MHC‐II molecule.

In conclusion, we believe that NetMHCII and NetMHCIIpan can be used to improve MHC‐II binding predictions and reduce experimental costs for immunologists working within the field of epitope‐based vaccine design, and to improve our knowledge about the peptide–MHC interaction, a key event in the cellular immune response.