More than 1,900,000 published HLA-peptides MS data of mono-allelic or mixed-allelic cell lines which collectively expressed 20 high-frequency HLA-A and HLA-B allotypes are collected (16, 19, 29, 30). All these data are labeled in binary notation. Label=1 denotes MS-identified peptides (hits), whereas label=0 denotes peptides from the reference proteome (SwissProt) that were not detected via mass spectrometry.

The training datasets are peptides with the length from 8-mer to 11-mer, which are represented by a one-letter amino acid alphabet (a total of 20 distinct amino acids, namely ‘ACDEFGHIKLMNPQRSTVWY’). Such length range captures ~95% of all HLA class I-restricted peptides. To implement machine learning, the peptide sequences are vectorized by a one-hot encoding scheme. Peptides with multiple lengths (8-mer to 11-mer) were represented as fixed-length vectors by using a padded character ‘Z’. Each amino acid and the padded ‘Z’ are encoded as a one-hot vector (see Figure S1 for details). As a result, peptides are encoded as the fixed matrix of 11 rows (maximum length) by 21 columns (20 distinct amino acid alphabets and the padded character ‘Z’).

The collection of MS datasets shows a severe class imbalance. Overall, the total number of 0-labeled data is 1,866,484 which is 39 times as many as the 1-labeled counterparts. An extreme case can be found in HLA-A*02:07 datasets where the negative-labeled records are 72 times more than 1-labeled records. Such extreme imbalance influences the prediction of the machine learning model, inclined to show a better performance on the 0-labeled peptides (the majority) and a worse on the 1-labeled ones (the minority) (31). Thus, the class balance is adjusted via over-sampling and under-sampling procedures in preprocessing the training datasets. Briefly speaking, the under-sampling goes by removing the 0-labeled training data points at random, whereas the over-sampling duplicates the 1-labeled data points. Table 1 shows the proportions of over-sampling and under-sampling on different HLA alleles.

Usually, the Convolutional Neural Network (CNN) consists of convolutional layers, pooling layers and fully connected (dense) layers. In this study, an advanced CNN which is inspired by the inception module from GoogLeNet is used (32, 33). Three parallel convolutional sections with eight two-dimensional convolutional kernels for each were constructed to maximize the feature extraction (see Figure S2 for details). The output of three convolutional layers connects to a flattened matrix and is delivered to the fully-connected layers which contain 100 hidden nodes. The output layer displays the results of binary classification by two nodes where a tested peptide is classified as binding or not binding to HLA.

The model is implemented with Tensorflow (v. 1.14.0) and trained by Adam optimization algorithm with standard parameters on an NVIDIA GeForce RTX 2080 Ti GPU. Instead of the frequently-used activation function Rectified Linear Unit (ReLU), the advance function of Leaky ReLU (α=0.2) is applied to activate the model and the “drop-out” and “early stopping” schemes are introduced to avoid overfitting.

The peptides of the MS dataset are randomly split into training sets, validation sets and test sets, and all three sets have approximately the same distribution of 1-labeled and 0-labeled peptides. The validation sets are used only for early stopping. The training sets are used to perform feed-forward and backpropagation and the test sets are used to evaluate performance via AUC.

To benchmark the APPM and other HLA-peptide predictors, we collected HLA-bound peptides MS datasets from other studies that use cell lines to express a single HLA allele (34, 35). From these MS-identified peptides (hits), we generated non-binders (decoy sets) by sampling unobserved peptides from the same proteins through the Uniprot human reference proteome (UP000005640_9606) as previously described (36). For each MS-identified peptide, we randomly selected 99-time decoy peptides of four different lengths (8, 9, 10, 11), and the number of each length is the same. The rationale for the 99-fold bias is that for a sample of peptide fragments from an organism, it is commonly considered that approximately 1%∼2% of the fragments will bind to MHC receptors (37). After removing the peptides appearing in the model training data and the duplicate sampled from different proteins, we obtained a mono-allelic benchmark dataset.

Sensitivity, also called recall, was calculated as:

Specificity was calculated as:

Positive predictive value, also called precision, was calculated as:

To obtain a consensus driver mutations list, we download the driver-mutations dataset processed and compiled by TCGA MC3 and driver working group (https://gdc.cancer.gov/about-data/publications/pancan-driver) (38, 39). The driver-discovery dataset was derived from a compiled MAF file of 9079 TCGA samples across 33 different cancer types (syn7824274, https://gdc.cancer.gov/about-data/publications/mc3-2017). Based on sequencing and structure analyses, we ultimately selected 3,437 cancer driver mutations as the consensus list were identified by ≥ 2 approaches from CTAT-population, CTAT-cancer, or structural clustering (see Supplementary File 4 ).

For each driver mutation, we extract 8-11mers candidate peptides that contain the driver specific mutant amino acid for neoantigen screening. For instance, the extracting procedure of 9-mer candidate peptides is described as follows ( Figure S3 ). Firstly, we extracted a 17-mer peptide from the protein sequences, where the mutant amino acid was placed in the center with eight upstream and downstream wild amino acids as flanks. Secondly, by using the sliding window protocol, a 9 amino acid size window was slid N (N = 9) times to obtain 9-mer peptides. Briefly speaking, the mutant amino acid serves as the end point of the first 9-mer peptide. This 9-mer sliding window moves along the 17-mer fragment until the mutated point becomes the starting point of the 9-mer. Peptides with other lengths are treated in the same way.

We aimed to improve the precision and specificity of the HLA-peptide prediction approaches through a novel tool that has been trained on improved training data and a new supervised machine learning model. HLA-Peptides of MS data were eluted by immunoprecipitation of HLA molecules and then identified by liquid chromatography-tandem mass spectrometry (LC-MS/MS) (40, 41). Compared with in vitro binding affinity assays, MS data directly profiles peptides that are actively presented by cells or tissues (42). We collected publicly available HLA-peptides MS data from 16 mono-allelic HLA-A and HLA-B cell lines genetically engineered to express a single HLA allele and from B lymphocytes or cancer cell lines expressing multiple HLA complex alleles (16, 19, 29, 30, 43). These MS data consist of 20 high-frequency HLA-I alleles. We split the datasets into three sets: training, validation and testing sets (Methods). Owing to so many negative peptides (from reference proteome), we apply the over-sampling and under-sampling scheme, which neutralizes the substantial fraction of the imbalance issue.

Using these public HLA-peptides MS data, we build a convolutional neural network (CNN) framework to predict HLA-I presentation, a form of deep learning that excels at handling general sequence data such as amino acid sequences ( Figure 1 ) (28). The model has three parallel convolutional modules, each consisting of eight two-dimensional convolutional layers, which preserved HLA class I-peptide binding features.

To estimate the predictive performance of APPM, we first compared the prediction results of APPM with the IEDB recommended method (2020.04) (NetMHCpan4 EL (44), the state-of-the-art class I binding predictors available at http://tools.iedb.org/mhci/) in terms of PPV. We compiled a benchmark using published MS data from cell lines genetically engineered to express a single HLA-I allele. In this mono-allelic benchmark, the MS-identified peptides are true positives where length-matched amino acid fragments from the same protein as negative peptides (decoys). For each paired HLA allele and peptide, NetMHCpan4 EL produced a binding score and percentile ranks. Using the recommended threshold of the percentile rank (top 2% ranks are considered binders), we obtained the average specificity and positive predictive value (PPV) of 0.97 and 0.22 for NetMHCpan4 EL ( Supplementary File 1 ).

When tested on the same data, APPM outperformed NetMHCpan4 EL with the specificity of 0.99 and PPV of 0.40. The improvement in reducing false positives rates was substantial, with an average of 80% increase in PPV ( Figure 2A ). For the 20 frequent haplotypes of HLA class I, APPM only exhibited a slightly lower PPV than NetMHCpan4 EL on HLA-A*02:01, but presented higher PPV for the rest of 19 HLA haplotypes, particularly with more than one fold of increase for HLA-A*02:03, HLA-A*29:02, HLA-A*32:01 and HLA-B*40:01 ( Figure 2B ), suggesting the advantage of our algorithm.

Interestingly, a low overlap rate (19%) is observed between APPM and NetMHCpan4 EL for the false-positive peptides ( Figure 3A ), probably due to the different prediction mechanisms. In this case, we hypothesized that the prediction performance could be improved by combining these two predictive approaches. We redefined the predictive results: only peptides identified positively in both methods are regarded as positives. Using the combined predictions, we obtained the PPV of 0.51 ( Figure 3B ), which is significantly higher than that of both APPM and NetMHCpan4 EL ( Figure 3C , p = 0.013, t-test and Figure 3D , p < 0.001, t-test), without significant decrease of sensitivity ( Figure 3E , p = 0.1, ANOVA). These results suggested that the combined predictions from different algorithms can improve the positive rate for neoantigen selection, which is consistent with previous studies (45, 46).

To illustrate the binding characteristics of HLA-I alleles with peptides, we draw allele-specific presentation motifs for 20 HLA-I alleles (see Supplementary File 2 for motifs of all alleles). Consistent with previous studies (17, 19, 47), these motifs revealed the dependence of HLA presentation on each sequence position for peptides of multiple lengths 8-11 ( Figure 4A ). For example, the anchor residues of 9mer are amino acid at position 2 (refer as P2, a similar abbreviation for other positions) and P9, while 11mer at P2 and P11.

In contrast to previous work (48), some distinct HLA alleles have similar presentation motifs. For instance, HLA-A*02:01 and HLA-A*02:03 have the same binding specificity, meaning the pockets preferentially bind to bind the peptides with leucine at P2 and valine/leucine at the last position. Likewise, HLA-A*03:01 and HLA-A*11:01 presented lysine at the last position, while HLA-B*40:01, HLA-B*44:02, and HLA-B*44:03 prefer to deliver peptides with glutamate at P2 ( Figure 4B ).

Moreover, we analyzed the amino acid properties of anchor residues of 20 HLA alleles and refined their binding character: these binding peptides enriched in hydrophobic amino acids at anchor residues. It is consistent with the known preference of HLA-I binding and presentation (23, 49). We also explored the whole preference of amino acid properties among HLA-A and HLA-B molecules on anchor residues ( Figure 4C ). Besides the common preference of hydrophobic amino acids, HLA-A alleles prefer to bind basic and polar amino acids, while the HLA-B alleles prefer acidic amino acids.

It is considered that the quality rather than the quantity of neoantigens may lead to a robust and durable response to immunotherapy (50). Most of the putative neoantigens are considered as the product of passenger rather than driver mutations, and their loss through chromosomal instability during tumor evolution may be readily tolerated. Therefore, targeting driver-mutation-neoantigens could manifest durable anti-tumor responses and may reduce the resistance to neoantigen therapies.

We applied the combining approach of APPM and NetMHCpan4 to predict neoantigens derived from oncogenic driver mutations. The consensus driver-mutation list was compiled and discovered by The Cancer Genome Atlas (TCGA) Multi-Center Mutation Calling in Multiple Cancers (MC3) working group and driver working group among 9079 samples across 33 cancer types (38, 39). For a total of 3,437 missense driver mutations, we identified ~ 16,000 putative neoantigens in the context of 20 high-frequency HLA alleles ( Supplementary File 3 ).

Among these driver mutations, only 15% (513/3437) do not yield putative neoantigens, while the products of the other could be bound and presented by these HLA alleles. We identified 36 high-frequent shared putative neoantigens derived from eight oncogenic driver mutations with more than 1% coverage of multiple cancer patients in the 9079 TCGA cohort ( Table S1 ), e.g. HLA-A*03:01_KIGDFGLATEK from BRAF_p.V600E with 5.60% (508/9079) in Pan-Cancer. Besides, we also found tumor-specific shared potential neoantigens with over 10% frequency in a given cancer type ( Table S2 ). For example, HLA-B*15:01_IIIGCHAY from IDH1_p.R132C with 11.76% (4/34) in CHOL. Importantly, the immunogenicity of some shared putative neoantigens we identified has been confirmed experimentally ( Table 2 ) (51). For instance, VVVGAGDVGK from KRAS_p.G13D has been shown to be immunogenic in the context of the HLA-A*03:01 allele. Overall, these putative shared driver-mutation-neoantigen pools provide a potential list of targets for off-the-shelf immunotherapy.