NetMHCIIphosPan: A Machine Learning Tool for Predicting HLA Class II Antigen Presentation of Phosphorylated Peptides

HLA-II MS Data Sets

In this study, we employed HLA class II (HLA-II) phospholigands derived from various MS-based immunopeptidomics experiments. The PEAKS DeepNovo Peptidome, version 11.5 (Bioinformatics Solutions, Waterloo, Canada), hereafter referred to as PEAKS, was used to search the raw MS data. In the final list of identified peptides, only peptides found by database search were included in the study. That is, peptides identified by homology search (peptides with mutations) or de novo sequencing were excluded.

Raw MS data from 20 single-allelic (SA) cell lines from Abelin et al.²⁶ and 23 multiallelic (MA) cell lines from Racle et al.²⁵ were used in this search. Sequential HLA purification experiments were performed on 12 of the 23 multiallelic cell lines, first with an HLA-DR monoclonal antibody (termed “Racle_DR” cell lines) and then with a pan-HLA-II antibody (termed “Racle” cell lines), yielding different subsets of ligands. Of the remaining, 10 multiallelic cell lines were only processed with a pan-HLA-II antibody, and one was only processed with an HLA-DR monoclonal antibody. This resulted in a total of 35 subsets of ligands derived from multiallelic cell lines.

Additionally, we included MS data sets from previous studies by Nilsson et al.,^5,6 Kaabinejadian et al.,⁴ and van Balen and colleagues.²⁷ These data sets that were derived from either genetically engineered single-allelic (SA) or homozygous cell lines included 16 HLA-DQ, 14 HLA-DP, and 11 HLA-DR data sets from which HLA class II was purified using locus-specific monoclonal antibodies, SPVL3, B7/21, and L243, respectively. This allowed us to expand our search to cover the HLA-DP and HLA-DQ as well.

The identified sequences were searched against the Swiss-Prot Homo sapiens (Taxon identifier 9606) database, including 20,395 entries, accessed on 12 May 2021 (https://www.uniprot.org). Precursor mass error tolerance of 10 ppm and fragment mass error tolerance of 0.02 Da were used. The enzyme was set to none and HCD (High-energy collisional dissociation) was selected as the fragmentation method. Variable post-translational modifications (PTM) including acetylation (N-term), deamidation, pyroglutamate formation, oxidation, sodium adducts, phosphorylation (STY), cysteinylation and citrullination were included in the database search. The maximum number of PTMs per peptide was set at 3 and the peptide length was restricted to 7 to 30 amino acids. Identified peptides were further filtered at a false discovery rate (FDR) of 5% using the PEAKS decoy-fusion algorithm, in line with earlier work,²⁴ to allow for an increased peptide yield. Default settings were used for the rest of the parameters including confident amino acid threshold and DeepNovo Score. For data acquired using the Triple-TOF Sciex 5600 instrument,^4–6 CID (Collision-induced dissociation) was selected as the fragmentation method, and precursor mass error tolerance was set at 30 ppm. Other parameters remained the same as above, and only peptides found by database search (DB Search) were included in the study.

In addition to the PEAKS analysis, the raw files from the Racle data sets were reanalyzed using FragPipe v23.1, incorporating MSFragger v4.3 and Philosopher v5.1.2.^18,28,29 MSFragger searches were performed using the same Swiss-Prot Homo sapiens database (accessed 12 May 2021) to ensure comparability with PEAKS output. Precursor mass tolerance was set to ±10 ppm and fragment mass tolerance to 0.02 Da. The enzyme specificity was set to “nonspecific”. Variable modifications were oxidation (M), phosphorylation (STY), cysteinylation (C), and protein N-terminal acetylation, with a maximum of 2 variable modifications per peptide. Peptide lengths were restricted to 7 to 30 amino acids. Identified peptides were filtered with sequential 1% FDR at PSM- and protein-levels using Philosopher’s implementation of PeptideProphet. MS1-level peptide quantification was performed using IonQuant,³⁰ with default label-free quantification (LFQ) settings, including MaxLFQ, match-between-runs (MBR), and intensity normalization across runs. To ensure methodological consistency and enable a fair comparison with FragPipe, we conducted a second PEAKS search using the same search parameters.

All identified phosphopeptide sequences of length 12–21 amino acids are listed in Supplementary Table S1.

Training Data Sets

For all trained methods, negative data were generated as random natural phosphopeptides sampled for each phospholigand data set, covering the same peptide length range as the latter, but with equal numbers of random peptides for each peptide length. More precisely, this number of random peptides was set to five times the number of ligands for the most prevalent peptide length. In this manner, the algorithm can learn the difference in peptide length preferences for ligands, as described before.^31,32 The random phosphopeptides were sampled from the human phosphoproteome downloaded from PhosphoSitePlus,³³ a comprehensive MS database of post-translational modifications (PTMs).

Several methods were trained on different data sets. Models M1, M2 (and its variants M2′ and M2″) were trained on the original data from Solleder et al.²⁴ Model M3 was trained on the same original MS data reanalyzed using the workflow defined here. Models M4-M6 were trained using an extended data set, expanding the data set from M3 with phospholigands obtained from earlier MS studies (Saghar_DR/DP/DQ^4–6 and PvanBalen_DP²⁷) reanalyzed using the workflow defined here. M4 was trained using only the extended phospholigand data set. For M5 and M6, the training data was extended by unmodified MS eluted ligands derived from human cell lines or alleles, also from the NetMHCIIpan-4.3 training data set. These data sets consist of single-allelic (SA) cell lines (101,449 peptides from 70 human cell lines covering 52 HLA-II alleles) and multi-allelic (MA) cell lines (44,956 peptides from 37 human cell lines covering 52 HLA-II alleles). For all models, a set of unmodified peptide binding affinity (BA) data extracted from the training data set of the NetMHCIIpan-4.3 method (125,985 peptides covering 72 HLA alleles) was included in the training data. For the models M5 and M6, phospholigands and random natural phosphopeptides were upscaled (i.e., duplicated) 10 times when incorporating them into the training data set.

Data Partitioning

As described below, models were trained and evaluated using 5-fold cross-validation. Before generating the five training data partitions, phosphosites represented as “sty” (where “s” stands for pS, “t” for pT, and “y” for pY) were translated into their corresponding uppercase “STY” amino acids. Next, a common motif algorithm, inspired by the original Hobohm 1 algorithm, was used to cluster peptides sharing any identical 9-mer subsequences.³⁴ This ensures that nested peptides (ignoring phosphosites) are placed in the same partition. After partitioning, the phosphosites were back-translated to “sty”. The final training data set is available at https://services.healthtech.dtu.dk/services/NetMHCIIphosPan-1.0/, “Supplementary material” tab.

Model Training

Models were trained using a 5-fold cross-validation scheme, employing the NNAlign_MA algorithm.³⁵ The model architecture was set to 40 and 60 hidden neurons, and each architecture was trained with 10 random weight initializations per cross-validation fold, resulting in a final ensemble consisting of 100 (5 × 2 × 10) network models. All models were trained using backpropagation with stochastic gradient descent, for 400 epochs, without early stopping, and a fixed learning rate of 0.05. Peptides were represented using sparse sequence encoding with an extended amino acid alphabet: ARNDCQEGHILKMFPSTWYVsty, where “sty” represents the phosphorylated amino acids pS, pT, and pY as described above. The motif length, or binding core, was set to 9 amino acids. Only SA data was included in the training for a burn-in period of 30 epochs, followed by training cycles including both SA and MA data. During the first 10 training iterations, only ILVMFYWRKsty amino acids were allowed at position 1 of the binding core. Models were trained including peptide flanking region (PFR) amino acid composition encoding, PFR length encoding, and peptide length encoding as described earlier.³² If required, models were trained allowing for peptide inversion as described in ref 6.

Independent Benchmark Data Sets

External phospholigand data sets were gathered to benchmark the NetMHCIIphosPan method against other state-of-the-art tools, such as MixMHC2pred-1.3.²⁴ Two distinct class II immunopeptidomics data sets were selected, both acquired in the laboratory of Nicola Ternette. The first data set consisted of unpublished data generated on Orbitrap Fusion Lumos Tribrid Mass Spectrometer from fresh frozen tissue surgically resected from 14 subjects diagnosed with esophageal adenocarcinoma (EAC), tissue consisting of both tumor and normal adjacent esophagus and stomach. The second data set consisted of previously published class II immunopeptidomes data from monocyte-derived dendritic cells (MDDCs) generated on a Q-Exactive HF-X Hybrid Mass Spectrometer from 5 HLA-DRB1-heterozygous healthy blood donors.³⁶ Here, each donor MDDC was pulsed with a recombinant SARS-CoV-2 S protein vaccine immunogen candidate. Both tissue and MDDC class II immunopeptidomes were generated using the anti-HLA-DR monoclonal antibody (L243). In the case of the “Esophageal tumor” data set, the HLA typing only covers DRB1 alleles (and HLA-DQ alleles, though these were not employed for making predictions). For these data sets, DRB3/4/5 alleles were inferred with the HLAAssoc-1.0 tool.³⁷ Regarding the “Spike” data set, full HLA typing is available for these samples.

Preparation of Peptides from Tissues and Cells

Tissues were placed in a lysis buffer (0.5% (v/v) IGEPAL 630, 50 mM Tris pH 8.0, 150 mM NaCl, and 1 tablet complete Protease Inhibitor Cocktail EDTA-free (Roche) per 10 mL buffer) in bead beater tubes for soft tissues (Precellys Evolution, Bertin Technologies, France), bead beating was performed with 5 cycles at 7200 rpm for 20 s, separated by 20 s pauses followed by mixing for 30 min at 4 °C to obtain complete lysis. Lysates were cleared by centrifugation initially at 500 × g for 10 min and supernatants were removed and centrifuged at 21,000 × g for 45 min at 4 °C to complete clarification. Peptide-HLA complexes were immunoprecipitated with L243 mAb cross-linked protein G sepharose beads at 4 °C overnight. The resin was washed consecutively in Tris buffer 50 mM, pH 8.0 containing 150 mM NaCl, 5 mM EDTA; 150 mM NaCl; 450 mM NaCl; no addition. Peptides were eluted twice with 500 μL 10% acetic acid and dried in a speed vac before cleanup by HPLC as described in ref 36.

LC-MS Data Acquisition

The SARS-CoV-2 data set was generated as described,³⁶ using an Ultimate 3000 RSLC nano system (Thermo Scientific) coupled to a Q-Exactive HF-X Hybrid Mass Spectrometer (Thermo Fisher Scientific). For the EAC tissue, purified peptides were initially trapped and separated; spectra were acquired using an Ultimate 3000 RSLC nano system (Thermo Scientific) coupled to an Orbitrap Fusion Tribrid Mass Spectrometer (Thermo Fisher Scientific) with a nanoelectrospray ion source. Separation was achieved on a PepMap C18 EASY-Spray column (3 μm particle size, 75 μm × 50 cm) (Thermo Scientific) in buffer A (0.1% formic acid) using a linear gradient of 2%–35% buffer B (acetonitrile and 0.1% formic acid) at a flow rate of 300 nL/min over 60 min. The ion transfer tube temperature was 305 °C and full MS spectra were acquired (300 to 1500 m/z) in the Orbitrap at 120,000 resolution with an automatic gain control (AGC) target of 40,000. Precursors were isolated using an isolation width of 1.2 amu and MS2 scans were acquired at maximal speed (120 ms with an AGC target of 300,000). Higher-energy collisional dissociation (HCD) with NCE of 28 was used to fragment peptides with a charge state of 2–4, for singly charged peptides NCE of 32 was applied with a lower priority in the decision tree. Fragments were analyzed in the Orbitrap at 30,000 resolution.

Phosphopeptide Data Analysis and Benchmark Data Set Generation

Both data sets were analyzed with the PEAKS DeepNovo Peptidome 11.5 platform, applying the same search parameters and settings mentioned before for the training data sets (refer to the section entitled “HLA-II MS Data Sets”). The list of identified peptide sequences can be found in Supplementary Table S2.

As all external evaluation data originated from MS-based experiments, each data set was enriched with random natural phosphopeptides in the same manner as for the training data. For the independent “Esophageal tumor” data set, normal and tumor samples from each subject were combined into a single data set. Only data sets with at least 10 positives were retained, resulting in a benchmark comprising 2058 phospholigands (available at https://services.healthtech.dtu.dk/services/NetMHCIIphosPan-1.0/, “Supplementary material” tab).

Performance Evaluation

When assessing the performance of NetMHCIIpan-4.3 on our method’s training data set and when evaluating the performance of NetMHCIIpan-4.3 and the method developed here on the external data sets, the following strategy was employed.

For NetMHCIIpan-4.3, which was trained exclusively on unmodified ligands, three “phospho-equivalent” data sets were generated: (1) phosphosites were replaced with glutamic acid (“E”) due to its similar charge, (2) phosphosites were substituted with amino acids of similar charge (pS → D, pT → E, pY → Y) (“DEY”), or (3) phosphosites were reverted to their unmodified forms (“STY”). Next, the overlap between the corresponding “phospho-equivalent” data set and the NetMHCIIpan-4.3 training data set was computed using the common motif algorithm.³⁴ For peptides with a 9-mer overlap with peptides placed in a given NetMHCIIpan-4.3 training partition, predictions were made using the neural network ensemble that was not trained on that specific partition (NetMHCIIpan-4.3 was trained using a 5-fold cross-validation scheme, which implies that 20 out of the total 100 networks were used to make predictions). If no overlap was found, a random neural network ensemble (also consisting of 20 networks) was selected for predictions.

To evaluate our method, which was trained on both unmodified and phosphopeptides, phosphopeptides in the external data sets were converted from “sty” to “STY”. Then, the same evaluation procedure described before was followed: predictions were generated using the neural network ensemble trained on orthogonal data (as defined in the “converted” peptide space) to the evaluation data to prevent artificially inflated performance values.

Performance Metrics

The model’s ability to distinguish ligands from random peptides was assessed using four performance metrics: the ROC-AUC (Area Under the Receiver Operating Characteristic Curve), the ROC-AUC0.1 (ROC-AUC integrated up to a False Positive Rate of 10%), the PR-AUC (Area Under the Precision-Recall Curve), and PPV (Positive Predictive Value). Note that ROC-AUC0.1 values were not normalized using the McClish standardization; therefore, a random performance for this metric does not correspond to 0.5. The PPV was computed as the proportion of true positives among the top N ranked peptides, where N corresponds to the number of ligands (positive instances) for a given HLA or cell line. In 5-fold cross-validation, the data was partitioned into five subsets. For each fold, models were trained on four-fifths of the data and evaluated on the remaining one-fifth, which served as the test set. This procedure was repeated five times, with each subset used once as the test fold. The resulting five subsets of test predictions, which together covered the entire data set were concatenated, and all performance metrics were computed on the combined set of “raw” prediction scores. For the independent benchmark, a similar approach was applied. For each peptide, models from a single training fold identified from a 9-mer overlap (with phosphosites being translated into their corresponding uppercase “STY” amino acids) to the training partitions were used for the predictions. In case no overlap was found, models from a random training fold were used. Next, metrics were similarly calculated using “raw” prediction scores generated for each HLA allele or cell line. All metrics were reported only for data sets or cell lines containing at least 10 positive instances.

Statistical Rationale

The corresponding statistical tests employed to compare model performances or other results are described in detail in the legend of each of the figures. In all cases, p-values less than 0.05 are considered to be statistically significant. In Figure 8, pairwise comparisons of model performances were made using one-tailed binomial tests. For each comparison, 10,000 resamples of the entire predicted benchmark data set were generated, allowing for replacement (each subsample should contain at least 1% of the total original positive instances). Next, the p-value for one method outperforming the other was estimated as the fraction of resampled data sets where the other method had superior performance.

Data Visualization

Sequence logos were made using the R “ggseqlog” package developed by Wagih et al.³⁸ All the remaining plots in this study were made in R using the “ggplot2” package.

Data Overview: Improved Phosphopeptide Identification

As a starting point, we extracted the original HLA-II phosphorylated ligand data set employed in Solleder et al.²⁴ In this work, the authors reprocessed high-quality MS samples from Abelin et al.²⁶ and Racle et al.²⁵ using oxidation (M), N-terminal acetylation and phosphorylation (STY) as PTMs in their spectral searches with the MaxQuant software (version 1.5.5.1), and from these HLA-II phospholigandome data characterized the binding motifs for more than 30 HLA-II molecules. A main challenge, described in the work by Solleder et al., concerning this HLA-II phospholigandome data set, is the presence of co-eluted contaminants or incorrectly identified peptides resulting in approximately half of the data being discarded by the authors.

We speculated that a substantial source of the seeming low quality of data generated in the study by Solleder et al. stems from a rather old version of the MaxQuant software used in the study and hypothesized that the issue could at least in part be corrected by using a more recent method for MS spectral annotation. Many such novel tools exist, and we do not intend to suggest an optimal choice. Rather, we simply picked one such tool, PEAKS, for the task given our in-house expertise and experience with its usage. Based on a series of pilot studies, we define a workflow tailored to discover PTM peptides with properties aligning with MHC class II antigen presentation (peptide length distribution and MHC binding potential as defined from models trained on unmodified peptides). Using this workflow, we next reprocessed the raw MS samples from the Solleder et al. study. Further, we included other MS samples in our searches (Saghar_DR/DP/DQ^4–6 and PvanBalen_DP²⁷) to expand the allele coverage to better characterize HLA-DP and DQ binding motifs (refer to Materials and Methods for more details on data sets). Using this approach, we were able to increase peptide numbers and HLA coverage of the data in comparison to the original publication. Figure 1 provides an overview of the data sets reprocessed with MaxQuant by Solleder et al. and the same data sets reanalyzed with PEAKS (refer to Supplementary Figure S1A–C for the statistics over the full PEAKS data sets). Overall, the MaxQuant and PEAKS data sets display similar peptide length distributions preserving the peak at 15 amino acids long (Figure 1A), although the former is shifted toward shorter peptide lengths. The peptide length distribution of unmodified HLA-II ligands falls between the distributions of these two data sets. The distribution of phosphosites in the two data sets is pS (55%/50%), pT (34%/36%), and pY (11%/14%) (for MaxQuant/PEAKS data sets, Figure 1B), in line with the proportions of observed phosphosites in human proteins reported in PhosphoSitePlus³⁹ (pS: 58%; pT: 25%; pY: 17%). Most of the phosphorylated ligands (85%/88%) contain only one phosphosite, while a small fraction contain two (13%/11%), and almost none contain three (2%/1%) phosphosites (for MaxQuant/PEAKS data sets, Figure 1C). PEAKS recovers 33%, 246%, and 118% more phosphopeptides for the Abelin, Racle (data generated with pan class II antibody in the immunoprecipitation), and Racle_DR (data generated with a DR-specific class II antibody in the immunoprecipitation) data sets, respectively, in comparison to the original MaxQuant analysis (Figure 1D). Considering these three data sets reanalyzed in common, Abelin, Racle, and Racle_DR, the overlap between phosphopeptides identified by PEAKS and MaxQuant was limited. In the Abelin data set, only 15 phosphopeptides were identified by both tools, corresponding to 2.0% of MaxQuant and 1.0% of PEAKS identifications. For the combined Racle and Racle_DR data sets, 499 overlapping peptides were observed, representing 15.1% of MaxQuant and 4.8% of PEAKS identifications. Across these three data sets, only 3.8% of the total unique identified phosphopeptides overlap between the two approaches. Taking all data sets into account, 19,019 phosphopeptides were identified with the PEAKS software, while 4,632 phosphopeptides were identified with MaxQuant (including only 12 to 21-mers in both cases).

We acknowledge that the comparison is limited by the outdated version of MaxQuant used in the earlier study. We will return to this later in the manuscript, where we include a study using a third state-of-the-art algorithm for immunopeptidome identifications, and here demonstrate a high consistency to the PEAKS results. In summary, these results suggest that using the workflow defined here results in a major increase in peptide yields of identified phospholigands compared to that of earlier workflows.

Enhanced Predictive Models for Phospholigand Binding

Next, to investigate the quality of the data sets and assess to what degree each could lead to an understanding of the HLA-II phospholigandome, models were trained either with the original data from the Solleder et al. study (Abelin, Racle, and Racle_DR data set reanalyzed using MaxQuant) or with data from the same EL data sets reanalyzed using the workflow defined here.

All developed models were trained using the NNAlign_MA machine-learning framework.³⁵ In short, the NNAlign_MA algorithm enables simultaneous training and motif deconvolution of multi-allelic (MA) data, i.e., MS eluted ligands from cell lines expressing multiple HLA-II alleles. This aspect is crucial for this work since a substantial amount of the phosphorylated HLA-II ligands employed in the training of these methods originate from multi-allelic data, which implies they have no predetermined HLA-II allelic restriction. Indeed, this becomes even more relevant for this study as the HLA-II binding preferences of phosphorylated peptides have been very scarcely characterized and may differ from the HLA-II binding motifs described for their unmodified counterparts.

Models were trained using a 5-fold cross-validation scheme. To prepare the training data sets, random negatives were sampled from the human phosphoproteome deposited in PhosphoSitePlus,³³ and peptides were grouped using a common motif algorithm, ensuring no 9-mer overlap between partitions³⁴ (for details on Training Data Sets, Data Partitioning or Model Training refer to Materials and Methods). Details related to the number of phospholigands (positive instances), random phosphopeptides (negatives) and the number of cell lines included in the training data set of these models, are given in Supplementary Table S3. For each fold, the model was trained on four-fifths of the data and evaluated on the remaining one-fifth. This process was repeated five times, resulting in five sets of test set predictions that, together, covered the entire data set. These predictions were concatenated, and metrics were computed on the combined set of prediction scores. We report four complementary performance metrics: ROC-AUC, ROC-AUC0.1 which corresponds to ROC-AUC integrated up to a 10% FPR, Precision-Recall AUC (PR-AUC), and Positive Predictive Value (PPV) computed over the top N predictions, where N corresponds to the number of HLA-II ligands for each data set. ROC-AUC0.1 focuses on the high-specificity region of the ROC curve, which is critical when minimizing false positives, while PPV directly quantifies the model’s ability to prioritize true HLA-II ligands among top-ranked predictions, where experimental validation efforts are concentrated. Note that the ROC-AUC0.1 values reported are not normalized and that a random performance for this metric here corresponds to 0.05. Finally, it is relevant to mention that the NNAlign_MA framework requires both eluted ligand (EL) and binding affinity (BA) data to perform training. Therefore, in this case, the phospholigand training data sets were combined with more than 125,000 unmodified peptides measured in BA assays from the training data set of the NetMHCIIpan-4.3 method.

Figure 2 illustrates the 5-fold cross-validation performances of three methods, two of them, M1 and M2, trained with the MaxQuant data sets, and one trained with the PEAKS data sets (M3). For “MaxQuant_w_filter” M2 method, we applied the same filtering strategy as Solleder et al. using MixMHC2pred-2.0⁴⁰ predictions replacing phosphosites by glutamic acid to filter the data and discard potential contaminants or misidentified peptides (%rank >10). Supplementary Figure S2 expands this analysis including the performances of two other methods, trained on the MaxQuant data sets and filtered according to NetMHCIIpan-4.3 predictions.

We notice that the different methods are trained and evaluated on data sets of varying sizes, which complicates direct comparison. Nonetheless, the results presented in Figure 2 and Supplementary Figure S2 show that, although filtering of potential noise improved the performance of the models trained and evaluated on the MaxQuant data sets, the model trained on the full, unfiltered PEAKS data sets consistently exhibit a superior performance. This enhanced predictive power is also reflected in the overall lower proportion of predicted non-binders (phosphopeptides with %rank >20) within the PEAKS data sets, in relation to the MaxQuant and MaxQuant_w_filter data sets (refer to Figure 2E). This robustness further contributes to the superior performance even in cell lines with very few data points. For instance, the M3 model achieves an average ROC-AUC0.1 of 0.56 on the 10 data sets with the fewest positive examples (average of 130 positive data points per data set), while the M1 model achieves an average ROC-AUC0.1 of only 0.25 on the corresponding 10 data sets with the most positive examples (average of 183 positive data points per data set) (see Supplementary Table S4). This result suggests that the superior performance of the PEAKS-trained model is attributed not only to the increased data quantity but also a higher data quality.

Integrating Unmodified and Phospholigand Data

In light of these results, and in line with conclusions from earlier work on HLA class I,⁴¹ we next investigated whether the predictive performance of the M3 model could be improved by expanding the training set to include unmodified ligands. To achieve this aim, we constructed a new training data set by combining phospholigands with over 140,000 unmodified HLA-II ligands, extracted as a subset of the NetMHCIIpan-4.3 training data set. This combined set also included phospholigands derived from the additional van Balen and Saghar data sets described in Figure 1. To balance the number of unmodified ligands and phospholigands included in the training data set, phospholigands were upscaled by a factor of 10 (refer to Methods for more details). Furthermore, we also allowed for the inverted peptide binding described earlier for unmodified peptides for a subset of HLA-DP molecules.^6,40,42 For these molecules, the ligand repertoire has been shown to cluster into two mirror-symmetric motifs, consistent with peptide binding occurring either in the canonical N → C orientation or in a reversed C → N orientation. Supplementary Table S3 details the number of phospholigands, negative phosphopeptides, and cell lines used for model training.

Figure 4A–D shows the 5-fold cross-validation performance of the model M4, trained on phospholigands only, the model M5, trained on phospholigands and unmodified ligands, and the model M6, trained on the same data but including the option for the inverted binding motif. In terms of the ROC-AUC, ROC-AUC0.1, PPV and PR-AUC metrics, the M5 and M6 methods demonstrated close to equal performance and significantly outperformed the M4 method (one-tailed binomial tests, for the ROC-AUC all P < 0.05, and ROC-AUC0.1, PPV and PR-AUC all P < 0.0001). Next, we repeated the analysis to evaluate the similarity between unmodified ligand and phospholigand binding motifs for the studied methods, also based on model cross-validation predictions. As before, only HLA-II molecules with more than 30 assigned ligands were included. The results of this analysis are presented in Figure 4E, comparing the three methods across different loci. For HLA-DR, the distribution of PCC values is highly similar for all three methods. For the HLA-DQ alleles, there is a marked difference between the PCC distributions of the two models trained including unmodified data to the model trained without. Here, the M5 and M6 (including inversions) methods both display a median PCC > 0.77, while for the M4 method this value is 0.60. Finally, for HLA-DP, the M6 model shows higher PCC values compared to the two other methods, with a median PCC of 0.91. Inspecting sequence logos representing phospholigand binding motifs for selected HLA-DP and HLA-DQ molecules (refer to Supplementary Figure S5), it becomes clear that the motifs of the M5 and M6 methods mimic more closely the binding preferences found in the unmodified ligands (in this case represented by the NetMHCIIpan-4.3 logos). In particular, the sequence logos for the two exemplary HLA-DP molecules illustrate that the concordance between unmodified and phosphopeptide binding motifs is maximized for the M6 model, which allows for the inverted peptide binding mode, resulting in the consistent placement of “K” at P1 instead of P9. Supplementary Figure S6A shows that the presence of the inverted binding mode for phospholigands follows the same rules as earlier described for unmodified ligands. That is, this binding mode is almost solely present for HLA-DP molecules, and mainly for molecules sharing the HLA-DPA1*02:02/01 alleles (see Supplementary Figure S6B).

In summary, these results demonstrate that the inclusion of the unmodified ligands and the option for the inverted binding mode have aided the deconvolution of complex multiallelic phospholigand data sets, especially enhancing the characterization of HLA-DP and HLA-DQ binding motifs, which contributes to the overall moderate improvement in model performance. Moving forward, the M6 method will be referred to as NetMHCIIphosPan, our method of choice.

Benchmarking against the State-Of-The-Art Method: NetMHCIIpan-4.3

Next, we benchmarked the NetMHCIIphosPan method against its closest predecessor, NetMHCIIpan-4.3. To compare with the earlier developed method, solely trained on unmodified ligands, three “phospho-equivalent” training data sets were generated: (1) phosphosites were replaced with glutamic acid (“E”) due to its similar charge, (2) phosphosites were substituted with amino acids of similar charge (pS → D, pT → E, pY → Y) (“DEY”), or (3) phosphosites were reverted to their unmodified forms (“STY”).

Figure 5 shows that NetMHCIIphosPan significantly (all P < 0.0001) surpassed NetMHCIIpan-4.3, evaluated over the different “phospho-equivalent” data sets across all metrics (for performance evaluation details refer to Methods section).

Concordance between Unmodified and Phosphorylated Peptide Binding Motifs

Following, we assessed if the allele-specific binding motif preferences captured by NetMHCIIphosPan aligned with the motifs present in the training data. This analysis is essential to assess to what degree the model has learned the underlying motifs. Here, we compared logos for the 25 HLA-II molecules with the highest number of assigned phospholigands during training constructed either from the 5-fold cross-validation NetMHCIIphosPan predictions (excluding phospholigands with predicted %rank >5 and without phosphosites within the binding core in order to focus the analysis on peptides with very high antigen presentation likelihood only) or from the top 200 (0.1%) NetMHCIIphosPan predictions from a set of 200,000 random phosphopeptides (refer to NetMHCIIphosPan logos in Supplementary Figure S7). As before, we evaluated the pairwise similarity between unmodified and phospholigand logos for each individual HLA-II molecule by computing the Pearson correlation coefficient (PCC) between their corresponding flattened PSFMs. Supplementary Figure S8 shows the distribution of PCCs, with a median value of 0.90, confirming that the motifs captured by the model to a high degree align with the one found in the phospholigands training data sets.

Subsequently, we created sequence logos and visually compared the phospholigand HLA-II binding motifs learned by our method with the corresponding NetMHCIIpan-4.3 motifs for unmodified ligands. To achieve this aim, we, like above, predicted two sets of 200,000 random phospho- and unmodified peptides each and selected the top 200 (0.1%) binders. Figure 6 presents sequence logos for ten alleles selected from the top 25 with the highest number of assigned phospholigands during NetMHCIIphosPan training. Supplementary Figure S7 extends this analysis to include all top 25 HLA-II molecules, as previously described. This analysis again confirms the high concordance between the binding motifs of phosphorylated and unmodified ligands. Focusing on the phospholigand motifs, the logos reveal that phosphosites (shown in violet) are distributed throughout the binding motifs, though they predominantly localize to canonical HLA-II anchor positions.

Insights into Phosphosite Localization Preferences

Given this last observation, we next investigated in more detail if phosphosites shared a specific localization preference within the HLA-II ligands. This was evaluated first within the phospholigand with respect to the NetMHCIIphosPan predicted 9-mer binding core, and second within the predicted binding core with respect to the 9 possible binding positions. For the first analysis, we calculated the relative frequency of appearance of the phosphosite within the binding core or in its corresponding left or right peptide flanking regions (PFRs), averaged across peptides 13 to 17 amino acids long. This peptide length selection covers 69% of total phospholigands (refer to Figure 1A) and was intentionally defined to be symmetric around 15-mers to compute an overall representative mean relative frequency. Taking natural random “non-binding” phosphopeptides (peptides with predicted %rank >20) as a baseline, we observe (see Supplementary Figure S3) that all three categories of phospholigands (%rank: [0,1], (1,5], (5,20]) behave similarly to this baseline, indicating that there is no particular preference with respect to the phosphosite localization between the binding core and PFRs.

Following this, we studied if any of the 9 binding core positions were enriched or depleted in phosphorylated amino acids. HLA-II binding predictions were computed for 200,000 random natural phosphopeptides for HLA molecules with 100 or more phospholigands assigned during model cross-validation (in total 50 HLAs). For each of the selected HLAs, the top 2.5% scoring predictions were taken as binders (top 5000 peptides per HLA), and an equal number of randomly sampled natural phosphopeptides with a predicted % rank >20 were considered non-binders.

Figure 7 shows the log2 ratio between the relative frequency of binders vs non-binders for each of the 9 binding core positions, discriminating by HLA loci/subtype and phosphorylation type. Overall, there are very limited consistent phosphosite enrichment/depletion signals. One exception being a few DRB1 alleles that are suggested to have enriched preference for pY at P1. In terms of depletion, most such observations appear to be located at the conventional class II anchor positions. Key examples of such marked HLA anchor position depletions include both pS and pT at P1 for almost all DRB1, DRB3/4/5, and HLA-DP molecules, and pY at P1 for most HLA-DP molecules, pY at P6 and P9 for DRB1 molecules, pY at P9 for HLA-DP and pY at P6 and P9 for HLA-DQ molecules.

Overall, these observations, together with the HLA binding motif analyses discussed earlier, suggest that phosphosite location is to a high degree aligned with the preferences described for unmodified S/T/Y residues. This is further supported by Supplementary Figure S9, which compares the log₂ binder/non-binder ratios derived from random natural phosphopeptides with those obtained from random natural unmodified peptides across all nine core positions, stratified by HLA-II loci and residue type. When considering phosphosite depletion (log₂ ratio ≤ −1), a strong concordance with unmodified peptides is observed, with 92% (192/208) of depleted phosphosites matching depleted S/T/Y positions in unmodified peptides, of which 91% (175/192) are located at canonical HLA-II anchor positions (P1, P4, P6, and P9). These results indicate that phosphorylated S/T/Y residues are disfavored at certain anchor positions largely since their unmodified counterparts are themselves not preferred at these same positions. In contrast, among allele–position combinations showing phosphosite enrichment (log₂ ratio ≥ 1), only 58% (7/12; with no overlap observed for HLA-DP or HLA-DQ) coincide with core positions enriched for unmodified S/T/Y residues, and all such coincident cases occur at anchor positions. Together, these observations indicate that, within the current prediction framework, the impact of phosphorylation on HLA-II binding is predominantly constrained by the underlying enrichment/depletion patterns defined for the corresponding unmodified residues.

To further examine the impact of phosphorylation on peptide binding, we compared NetMHCIIphosPan cross-validation predictions for phospholigands in the training data and their unmodified counterparts (Supplementary Figure S10). As in the earlier analyses, phospholigands with a predicted %rank >20 were excluded as likely HLA-irrelevant coimmunoprecipitated contaminants. Among phosphopeptides that also appeared in unmodified form in the training data set (10% of total), 95% retained the same predicted binding core, while only 5% were assigned an alternative one. For the remaining 90% without matched unmodified forms, 91% conserved their predicted binding core, and 9% exhibited a core shift. Thus, in the vast majority of phospholigands (91%), the predicted binding core remains conserved upon dephosphorylation, and for this subset, changes in binding likelihood were generally small, as indicated by the narrow distribution of delta %rank values centered around zero (see Supplementary Figure S10). In contrast, phospholigands with a shifted binding core exhibited delta %rank values of larger magnitude, although still modest in most cases. In both groups, a prevalent trend toward increased %rank upon dephosphorylation was observed, reflected by the density shift toward positive values in the delta %rank distributions, indicating that removal of the phosphate group tends to reduce predicted binding likelihood. Focusing on the most relevant subset of phospholigands, those predicted as binders by our method (% rank ≤5), which account for 74% of the training data set (14,046 out of 19,019), only 2.2% exceeded a %rank of 10, and just 0.37% surpassed 20 upon dephosphorylation.

In summary, these results indicate that phosphorylation does not substantially alter either the predicted 9-mer binding core or the binding likelihood, reinforcing our previous observation that phosphosite localization preferences coincide with those of unmodified S/T/Y residues. While dephosphorylation may moderately weaken phosphopeptide–HLA-II interactions in some cases, according to the developed model, the overall impact is limited.

A Comment on Phosphosite Localization Certainty

Our analyses thus far have not accounted for the ambiguity in phosphosite localization, as indicated by the AScore (or Ambiguity Score) parameter computed by PEAKS. A more detailed examination of this aspect is provided in the Supplementary Note. Briefly, in this examination, we observed that model performance does not improve when exclusively trained on high-confidence site data (ligands with AScore >13), suggesting limited sensitivity to localization certainty as annotated by PEAKS. Additionally, attempts to increase the accuracy of phosphosite localization annotation either by using phosphosite prediction tools, such as NetPhosPan⁴³ or PhosphoLingo,⁴⁴ or by overlapping with public phosphoproteomics databases demonstrated, in all cases, a very limited recovery (i.e., only a small fraction of the observed phosphosites in the identified phospholigands were either predicted to be phosphorylated or matched to entries in the phosphoproteomics databases). In the case of phosphosite predictions, sequence composition analysis revealed biases in phosphosite prediction tools, which tend to favor cytosolic or nuclear kinase motifs. Therefore, these tools are of limited utility for our purposes, as most of the identified HLA-II phospholigands are derived from the endolysosomal pathway (as predicted by DeepLoc-2.0⁴⁵). Together, these findings reveal that current phosphorylation prediction tools and databases are restricted by biases and specific data set characteristics, emphasizing the need for tailored approaches to achieve accurate phosphosite localization assignment.

NetMHCIIphosPan Predictions on Phosphopeptides Identified by FragPipe

To further evaluate the generalizability of NetMHCIIphosPan and investigate the limited overlap between phospholigands identified by PEAKS and MaxQuant, we reanalyzed the Racle and Racle_DR data sets using FragPipe (v23.1, MSFragger v4.3, Philosopher v5.1.2;^18,28,29 for more details refer to Methods). These data sets were selected for analysis as they yielded the highest number of phosphopeptide identifications among all samples, as previously shown in Figure 1. We further reran the PEAKS analyses with adjusted search parameters, including a 1% FDR and a reduced number of variable PTMs, to allow for a better comparison between the two analyses (for details refer to Methods). Comparing the sets of unmodified peptides identified by the two methods demonstrates a very high concordance between the two approaches both in terms of peptide yield (242,247 and 261,781 unique peptides identified by PEAKS and FragPipe, respectively) and overlap of 86% (with respect to the minority set; see Supplementary Figure S11A). However, when focusing on the phosphopeptides, the concordance between the two methods is greatly reduced. Here, the overlap is of 29% (with respect to the minority set), and the overall peptide yield is substantially lower when considering FragPipe compared to PEAKS, with 3320 and 5407 unique phosphorylated peptides identified by FragPipe and PEAKS, respectively, see Supplementary Figure S11B. When phosphorylated residues are reverted to their unmodified forms, thereby disregarding both the number of phosphorylations and the exact phosphorylation sites when multiple S/T/Y acceptor sites are present within the same sequence, only a modest increase in overlap is observed (Supplementary Figure S11C).

To further assess the accuracy of the phospholigands identified by the two methods, we stratified the peptides into two subsets: peptides identified by both methods (“overlapping”), or peptides identified only by one of the two methods (“non-overlapping”), and next extracted prediction scores for the peptides using NetMHCIIphosPan. For peptides identified exclusively by FragPipe, predictions were generated as an external evaluation set, since these peptides were not included in model training. For peptides identified by PEAKS, which were used in model training, predictions were obtained through 5-fold cross-validation. As shown in Supplementary Figure S11D, overlapping phosphopeptides displayed consistently high predicted presentation likelihoods, with close to 77% of the peptides scoring at %rank ≤1. Among non-overlapping peptides, both methods display a comparable performance with only 6–8% of the peptides being assigned a rank >20%. In terms of rank values, PEAKS displays a small advantage over FragPipe, with 75% of the peptides being predicted with a rank ≤5%. In comparison, this value is 67% for FragPipe. Returning to the peptides identified using PEAKS with the search criteria applied to extract the training data for NetMHCIIphosPan (among other things relaxed to an FDR of 5% and an increased number of PTMs (refer to Methods and Figure 1)), we show in Supplementary Figure S11E, the distribution of percentile rank values for the two peptide subsets; the subset identified using the search criteria with the FDR of 1%, and the subset only identified using the relaxed settings. The results of this demonstrate a very similar %rank distribution of the two subsets, with less than 5% of the peptides in both data sets being assigned a rank >20%, and more than 80% of the peptides being assigned a rank value less than 5%.

Together these results suggest that despite a rather low concordance between the identified sets of phosphorylated peptides between the FragPipe and PEAKS methods, the HLA binding patterns contained within the different data sets is consistent, allowing a method trained on the PEAKS data to accurately predict peptides identified by FragPipe. Further, the similarities between the percentile rank values for the peptides identified using an FDR of 1% or 5%, supports that the impact of the increased false discovery imposed by the relaxed FDR rate is very minor if at all present.

Benchmarking against Other Earlier Methods: MixMHC2pred-1.3

As a final test, the NetMHCIIphosPan method was benchmarked against MixMHC2pred-1.3²⁴ employing an independent immunopeptidome data set of 2058 phospholigands obtained from esophageal adenocarcinoma (EAC) samples (11 cell lines) and monocyte-derived dendritic cells pulsed with a recombinant SARS-CoV-2 Spike protein vaccine candidate (5 cell lines). Note that the overall peptide yield from the EAC samples was poor given their low HLA-II expression (for details on this data set refer to Materials and Methods). As explained before, we evaluated the performance with the ROC-AUC, ROC-AUC0.1, PPV and PR-AUC metrics. Note that the external data sets only cover DRB1 and a few DRB3/4/5 alleles. Also, it is important to stress that, in contrast to NetMHCIIphosPan, MixMHC2pred-1.3 can only make predictions for alleles present in its training data set. To accommodate for this limitation, the benchmark was limited to consider only the HLAs available in MixMHC2pred-1.3. Of note, when evaluating performance using this external benchmark, we ensured no overlap between the training data set of our model (or NetMHCIIpan-4.3) and the independent evaluation data (for more details refer to Materials and Methods). In contrast, MixMHC2pred-1.3 was evaluated without controlling for such overlap, potentially resulting in inflated performance estimates due to data leakage.

Figure 8 presents the benchmark results, showing that NetMHCIIphosPan outperformed MixMHC2pred-1.3 across all metrics for the two external data sets (all P < 0.0001; for performance evaluation details refer to Methods section). Remarkably, performance values for the EAC data set remain low, particularly for ROC-AUC0.1, PR-AUC and PPV, compared to the cross-validation results (Figures 2 and 4), suggesting a lower quality for this data set compared to the training data sets. Supplementary Figure S12 extends the comparison shown in Figure 8 by including prediction results from NetMHCIIpan-4.3, evaluated on three “phospho-equivalent” data sets. Additionally, Supplementary Figure S13 reports prediction performance across the full allele set, excluding MixMHC2pred-1.3, which supports only a limited set of alleles, as previously noted.

The NetMHCIIphosPan method is made available as a web server and stand-alone tool at https://services.healthtech.dtu.dk/services/NetMHCIIphosPan-1.0. The tool has functionality like the other NetMHC-tools and allows for the submission of either phosphopeptides or proteins with annotated phosphosites, making predictions for any HLA class II molecule with a known sequence (for details refer to the web server instructions). The web server implementation returns predictions based on the full ensemble of 100 networks constituting the final NetMHCIIphosPan predictor.