NetMHCIIphosPan: a machine learning tool for predicting HLA class II antigen presentation of phosphorylated peptides

Abstract

Phosphorylated peptides presented by human leukocyte antigen (HLA) class II molecules play pivotal roles in immune regulation, yet their characterization and prediction remain challenging due to data noise and limited HLA coverage. Here, we introduce NetMHCIIphosPan, a prediction method for HLA-II antigen presentation of phosphorylated peptides, developed using mass spectrometry (MS)-based immunopeptidomics datasets. Employing a refined peptide identification workflow, we reanalyzed earlier HLA-II phospholigand datasets and trained predictive models, achieving superior performance compared to models trained on original data. Binding motif analysis revealed that HLA-specific preferences for phospholigands closely aligned with those of unmodified ligands. Incorporating unmodified ligands into training further enhanced predictive accuracy, particularly for HLA-DP and HLA-DQ molecules. NetMHCIIphosPan outperformed existing tools, such as NetMHCIIpan-4.3 and MixMHC2pred-1.3, for prediction of HLA antigen presentation of phosphorylated peptides demonstrating robustness and utility. This work establishes NetMHCIIphosPan as a state-of-the-art tool for understanding the HLA-II phospholigandome, with potential applications in immunotherapy and vaccine design.

Keywords: HLA class II; antigen presentation; immunoinformatics; immunopeptidomics; machine learning; neural networks; phosphorylation; post-translational modification; predictive methods.

Fig. 1.. Overview of the phospholigand datasets identified by MaxQuant and PEAKS software.

(A) Peptide length distribution. In grey, the distribution of unmodified HLA-II ligands lengths from the training dataset of the NetMHCIIpan-4.3 method (“UM_NetMHCIIpan”). (B) Phosphorylation type prevalence, considering all phosphosites found in modified peptides. (C) Distribution of the number of phosphosites per modified peptide. (D) Number of peptides identified in each of the analyzed datasets. Note that the PvanBalen and Saghar datasets were not analyzed using MaxQuant, and MaxQuant values are hence not available for these datasets in D.

Fig. 2.. Five-fold cross-validation performance of the methods trained with phospholigands identified by MaxQuant and PEAKS.

M1 refers to the method trained and evaluated on the complete MaxQuant dataset, M2 refers to the method trained and evaluated on the MaxQuant dataset, pre-filtered according to MixMHC2pred-2.0 predictions (%rank≤10), and M3 refers to the method trained and evaluated on the complete PEAKS dataset. Performance metrics included are (A) ROC-AUC, (B) ROC-AUC 0.1, (C) PPV and (D) PR-AUC. ROC-AUC refers to the area under the Receiver Operating Characteristic curve; ROC-AUC 0.1 is the partial AUC integrated up to a FPR of 10%; PR-AUC is the area under the Precision-Recall curve; and PPV (Positive Predictive Value) is calculated as described in Materials and Methods. Performance is evaluated in a per-dataset manner from the concatenated phosphopeptide test sets. (E) Proportion of predicted non-binders (%rank>20) according to each of the methods. The center line inside the box indicates the median value of the plotted metric and the box covers the interquartile range. The whiskers extend to, at most, 1.5-fold of the interquartile range. The individual data points in the jitter plots correspond to different cell lines. Only cell lines with 10 or more positive instances for all trained methods were included in the analysis (N=31 cell lines). Y-axis ranges differ between panels A and B-D.

Fig. 3.. Concordance between unmodified and phospholigand binding motifs for HLA-II molecules with 30 or more peptides assigned during cross-validation.

(A) PSFMs (position-specific frequency matrices) were computed to represent HLA-II binding motifs for unmodified and phosphopeptides, based on cross-validation predictions of the different methods. The three methods included are M1, M2 and M3 described in Fig. 2. In the case of phosphopeptides, phosphosites were ignored and frequencies per position were renormalized to sum up to 1. For unmodified peptides, predictions from NetMHCIIpan-4.3 were employed. To compare binding motifs for a given allele, PCCs (Pearson Correlation Coefficient) were calculated between flattened unmodified and phosphopeptide PSFMs. The violin plots show the distribution of PCCs for each of the trained methods (each data point corresponds to a different HLA-II molecule, the diamond indicates the median PCC, and, above each plot, the analyzed number of HLA-II molecules is noted). (B) Sequence logos for a subset of four selected DRB1 alleles, generated from cross-validation predictions as in (A). The first row corresponds to logos illustrating phospholigand (“PH”, outlined in violet) binding motifs for 4 HLA-II molecules according to the M2 method, while the second row shows the same for the M3 method. The third row shows logos representing unmodified (“UM”, outlined in grey) HLA-II ligand binding motifs according to NetMHCIIpan-4.3. Each logo displays in parenthesis the number of peptides employed for its generation. In the case of the M3 and NetMHCIIpan-4.3 methods, logos were constructed by randomly subsampling 200 peptides per allele from the cross-validated test set predictions. Amino acid color coding: basic in blue (RHK), acidic in red (DE), polar in pink (STNQYC), hydrophobic in dark grey (GPLIVAMFW), and phosphosites in violet (sty).

Fig. 4.. Improvement of the baseline method trained with phospholigands by adding unmodified ligands to the training dataset and allowing for peptide inversion.

Five-fold cross-validation performance of the M4 method, trained only with phospholigands (“PH”), the M5 method, trained with both phospholigands and unmodified ligands (“PH+UM”), and the same method, M6, trained with peptide inversion (“PH+UM (w_inv)”; N=96 cell lines). Performance metrics included are (A) ROC-AUC, (B) ROC-AUC0.1, (C) PPV and (D) PR-AUC. ROC-AUC refers to the area under the Receiver Operating Characteristic curve; ROC-AUC 0.1 is the partial AUC integrated up to a FPR of 10%; PR-AUC is the area under the Precision-Recall curve; and PPV (Positive Predictive Value) is calculated as described in Materials and Methods. Performance evaluation, computed over the concatenated phosphopeptide test sets, and plots were made as in Fig. 2. One-tailed binomial tests were performed to compare the performance of the methods, with arrowheads indicating the direction of the tests (* P<0.05, ** P<0.01, and **** P<0.0001; ns, not significant). Y-axis ranges differ between panels A and B-D. (E) Concordance between unmodified and phospholigand binding motifs, based on model cross-validation predictions, for HLA-II molecules with 30 or more peptides assigned discriminated by loci. (plots were made as in Fig. 3A; above each plot, the number of analyzed HLA-II molecules is indicated).

Fig. 5.. Five-fold cross-validation performance of NetMHCIIphosPan and NetMHCIIpan-4.3 on training dataset phospholigands.

Performance metrics included are (A) ROC-AUC, (B) ROC-AUC0.1, (C) PPV and (D) PR-AUC. ROC-AUC refers to the area under the Receiver Operating Characteristic curve; ROC-AUC 0.1 is the partial AUC integrated up to a FPR of 10%; PR-AUC is the area under the Precision-Recall curve; and PPV (Positive Predictive Value) is calculated as described in Materials and Methods. Performance is evaluated in a per-dataset manner from the concatenated test sets (N=96 cell lines, for details see Materials and Methods). The center line inside the box indicates the median value of the plotted metric and the box covers the interquartile range. The whiskers extend to, at most, 1.5-fold of the interquartile range. The individual data points in the jitter plots correspond to different cell lines. One-tailed binomial tests were performed to compare the performance of the NetMHCIIphosPan method against the NetMHCIIpan-4.3 method (different “phospho-equivalent” datasets) (**** P<0.0001). Y-axis ranges differ between panels A and B-D.

Fig. 6.. Sequence logo representations of HLA-II binding preferences for unmodified and phosphorylated peptides.

Ten alleles were selected from among the 25 with the highest number of assigned phospholigands during NetMHCIIphosPan training. To construct the motifs two sets of 200,000 phospho- and unmodified peptides each were randomly sampled and predicted with NetMHCIIphosPan and NetMHCIIpan-4.3 respectively. The top 0.1% scoring peptides were used to construct the logos (all peptides had a predicted %rank≤5, and phosphopeptides without phosphosites in the binding core were excluded). Amino acid color coding: basic in blue (RHK), acidic in red (DE), polar in pink (STNQYC), hydrophobic in dark grey (GPLIVAMFW), and phosphosites in violet (sty).

Fig. 7.. Phosphosite enrichment/depletion within the NetMHCIIphosPan predicted 9-mer binding core, discriminated by allele (columns) and phosphorylation type (rows).

HLA-II binding predictions were computed for 200,000 random natural phosphopeptides for alleles with 100 or more phospholigands assigned during model cross-validation (in total 50). For each of the selected alleles, the top 2.5% scoring predictions were taken as binders, and an equal number of randomly sampled phosphopeptides with a %rank>20 were considered non-binders. The y-axis shows the log2 ratio between the relative frequency of binders vs. non-binders for each of the 9 binding core positions. Each data point in the plot corresponds to a different allele and, data points colored in blue indicate an enrichment of binders relative to non-binders (the blue dashed line corresponds to a log2 ratio of 1) whereas data points in red show the opposite, a depletion of binders relative to non-binders (the red dashed line corresponds to a log2 ratio of −1). Supplemental Table S5 includes the computed log2(Ratio) for each binding core position, HLA-II molecule, and phosphorylation type.

Fig. 8.. Performance of NetMHCIIphosPan and MixMHC2pred-1.3 methods on two phospholigand external evaluation datasets covering HLA-DR alleles.

Performance metrics included are (A) ROC-AUC, (B) ROC-AUC0.1, (C) PPV and (D) PR-AUC. ROC-AUC refers to the area under the Receiver Operating Characteristic curve; ROC-AUC 0.1 is the partial AUC integrated up to a FPR of 10%; PR-AUC is the area under the Precision-Recall curve; and PPV (Positive Predictive Value) is calculated as described in Materials and Methods. Performance is evaluated in a per-cell line manner from the concatenated test sets. Metrics were computed for cell lines with at least 10 phospholigands (N=5 cell lines for “Spike_DR” and N=11 cell lines for “Esophageal_Tumor_DR”). The center line inside the box indicates the median value of the plotted metric and the box covers the interquartile range. The whiskers extend to, at most, 1.5-fold of the interquartile range. The individual data points in the jitter plots correspond to different cell lines. A one-tailed binomial test was employed for model performance comparison, and p-values were estimated using a bootstrap analysis, for more details refer to Materials and Methods. (**** P<0.0001). Y-axis ranges differ between panels A and B-D.