INeo-Epp: A Novel T-Cell HLA Class-I Immunogenicity or Neoantigenic Epitope Prediction Method Based on Sequence-Related Amino Acid Features.

In silico T-cell epitope prediction plays an important role in immunization experimental design and vaccine preparation. Currently, most epitope prediction research focuses on peptide processing and presentation, e.g., proteasomal cleavage, transporter associated with antigen processing (TAP), and major histocompatibility complex (MHC) combination. To date, however, the mechanism for the immunogenicity of epitopes remains unclear. It is generally agreed upon that T-cell immunogenicity may be influenced by the foreignness, accessibility, molecular weight, molecular structure, molecular conformation, chemical properties, and physical properties of target peptides to different degrees. In this work, we tried to combine these factors. Firstly, we collected significant experimental HLA-I T-cell immunogenic peptide data, as well as the potential immunogenic amino acid properties. Several characteristics were extracted, including the amino acid physicochemical property of the epitope sequence, peptide entropy, eluted ligand likelihood percentile rank (EL rank(%)) score, and frequency score for an immunogenic peptide. Subsequently, a random forest classifier for T-cell immunogenic HLA-I presenting antigen epitopes and neoantigens was constructed. The classification results for the antigen epitopes outperformed the previous research (the optimal AUC = 0.81, external validation data set AUC = 0.77). As mutational epitopes generated by the coding region contain only the alterations of one or two amino acids, we assume that these characteristics might also be applied to the classification of the endogenic mutational neoepitopes also called "neoantigens." Based on mutation information and sequence-related amino acid characteristics, a prediction model of a neoantigen was established as well (the optimal AUC = 0.78). Further, an easy-to-use web-based tool "INeo-Epp" was developed for the prediction of human immunogenic antigen epitopes and neoantigen epitopes.

1. Introduction

An antigen consists of several epitopes, which can be recognized either by B- or T-cells and/or molecules of the host immune system. However, usually only a small number of amino acid residues that comprise a specific epitope are necessary to elicit an immune response [1]. The properties of these amino acid residues causing immunogenicity are unknown. HLA-I antigen peptides are processed and presented as follows: (a) cytosolic and nuclear proteins are cleaved to short peptides by intracellular proteinases; (b) some are selectively transferred to the endoplasmic reticulum (ER) by the TAP transporter, and subsequently are treated by endoplasmic reticulum aminopeptidase; and (c) antigen-presenting cells (APCs) present peptides containing 8-11 AA (amino acid) residues on HLA class I molecules to CD8+ T-cells [2]. Researchers can now simulate antigen processing and presentation by computational methods to predict binding peptide-MHC complexes (p-MHC). Several types of software systems have been developed, including NetChop [3], NetCTL [4], NetMHCpan [5], and MHCflurry [6]. However, despite that the binding to MHC molecules of most peptides is predicted, only 10%~15% of those have been shown to be immunogenic [7-10]. For neoantigens, the result was approximately 5% (range: 1%-20%) due to central immunotolerance [11, 12]. As a result, the cycle for vaccine development and immunization research is extended. Here, we aim to develop a T-cell HLA class-I immunogenicity prediction method to further identify real epitopes/neoepitopes from p-MHC to shorten this cycle.

Many experimental human epitopes have been collected and summarized in the immune epitope database (IEDB) [13], which makes it feasible to mathematically predict human epitopes. However, there still exist two limitations: (i) a high level of MHC polymorphism produces a severe challenge for T-cell epitope prediction and (ii) there is an extremely unequal distribution of data to compare epitopes and nonepitopes. It is not conducive to analyze the potential deviation existing in TCR recognition owing to the presentation of different HLA peptides. A general analysis of all HLA-presented peptides, ignoring the specific pattern of TCR recognition of individual HLA-presented peptides, may result in a lower predictive accuracy.

With the advances in HLA research, Sette and Sidney [14] classified, for the first time, overlapping peptide binding repertoires into nine major functional HLA supertypes (A1, A2, A3, A24, B7, B27, B44, B58, and B62). In 2008, Sidney et al. [15] made a further refinement, in which over 80% of the 945 different HLA-A and B alleles can be assigned to the original nine supertypes. It has not been reported whether peptides presented by different HLA alleles influence TCR recognition. Hence, we collected experimental epitopes according to HLA alleles and assumed that epitopes belonging to the same HLA supertypes have similar properties.

Moreover, screening for endogenic mutational neoepitopes is one of the core steps in tumor immunotherapy. In 2017, Ott et al. [16] and Sahin et al. [17] confirmed that peptides and RNA vaccines made up of neoantigens in melanoma can stimulate and proliferate CD8+ and CD4+ T-cells. In addition, a recent research suggests that including neoantigen vaccination not only can expand the existing specific T-cells but also can induce a wide range of novel T-cell specificity in cancer patients and enhance tumor suppression [18]. Meanwhile, a tumor can be better controlled by the combination therapy of neoantigen vaccine and programmed cell death protein 1 (PD-1)/PD1 ligand 1 (PDL-1) therapy [19, 20]. Nevertheless, a considerable number of predicted candidate p-MHC from somatic cell mutations may be false positive, which would fail to stimulate TCR recognition and immune response. This is undoubtedly a challenge for designing vaccines against neoantigens.

In our study, based on HLA-I T-cell peptides collected from experimentally validated antigen epitopes and neoantigen epitopes, we aim to build a novel method to further reduce the range of immunogenic epitope screening based on predicted p-MHC. Finally, a simple web-based tool, INeo-Epp (immunogenic epitope/neoepitope prediction), was developed for prediction of human antigen and neoantigen epitopes.

2. Materials and Methods

The flow chart for “INeo-Epp” prediction is shown in Figure 1.

Figure 1

The flow chart for “INeo-Epp” prediction.

2.1. Construction of Immunogenic and Nonimmunogenic Epitopes

Peptides that can promote cytokine proliferation are considered to be immunogenic epitopes. However, nonimmunogenic epitopes may result from the following reasons: (a) p-MHC is truly unrecognized by TCR, (b) peptides are not presented by MHC (quantitatively expressed as rank(%) > 2, see Rank(%) Score (C24) for details), and (c) negative selection/clonal presentation is induced by excessive similarity to autologous peptides [21]. In this work, to further study the recognition preferences of T-cells, peptides with >2 rank(%) were regarded as not in contact with TCR, and sequences 100% matching the human reference peptides (ftp://ftp.ensembl.org/pub/release-97/fasta/homo_sapiens/pep/) were regarded as exhibiting immune tolerance. Hence, we removed these from the definition of nonimmunogenic peptides.

2.2. Construction of Data Sets: Epitopes, External Validation of Epitopes, and Neoepitopes

Antigen epitope data were collected from IEDB (linear epitope, human, T-cell assays, MHC class I, any disease was chosen). Data collection criteria accommodated for each HLA allele quantity > 50 and frequency > 0.5% (refer to allele frequency database [22]) (Table 1, check Table S1 for detailed information).

Table 1

Summary of IEDB epitope data.

HLA supertype	IEDB HLA data	Number		HLA allele frequency Asian/Black/Caucasian	Motif view
		Negative	Positive
A1	A01:01	811	103	0.154/0.046/0.164	1-2(ST)-3-4-5-6-7-8-9(Y)
	A26:01	83	19	0.041/0.014/0.030	1(DE)-2(ITV)-3-4-5-6-7-8-9(FMY)
A2	A02:01	1883	1580	0.049/0.123/0.275	1-2(LM)-3-4-5-6-7-8-9(ILV)-10(V)
A3	A11:01	196	174	0.139/0.014/0.060	1-2(IMSTV)-3-4-5-6-7-8-9(K)-10(K)
	A03:01	1400	169	0.063/0.083/0.139	1-2(ILMTV)-3-4-5-6-7-8-9(K)-10(K)
A24	A24:02	207	219	0.136/0.024/0.084	1-2(WY)-3-4-5-6-7-8-9(FIW)
	A23:01	1138	12	0.006/0.109/0.019	1-2(WY)-3-4-5-6-7-8-9-10(F)
B7	B35:01	63	248	0.062/0.068/0.055	1-2(P)-3-4-5-6-7-8-9(FMY)
	B07:02	523	244	0.034/0.005/0.0143	1-2(p)-3-4-5-6-7-8-9(FLM)
	B51:01	13	51	0.074/0.021/0.047	1-2(P)-3-4-5-6-7-8-9(IV)
B8	B08:01	317	195	0.036/0.037/0.114	1-2-3-4-5(HKR)-6-7-8-9(FILMV)
B27	B27:05	100	86	0.008/0.008/0.037	1(RY)-2(R)-3(FMLWY)-4-5-6-7-8-9
B44	B37:01	1036	10	0.034/0.005/0.014	—
	B40:01	67	65	0.022/0.012/0.052	—
	B44:02	73	66	0.008/0.020/0.095	1-2(E)-3-4-5-6-7-8-9(FIWY)
B58	B58:01	11	62	0.041/0.037/0.007	1-2(AST)-3-4-5-6-7-8-9(W)
B62	B15:01	3	70	0.016/0.010/0.060	1-2(LMQ)-3-4-5-6-7-8-9(FY)
Total		7924	3373
Remove negative rank(%) > 2		5123	3373
Remove negative human 100% similar		4943	3373

The external antigen epitope validation set was collected from seven published independent human antigen studies [23-29], consisting of 577 nonimmunogenic epitopes and 85 immunogenic epitopes (Table 2, S2 Table).

Table 2

External data included in validation set.

Publication time	PMID	Author	Nonepitopes	Epitopes
2013	23580623	Weiskopf et al.	477	42
2018	29397015	Luxenburger et al.	100	26
2018	30260541	Xia et al.	—	1
2018	30487281	Vahed et al.	—	4
2018	30518652	Khakpoor et al.	—	2
2018	30587531	Huth et al.	—	4
2018	30815394	Sekyere et al.	—	6
Total			577	85
Remove negative with rank(%) > 2 and HLA supertypes (not appeared in training set)			321	69

Here, we removed peptides for which HLA supertypes do not appear in the training set, because we assume peptides belonging to the same HLA supertypes to have similar properties. In the external validation set, some peptides bind to rare HLA supertypes. Their characteristics were not included in the training set. Hence, these peptides in the external validation data might lead to a classification bias.

The neoantigen data were collected from 11 publications [19, 30–39] and IEDB mutational epitopes, and 13 published data sets collected by Bjerregaard et al. in one publication [40] in 2017 (see Table 3, S3 Table for details) were also included.

Table 3

Neoepitope data included in this study.

Publication time	PMID	Author	Tumor type	Nonimmunogenic neoepitopes	Immunogenic neoepitopes	T-cell assay
2013-12	24323902	D. A. Wick et al.	Ovarian cancer	—	1	ELISPOT
2015-9	26359337	E. M. Van Allen et al.	Melanoma	—	18	Clinical benefit
2015-11	26752676	T. Karasaki et al.	Lung adenocarcinoma	—	4	—
2016-1	26901407	A. Gros et al.	Melanoma	12	14	ELISPOT
2016-5	27198675	E. Strønen, et al.	Melanoma	1134	16	CTL clone
2016-12	28405493	A. Nelde et al.	Lymphoma	—	2	ELISPOT
2017-6	28619968	X. Zhang et al.	Breast cancer	—	4	Flow cytometry
2017-10	29104575	M. Markus et al.	Melanoma	10	16	—
2017-11	29187854	A.-M. Bjerregaard et al.	Polytype	1874	42	ELISPOT et al.
2017-11	29132146	V. P. Balachandran et al.	Pancreatic	—	10	Flow cytometry
2018-5	29720506	T. Matsuda et al.	Ovarian cancer	—	3	ELISPOT
2018-12	29409514	K. Sonntag et al.	Pancreatic ductal carcinoma	—	3	Flow cytometry
2018-10	30357391	V. Randi et al.	—	6	35	—
Total				3030	168
Remove duplication				2837	164
Remove negative rank(%) > 2 and human 100% similar				1697	164

2.3. Construction of Potential Immunogenicity Feature

2.3.1. Calculation of Peptide Characteristics Based on Amino Acid Sequences

The formula for calculating peptide characteristics is shown in (1). P, P2, and P (N-terminal, position 2, C-terminal as anchored sites by default) are considered to be embedded in HLA molecules and have no contact with TCRs; therefore, they were not evaluated. where P is peptide, c is characteristic. P represents the characteristics of peptides, A represents amino acids, N represents the N-terminal in a peptide, C represents the C-terminal in a peptide, Pos represents the amino acid position in a peptide, and P represents characteristics of amino acids in peptides.

2.3.2. Frequency Score for Immunogenic Peptide (C22)

Amino acid distribution frequency differences between immunogenic and nonimmunogenic peptides at TCR contact sites (excluding anchor sites) were considered as a feature: where Pie+ represents immunogenic peptides, Pie− represents nonimmunogenic peptides. f′ represents amino acid frequency in the TCR contact position. Pie+(f′) represents the frequency of amino acids in immunogenic peptides at TCR contact sites.

2.3.3. Calculating Peptide Entropy (C23)

Peptide entropy [41] was used as a feature: where P represents peptide entropy. f represents amino acid frequency in the human reference peptide sequence. P represents the frequency in the human reference peptide sequence of amino acids in epitope peptides.

2.3.4. Rank(%) Score (C24)

HLA binding prediction was performed using NetMHCpan 4.0. Rank(%) provides a robust filter for the identification of MHC-binding peptides, in which rank(%) was recommended as an evaluation standard, rank(%) < 0.5 as strong binders, 0.5 < rank(%) < 2 as weak binders, and rank(%) > 2 as no binders.

2.4. Fivefold Cross-Validation, Feature Selection, Random Forests, and ROC Generation

The 5-fold cross-validation was implemented in R using the caret package [42] (method = “repeatedcv,” number = 5, repeats = 3). The feature screening results were generated in R using the package Boruta [43] (a novel random forest-based feature selection algorithm for finding all relevant variables, which provides unbiased and stable selection of important and nonimportant attributes from an information system). It iteratively removes the features which are proven by a statistical test to be less relevant than random probes. It uses Z score (computed by dividing the average loss by its standard deviation) as the importance measure, and it takes into account the fluctuations of the mean accuracy loss among trees in the forest. R package randomForest [44] was used for training data (the R language machine learning package caret provides automatic iteration selection of optimal parameters: mtry = 15 for antigen epitope and mtry = 14 for neoantigen epitope; the remaining parameters use default values). R package ROCR [45] was used for drawing ROC.

2.5. Web Tool Implementation

The front end of Ineo-Epp was constructed via HTML/JavaScript/CSS. The back end was written in PHP, connecting the web interface and Apache web server. A python script was used for calculating peptide characteristics and extracting mutation information. Models were built using R.

3. Results

Ultimately, 11,297 validated epitopes and nonepitopes with lengths of 8-11 amino acids were collected from IEDB. T-cell responses included activation, cytotoxicity, proliferation, IFN-γ release, TNF release, granzyme B release, IL-2 release, and IL-10 release. Seventeen different HLA alleles were collected (Figure 2(a)), and the detailed antigen length distribution is shown in Figure 2(b). Additionally, we collected the neoantigen data from 12 publications, including 2837 nonneoepitopes and 164 neoepitopes (Figure 2(c)), and the detailed neoantigen length distribution is shown in Figure 2(d).

Figure 2

Epitope/neoepitope peptide composition and amino acid length distribution. (a) Detailed data distribution of seventeen HLA alleles of antigen peptides, the proportion of each HLA allele (positive and negative) epitopes, and the corresponding HLA frequency in Asians, Blacks, Caucasians. (b) Proportion of antigen peptides with lengths of 8-11 AA. (c) Data distribution of HLA alleles of neoantigen peptides. (d) Proportion of neoantigen peptides with lengths of 8-11 AA.

The TCR contact position plays a crucial role in the analysis of immunogenicity, as TCRs might be more sensitive to some amino acids; the amino acid preference in the antigen epitope peptide and the antigen nonepitope peptide was further analyzed after excluding anchor sites (N-terminal, position 2, and C-terminal) (Figure 3). We found that TCRs tend to identify hydrophobic amino acids. For example, 3/4 hydrophobic amino acids (L, W, P, A, V, and M) occur more frequently in immunogenicity epitopes. Charged amino acids (e.g., D and K) are enriched in nonepitopes, whereas the rest of the charged amino acids (R, H, and E) show no difference. Based on the result in Figure 3, the amino acid distribution difference at the TCR contact sites was regarded by us as one of the immunogenicity features (i.e., Frequency Score for Immunogenic Peptide (C22)).

Figure 3

Antigen epitope amino acid distribution frequency in the TCR contact site of epitopes and nonepitopes. Frequency distribution of amino acids at TCR contact sites in antigen epitope and nonepitope peptides, and the amino acids below the dotted line are preferred by the epitope.

3.1. Classification Prediction Model for Antigen Epitopes

We constructed the features of peptides on the basis of the characteristics of amino acids (see Calculation of Peptide Characterstics Based on Amino Acid Sequences). All amino acid characteristics were selected from ProtScale [46] in ExPASy (SIB Bioinformatics Resource Portal). The 21 involved features are as follows: Kyte-Doolittle numeric hydrophobicity scale (C1) [47], molecular weight (C2), bulkiness (C3) [48], polarity (C4) [49], recognition factors (C5) [50], hydrophobicity (C6) [51], retention coefficient in HPLC (C7) [52], ratio hetero end/side (C8) [49], average flexibility (C9) [53], beta-sheet (C10) [54], alpha-helix (C11) [55], beta-turn (C12) [55], relative mutability (C13) [56], number of codon(s) (C14), refractivity (C15) [57], transmembrane tendency (C16) [58], accessible residues (%) (C17) [59], average area buried (C18) [60], conformational parameter for coil (C19) [55], total beta-strand (C20) [60], and parallel beta-strand (C21) [61] (see Table S4 for details). Also, Frequency Score for Immunogenic Peptide (C22), Calculating Peptide Entropy (C23), and Rank(%) Score (C24) were also taken into consideration. Together, 24 immunogenic features were collected, and all features were retained for antigen epitope prediction after screening using the R package Boruta. Compared with other characteristics, the frequency score for immunogenic peptide and rank(%) have higher impact, suggesting that they have more significant influence on antigen epitope classification (Figure 4(a)).

Figure 4

Feature selection in antigen epitopes and ROC curves of antigen epitope classification. (a) Peptide features: twenty-four features were screened, and we defined the features on the right of the dotted line as being effective. (b) Trained model: the line in blue represents antigen epitopes without screening; the line in green represents the selection with the deletion of the rank(%) > 2 nonepitope; the line in red represents the selection with the deletion of the nonepitopes 100% matching the human reference peptide sequence. (c) External validation: the ROC curves for the external verification set. The line in purple represents modeling using antigen epitopes without filtering, and the line in pink represents modeling using antigen epitopes removing nonepitopes with rank(%) > 2 and HLA for which supertypes did not appear in the training set.

The receiver operator characteristic (ROC) curve of models are shown in Figure 4. The fivefold cross-validation AUC was 0.81 in the prediction model for the antigen epitope (line in red, Figure 4(b)), and the externally validated (see Table 2) AUC was 0.75 (line in purple, Figure 4(c)). Here, we tried to remove peptides for which HLA supertypes did not appear in the training set from the externally validated antigen data, and the AUC, specificity, and sensitivity were increased to 0.78, 0.71, and 0.72, respectively (line in pink, Figure 4(c)). This, to some extent, verifies our conjecture about TCR specific recognition of different HLA alleles presenting peptides.

3.2. Classification Prediction Model for Neoantigen Epitopes

Neoantigens derived from somatic mutations are different from the wild peptide sequences. Therefore, some mutation-related characteristics were also taken into account. For instance, difference in hydrophobility before and after mutation (C25), differential agretopicity index (DAI, C26) [62], and whether the mutation position was anchored (C27). Finally, 27 features were selected for the neoantigen epitope prediction model. However, only 25 neoantigen-related features were retained after running Boruta, because C25 and C27 were removed. Also, rank(%) showed a marked effect (Figure 5(a)). In the fivefold cross-validation of the prediction model for neoantigen epitopes, AUC was 0.78 (Figure 5(b)).

Figure 5

Feature selection in neoantigen epitopes and ROC curves of neoantigen epitope classification. (a) Twenty-seven features were screened, and the 25 features on the right of the dotted line were reserved for modeling using a random forest algorithm. (b) ROC curves of neoantigen epitope classification.

3.3. Web Server for TCR Epitope Prediction

Based on the abovementioned validated features, we established a web server for TCR epitope prediction, named “INeo-Epp.” This tool can be used to predict both immunogenic antigen and neoantigen epitopes. For antigens, the nine main HLA supertypes can be used. We recommend the peptides with the lengths of 8-12 residues, but not less than 8. N-terminal, position 2, and C-terminal were treated as anchored sites by default. A predictive score value greater than 0.5 is considered as immunogenicity (positive-high), a score between 0.4 and 0.5 is considered as positive-low, and a score less than 0.4 is considered as negative-high. It is critical to make sure that the HLA-subtype must match your peptides (rank(%) < 2). Where HLA-subtypes mismatch, a large deviation of the rank(%) value may strongly influence the results. Additionally, the neoantigen model requires providing wild type and mutated sequences at the same time to extract mutation-associated characteristics, and currently only immunogenicity prediction for neoantigens of single amino acid mutations are supported. Users can choose example options to test the INeo-Epp (http://www.biostatistics.online/ineo-epp/neoantigen.php).

4. Discussion

Due to the complexity of antigen presenting and TCR binding, the mechanism of TCR recognition has not been clearly revealed. In 2013, Calis et al. [63] developed a tool for epitope identification for mice and humans (AUC = 0.68). Although mice and human beings are highly homologous, the murine epitopes may very likely cause limitations in identifying human epitopes. Inspired by J. A. Calis, our research here focused on human beings' epitopes and has been conducted in a larger data set.

By analyzing epitope immunogenicity from the perspective of amino acid molecular composition, we observed that TCRs do have a preference for hydrophobic amino acid recognition. For short peptides presented by different HLA supertypes, TCRs may have different identification patterns. The immunogenicity prediction based on all HLA-presenting peptides may affect the accuracy of the prediction results. That is, if the prediction could focus on specified HLA-presenting peptides, the results may improve. Therefore, in our work we used HLA supertypes to improve the prediction of HLA-presenting epitopes, including antigen epitopes and neoantigen epitopes, for a better recognition by TCRs. At present, neoantigen epitopes that can be collected in accordance with the standard for experimental verification are too few, the data of positive and negative neoantigens are unbalanced, and there is not enough data to be used for an external verification set. In the future, we will continue to refine and expand our training and verification datasets. Recently, Laumont et al. [64] demonstrated that noncoding regions aberrantly expressing tumor-specific antigens (aeTSAs) may represent ideal targets for cancer immunotherapy. These epitopes can also be studied in the future. Increased epitope data may also help empower the prediction of potentially immunogenic peptides or neopeptides.

5. Conclusions

Neoantigen prediction is the most important step at the start of preparation of a neoantigen vaccine. Bioinformatics methods can be used to extract tumor mutant peptides and predict neoantigens. Most current strategies aimed at and ended in presenting peptide predictions, and among the results of these predictions, probably only fewer than 10 neoantigens might be clinically immunogenic and produce effective immune response. It is time-consuming and costly to experimentally eliminate the false positively predicted peptides. Our methods as developed in this study and the INeo-Epp tool may help eliminate false positive antigen/neoantigen peptides and greatly reduce the amount of candidates to be verified by experiments. We believe that in the age of biological system data explosion, computational approaches are a good way to enhance research efficiency and direct biological experiments. With the development of machine learning and deep learning, we expect that the prediction of epitope immunogenicity will be continually improved.

In summary, this study provides a novel T-cell HLA class-I immunogenicity prediction method from epitopes to neoantigens, and the INeo-Epp can be applied not only to identify putative antigens, but also to identify putative neoantigens.

It needs to be stated here that we published the preprint [65] of this article in July 2019. This is a modified version.