Data retrieval of MHC class I binding and nonbinding peptides

A dataset of MHC class I binding and nonbinding peptides, with reduced similarity, was compiled for 22 molecules from the supplementary data provided by the MHCIPREDS web server, which is accessible at RED (https://web.archive.org/web/20130828192234/http://ailab.cs.iastate.edu/red/). The log-transformed MHC BA (LTMBA) values between 0 (no affinity) and 1 (very high affinity) were obtained using the relation, 1 – log (aff)/log (50,000), where aff is the experimentally measured BA in terms of half-maximal inhibitory concentration (IC₅₀) with nM (nmol/L) unit. An LTMBA threshold value of 0.426, equivalent to an IC₅₀ value of 500 nM, was used to classify binders and non-binders, which means peptides with LTMBA values ≥ 0.426 (IC₅₀ ≤ 500 nM) were classified as binders (Table 1).

Similar peptides were removed to generate a similarity-reduced cross-validation dataset for each MHC allele, and the data were randomly partitioned. The present study used a 5-fold cross-validation approach, i.e., the dataset was split into five subsets. Initially, the dataset was scanned for equal binders and non-binders based on LTMBA value ≥ 0.426 for each 22 MHC class I allele. Then, four out of five (4/5) of each dataset were used for the training, and 1/5 data points (peptides) were used for the blind evaluation (because none of the peptides in the evaluation set were included in the training set at any stage) as described by Nielsen et al. [31] (Supplementary file 1 in https://data.mendeley.com/datasets/dxz3dk3tcm/1).

Construction of MHC class I binding motifs (PSSM)

Initially, based on a training dataset of fixed length (9-mer) MHC class I binding peptide sequences (positive) and using a statistical model of sequence weighting methods available at EasyPred modeler (https://services.healthtech.dtu.dk/services/EasyPred-1.0/), 22 PSSM (9 × 20) were constructed, which represents the frequencies of residues observed for a position in a multiple sequence alignment (MSA) assuming no correlations exist between the different peptide positions [32]. The score S of a peptide to a motif is usually calculated as the sum of the log-odds ratio.

Where p_pa is the probability of finding amino acid a (a can be any of the 20 amino acids) at position p (p can be 1 to 9) in the motif, and q_a is the background frequency of amino acid a, log_k is the logarithm with base k. The scores are often normalized to half-bit by multiplying all scores by 2/log_k (equation 2). In half-bit units, the log-odds score S is calculated as:

The three sequence weighting methods, namely Henikoff & Henikoff 1/nr [26], Hobohm [25] clustering at 62% identity, and no clustering, were used to compensate for the over-representation among MSA, along with four different weights on pseudo counts (50, 100, 150, and 200) [33–35]. In the Henikoff & Henikoff 1/nr method [26], an amino acid a on position p in sequence k contributes a weight w_kp = 1/nr, where n is the number of different amino acids at a given position (column) in the alignment and r is the number of occurrences of amino acid a in that column. The weight of a sequence is then assigned as the sum of the weights over all positions in the alignment. However, in the Hobohm 62% clustering method [25], each peptide k in a cluster is assigned a weight w_k = 1/n_c, where n_c is the number of sequences in the cluster containing peptide k. When the amino acid frequencies are calculated, each amino acid in sequence k is weighted by w_k. The Henikoff & Henikoff method is as fast as the computation time increases linearly with the number of sequences. In contrast, in the Hobohm clustering algorithm, computation time increases as the square of the number of sequences. Additionally, the binding potential (score) of any peptide sequence (query) to a specific MHC allele was determined by aligning the corresponding PSSM with the peptide sequence and summing the scores that correspond to the residue type and position in the PSSM. To narrow down the probable binders from the list of scored and ranked peptides, a binding cut-off score was established that encompasses 85% of the peptide sequences in the training dataset (positive data) [28, 36].

Training of ANN for the forecast of MHC class I binding peptides

An ANN is a computer simulation of a system of interconnected processing units that can be trained to extract and remember a pattern present within a data set. It can subsequently recognize that pattern when presented with new data. In binding a peptide to the MHC molecule, the amino acids might compete for the space available in the binding groove. Therefore, the mutual information in the binding motif will allow the identification of higher-order sequence correlations. Neural networks with hidden layers are used to describe sequence patterns with such higher-order correlations (https://teaching.healthtech.dtu.dk/mpmbioinformatics/22801_04.pdf). Thus, the ANN simulations were performed for all 22 MHC class I molecules using the neural network options available in the EasyPred modeler [27]. Based on the other studies of MHC-binding peptides forecast, two ANN architectures, namely 180-1-1 and 180-2-1, were considered in the current study, which represent respective numbers of neurons in the input, hidden, and output layers, respectively [37, 38]. In these architectures, the input 9-mer peptide sequences were presented with 180 values (where individual amino acid was encoded as a binary string of length 20 with an exclusive position set to 0.9 and other positions set to 0.05). At the same time, the output from each neuron was transformed using the standard sigmoid function [39]. The training algorithm employed to generate the final network was the steepest descent method that learns from a training set of input-output pairs by modifying the network weight parameters such that the network generates a numerical value for each input close to the desired target output. Throughout the training process, the LTMBA value of a peptide served as a target value and was deliberated as a strong binder and a potent epitope for MHC class I alleles [40]. Nevertheless, the learning process utilized was error backpropagation [41]. In addition, the following parameters were used for ANN simulations: i) The top 80% of the training set was used to train the neural network, and the bottom 20% was used to stop the training to avoid fitting; ii) The learning rate was 0.005; iii) The maximum number of iterations was set to 300.

Evaluation parameters

The nonparametric performance measure, the area under the receiver operating characteristic curve (Aroc) and correlation coefficient (CC) values, was used to evaluate the predictive performance of the EasyPred modeler-based PSSM and ANN predictors. The receiver operator characteristic (ROC) curve is a plot of the true positive rate [TP/(TP + FN)] on the Y-axis versus the false positive rate [FP/(TN + FP)] on the X-axis for the entire range of the decision thresholds where, true positive (TP) is an experimentally proven binding peptide forecasted as a binder, false positive (FP) is an experimentally proven nonbinding peptide forecasted as a binder, true negative (TN) is an experimentally proven nonbinding peptide forecasted as a non-binder and false negative (FN) is an experimentally proven binding peptide forecasted as nonbinder [42, 43]. An Aroc value can be interpreted as the probability of distinguishing a true positive from a false positive. For the calculation of Aroc, peptides were classified into binders and non-binders at a cutoff value of 500 nM. This affinity threshold is associated with the most well-known T-cell epitopes. The values Aroc ≥ 0.90 indicate excellent, 0.90 > Aroc ≥ 0.80 good, 0.80 > Aroc ≥ 0.70 marginal, and Aroc < 0.70 poor predictions [42]. The CC is another widely used measure of the association between pairs of values (predicted versus experimental). It is calculated as:

where the value p_i is found using a prediction method of choice, and the a_i is the known corresponding target value. However, the overlined letters denote average values. The value of 1 corresponds to a perfect correlation and –1 to a perfect anti-correlation, and 0 value corresponds to a random prediction [42].

Generalization test and Epstein-Barr virus case study

The generalizability test was used to assess the ability of the PSSM and ANN models, trained for the allele HLA-A*0201, to generalize to the other 21 alleles. Training was performed on the same dataset, followed by testing on evaluation datasets of all 22 alleles. In addition to the above, the performance of the PSSM and ANN predictors was validated through observations revealed in the very recent study conducted by Wohlwend et al. [44] called the Epstein-Barr virus (EBV) dataset. This research indicates that the dataset comprises 11 distinct immunogenic epitopes that are restricted by 5 HLA class I alleles derived from EBV, as identified through IFNγ ELISpot assays. The study employs PSSM and ANN models to effectively identify both confirmed and new HLA class I epitopes from EBV, which have been experimentally validated through in vitro and ex vivo studies.

Prediction and relative rank measurement of PSSM and ANN predictors on NetCTLpan and SARS-CoV-2 dataset

The additional performance evaluation of the present EasyPred modeler-based PSSM and ANN predictors, along with recently published NetMHCpan 4.0 methods [eluted ligands (EL) and BA] (https://services.healthtech.dtu.dk/services/NetMHCpan-4.0/), was estimated on NetCTLpan and SARS-CoV-2 dataset in terms of comparative ranking of the reported ligands among all nonamers included in the source protein as described by Larsen et al. [45, 46]. This NetCTLpan dataset was gathered from the supplementary material of the NetCTLpan1.0 method (https://services.healthtech.dtu.dk/suppl/immunology/NetCTLpan.php) reported as SYFPEITHI 9-mer training and evaluation dataset [47, 48]. The dataset involved 413 sequences of naturally processed T cell epitopes restricted by 7 HLA class I molecules that were common between the present study and the NetCTLpan1.0 tool (Supplementary file 2 in https://data.mendeley.com/datasets/dxz3dk3tcm/1). At the time of conducting the present study, initially, no T cell epitope data were available for SARS-CoV-2, but there was a significant amount of information available on T cell epitopes for betacoronaviruses that cause similar diseases in humans, like SARS-CoV. Thus, we have compiled the SARS-CoV-2 dataset from the study of Grifoni et al. [49] that contains 11 T cell epitopes (9-mer) derived from surface glycoprotein (NCBI ID: QHD43416) and nucleocapsid phosphoprotein (NCBI ID: QHD43423) that showed 100% identity with SARS-CoV epitopes. Additionally, evaluation of PSSM and ANN predictors on the recent SARS-CoV-2 dataset of Gfeller et al. [50] was performed for prediction and % rank analysis of CD8+ T cell epitope in their source protein using allele-specific cutoff score.

Integrated forecast of CTL epitope processing

The CTL epitope processing forecast integrated the three significant forecast steps in the filtering approach, including MHC class I BA, TAP transport efficiency, and C-terminal proteasomal cleavage (Figure 1). The details are described below:

Display full size

Figure 1. Flow chart indicating the integrated forecast of cytotoxic T lymphocytes (CTL) epitopes based on the EasyPred modeler. MHC: major histocompatibility complex; TAP: transporters associated with antigen processing.

• 1.

  Amino acid sequences of a given target protein were parsed into 9-mer overlapping peptides. Peptides were then scored using EasyPred-based PSSM and ANN predictors, and values above the defined threshold (Table 2) for each 22 MHC class I allele were forecasted as binders. In the case of PSSM predictors, we calculated the threshold scores in terms of predicting peptides that bind with 85% of all epitopes in the training set because an established threshold associated with immunogenicity (i.e., IC₅₀ ≤ 500 nM) covers 80–90% of all immunogenic epitopes [51]. However, in the case of ANN predictors, the threshold value of an IC₅₀ less than 500 nM was considered for binder forecast.

• 2.

  Forecasted MHC class I binders were again scored by the quantitative weight matrix (9 × 20) for TAP BA in terms of –log IC₅₀ (pIC₅₀) described by Peters et al. [52] (Table S1) and/or Doytchinova et al. [53] (Table S2).

• 3.

  Peptides with scores less than the selected threshold value (default, pIC₅₀ < 4) were considered as TAP binders. The TAP transports peptides into the ER that are potentially N-terminally extended from the ligand that ends up in MHC. This means that the peptide that binds to MHC does not necessarily need to be a suitable substrate for TAP [52].

• 4.

  Subsequently, the proteasomal (constitutive or immuno) cleavage score of 12-mer overlapping peptide fragments generated from the target protein sequence was calculated along with recording of C-terminal position at the cleavage site (i.e., position of the sixth amino acid in the target protein) by using a quantitative cleavage weight matrix (12 × 20) published by Toes et al. [54] (Table S3 and S4).

• 5.

  The query 12-mer peptide with a score above a chosen threshold value (corresponding to 1–10%) was forecasted to be cleaved between the sixth and seventh amino acids [55].

• 6.

  Finally, MHC class I and TAP binding peptides identified through the above steps involving proteasomal (constitutive or immuno) cleavage site position at their C-terminal (steps 1–5) were considered CTL epitopes.

Construction and evaluation of PSSM and ANN predictors

Sequence weighting methods have been employed in creating a PSSM from frequencies of residues surveyed for a position in an MSA of MHC binders [32]. The present study utilizes three sequence weighting methods (Henikoff & Henikoff 1/nr, clustering at 62% identity, and no clustering) available at the EasyPred modeler (https://services.healthtech.dtu.dk/services/EasyPred-1.0/) for the construction of MHC allele-specific PSSM, including four corrections for low counts (weight on pseudo counts) 50, 100, 150, and 200 [34, 35]. As peptides in the evaluation dataset are not included in the training dataset, it is equivalent to a blind test [27]. For each MHC class I allele, only those PSSM were selected as predictors that gave maximum predictive performance evaluated in terms of Aroc and CC values (Table S5). The Aroc value operates independently of the anticipated scale, as it assesses the ranking of predictors and remains unaffected by the dataset’s composition, including varying ratios of binders and non-binders. The Aroc value provides a crucial metric for evaluating forecast quality, with a score of 0.5 indicating random forecasts and 1.0 representing perfect forecasts. Conversely, the CC value of one corresponds to perfect correlation, a value of zero indicates a random estimate, and a value of minus one signifies perfect anti-correlation. The Aroc value principally captures the probability that, given two peptides, one a binder and the other a non-binder, the forecasted score will be higher for the binder than the non-binder [43]. Furthermore, the MSA of each HLA class I allele binding motif was visualized by using the graphical depiction (sequence logo) method [56], where the height of a column of letters represents the information content (I) at that position. The height of each letter within a column is proportional to the frequency of the corresponding amino acid at that position and colored according to their physicochemical properties, such as acidic (DE)-red, basic (HKR)-blue, hydrophobic (ACFILMPVW)-black, and neutral (GNQSTY)-green (Figure S1). The similarity in peptide binding preferences is observed when comparing the MHC class I binding motif logo developed in the presented study with the logo stored in the MHC motif viewer. In the MHC motif viewer database, pairwise assessments of MHC binding motifs facilitate immediate analysis of epitope selection data in patient cohorts with HLA diversity [57]. For example, when comparing the binding motif of human (HLA-A*2402) to the chimpanzee (Patr-A*0701) allele, an unmistakable resemblance between the binding motifs of the two alleles was spotted, as noted by Sidney et al. [58]. A similar conserved motif for MHC class I peptide binding has also been shown between humans and rhesus macaques, as demonstrated by Dzuris et al. [59]. Consequently, the ability to differentiate various MHC binding specificities (motifs) has applications that range from the design of experiments for peptide binding assays to personalized medicine for a significant population, which includes the selection of peptides that are immunogenic in both humans and model organisms [60, 61].

Moreover, in the case of ANN simulations, two architectures (180-2-1 and 180-1-1) were initially used, and the input peptide sequences (binders with IC₅₀ ≤ 500 nM and rest non-binders) were offered in a conventional encoding, as described by Nielsen et al. [27]. Further, the network weights were renewed using gradient descent backpropagation algorithms. For a given peptide sequence of 9-mer, the ANN weights were renewed to lower the sum of squared errors between the forecasted and experimentally measured BA (target value). The training of the neural networks was performed using a five-fold cross-validation (Supplementary file 1 in https://data.mendeley.com/datasets/dxz3dk3tcm/1). For each of the five training and test subsets, a series of network training sessions is conducted, each utilizing two distinct hidden neurons (1 and 2) along with a single bin for balanced training. For every series, a single network exhibiting the highest test performance (measured by the highest CC and the lowest square error) was ultimately chosen as an ANN predictor (refer to Table S6). Subsequently, the predictive performance of the corresponding ANN network for each MHC allele was assessed based on Aroc (see Table S6). From these findings, it can also be deduced that increasing the number of hidden neurons within the ANN architecture does not have a significant impact on performance. Consequently, the Aroc value for the ANN predictor featuring a single hidden neuron was utilized for comparison with the PSSM predictor. It is crucial to note that, initially, a BA (IC₅₀) threshold for the MHC class I allele was established, where IC₅₀ values exceeding 500 nM were identified in competitive assays (involving an MHC and isolated peptides) and were compared with immunogenicity across various marker peptides and nonimmunogenic peptides [62]. This affinity threshold has been identified as being associated with the majority of known T-cell epitopes [63]. Nevertheless, this threshold value has also been demonstrated to be specific to MHC class I alleles [64]. An alternative way to assess the efficacy of peptide binding to MHC class I involves measuring the stability of peptide-MHC complexes over time, but affinity and stability do not rank the peptides in the same order within their source proteins and remain debatable [65, 66]. To compare the performance of ANN and PSSM predictors, the identical evaluation dataset from MHCIPREDS-IEDB was utilized. Consequently, upon comparison, the predictive performance of the ANN models was determined to be superior to that of the PSSM models for most of the MHC class I molecules (Figure 2). The ANN performance in terms of Aroc is maximum (0.97) for the alleles HLA-A*0202 and HLA-A*0203, as well as the minimum (0.56) for allele HLA-A*2902, whereas the maximum PSSM predictor performance in terms of Aroc is (0.93) for the allele HLA-A*0203 and the minimum (0.49) for allele HLA-B*4501 (Figure 2).

Display full size

Figure 2. Histogram of the predictive performance measured in terms of the Aroc value of PSSM and ANN predictors for 22 MHC class I alleles trained and evaluated on the MHCIPREDS-IEDB dataset. Aroc: area under the receiver operating characteristic curve; PSSM: position-specific scoring matrices; ANN: artificial neural networks; MHC: major histocompatibility complex.

From the above results, it is clearly stated that PSSM predictors, to a high degree, describe the binding motif of the corresponding HLA class I alleles, though they demonstrate lower Aroc than the ANN predictors. This is to be expected since the ANN can take higher-order sequence associations into account for a fixed length of MHC binding peptides [27]. Moreover, not only the size of the training and validation dataset but also the choice of specific algorithms can influence the efficiency of the forecast. To integrate PSSM-based MHC class I predictor in CTL epitope processing identification, we calculated the binding thresholds in terms of predicting peptides that bind with an IC₅₀ value less than 500 nM, an established threshold associated with immunogenicity for 85% of all epitopes in the training set [51]. These binding threshold scores for each of the 22 MHC class I molecules were determined for PSSM predictors to demarcate the range of putative binders among the top-scoring peptides (Table 2). Paul et al. [64] also revealed similar observations that different MHC molecules bind ligands at different (forecasted) binding thresholds (scores). If a single criterion for MHC class I alleles has to be deliberated, then it is preferable to select absolute BA, as in the case of the ANN predictor, where IC₅₀ ≤ 500 nM could be regarded as a reasonably good “universal” binding threshold. In a similar study, Bonsack et al. [67] confirmed the decisive threshold of IC₅₀ ≤ 500 nM for MHC class I binding peptides through in vitro validation. However, predictive efficacy is increased using allele-specific affinity thresholds.

Thus, based on these observations, the PSSM-based MHC class I binding forecasts were applied without rescaling, thereby preserving prospective fundamental biological differences between MHC class molecules. However, ANN-based MHC class I binding peptide forecasts were performed with IC₅₀ ≤ 500 nM (LTMBA value ≥ 0.426).

Evaluation of the EBV epitope dataset

The PSSM and ANN predictors successfully identified the EBV epitopes as revealed in the study conducted by Wohlwend et al. [44], above the chosen threshold (Table 2), except the epitope IACPIVMRY restricted by the HLA-B*1501 allele with a comparable BA score of 0.411 compared to the ANN predictor threshold (LTMBA value) of 0.426 (Table 3). This clearly indicates the real-life application of models in the identification of both established and novel HLA class I epitopes from EBV.

Generalizability evaluation

The generalization ability of the MHC class I binding prediction model is essential for epitope prediction, as there are many HLA alleles with inadequate data for training an allele-specific model [68]. Therefore, we have performed a detailed analysis of the performance of PSSM/ANN predictors trained on one allele and their ability to accurately predict other alleles in their evaluation datasets. Pan-specific algorithms can predict peptide binding to HLA alleles for which limited or even no experimental data are available [69]. In Figure 3, the Aroc values for the models trained on the HLA-A*0201 dataset are given for 22 different alleles. The ANN model excellently performs over the PSSM model for alleles of the HLA-A*02 type, but for the other alleles, the performance of both models is poor, except for HLA-B*4002, 4501, 5101. The prediction capabilities are good to marginal for some alleles, suggesting that cross-allele prediction is feasible in some cases. This may be due to MHC supertype classification systems that make clustered sets of HLA molecules with largely overlapping peptide repertoires. These classification systems normally depend on descriptions such as published motifs and/or analysed shared repertoires of binding peptides, etc. [70–73]. Generally, HLA-A and -B alleles are not clustered in the same supertype, but our PSSM/ANN predictor trained on HLA-A*0201 allele was able to make reasonable predictions for HLA-B*4002, 4501, and 5101 alleles.

Display full size

Figure 3. The cross-allele performance of the PSSM and ANN prediction models, trained on the HLA-A*0201 dataset and tested on the evaluation dataset of all 22 alleles. Aroc: area under the receiver operating characteristic curve; PSSM: position-specific scoring matrices; ANN: artificial neural networks.

Rank measure analysis of PSSM and ANN predictors with NetMHCpan 4.0 and NetCTLpan 1.1

CTL epitope processing forecast

The PSSM and ANN predictors of the EasyPred modeler described in the present study have been shown to perform best when high-sensitivity forecasts for CTL epitope identification are focused (Figure 4). Most of the MHC molecules achieve a sensitivity of 80% at the threshold determined by the score of the top 85% in the training set, while some alleles (HLA-A*0101, A*1101, and A*6801) approached the sensitivity of 100%. If focusing on optimal sensitivity, it was shown that the forecast algorithm should exclude both proteasomal cleavage and TAP forecasts, reducing the method to the MHC binding forecast alone (Figure 4). Whether this observation reflects actual biological aspects of the specificity overlap between the three pathway players or it simply occurs because the forecast of MHC class I affinity has gained accuracy during the recent years, the predictors for TAP transport efficiency and proteasomal cleavage have not much changed or been renewed [77].

Display full size

Figure 4. Sensitivity analysis of PSSM and ANN predictors, TAP-2003 [52] and TAP-2004 [53] matrices, as well as constitutive- and immuno-proteasomal cleavage matrices in identifying CTL epitopes. MHC: major histocompatibility complex; PSSM: position-specific scoring matrices; ANN: artificial neural networks; TAP: transporters associated with antigen processing; CTL: cytotoxic T lymphocytes.