Machine-Learning Prediction of Tumor Antigen Immunogenicity in the Selection of Therapeutic Epitopes

T cells can affect antitumor immune responses through recognition of tumor-specific antigens (TSA) presented by major histocompatibility complex (MHC) proteins. These peptides include tumor neoantigens, which are classically thought of as derived from mutation-containing proteins that generate novel immunogenic epitopes (1). Despite the ability of neoantigen therapeutic vaccines to promote tumor-specific T-cell responses in a number of preclinical models (2–4), clinical efficacy has yet to be demonstrated (5, 6). A significant challenge for translation of TSA therapies is the ability to select the subset of clinically relevant epitopes from all computationally predicted neoantigens. Many neoantigen prediction algorithms rely heavily on peptide/MHC binding affinity predictions to rank epitopes (7–14). Unlike murine preclinical models, where in vivo/ex vivo methods to further screen for immunogenicity can be applied (3, 15), no such benchtop prediction method for immunogenicity is currently available for humans. We have previously demonstrated in multiple murine models that the number of predicted neoantigens is much higher than the number of confirmed immunogenic neoantigens (15). Studies demonstrate that in some tumors, the number of predicted neoantigens is far greater than the number of immunogenic neoantigens that have been identified in mouse models (16, 17). As such, the development of an algorithm to predict the immunogenicity of neoantigen peptides (i.e., variant peptides predicted to bind MHC) would be valuable for screening predicted neoantigens for clinical application.

In addition to conventional single-nucleotide variant (SNV) neoantigens, studies have suggested the presence of tumor-specific mRNA splice variants (18, 19), expression of noncoding regions (20), and alternative ribosomal products (21–28), allowing an out-of-frame translation to occur outside the setting of an insertion/deletion (INDEL) mutation. An increasing need exists to define frequencies of predicted TSAs existing in an out-of-frame context, their clinical implications, and whether frame-filtering should be applied for computational neoantigen prediction. In the context of SNV tumor antigen prediction, allowing for out-of-frame calls may identify “pseudo-SNV” antigens (i.e., out-of-frame antigens that contain concurrent SNV mutations) with immunogenicity responses similar to what is observed in frameshift neoantigens. As a preliminary approach to identify both in- and out-of-frame neoantigens, we performed SNV tumor-antigen computational screening across all open reading frames, looking for (i) the correlates of immunogenicity for these predicted neoantigens and (ii) the capacity for out-of-frame epitopes to drive antitumor immunity.

Features associated with neoantigen immunogenicity remain unclear. Here, we have elucidated peptide-intrinsic features significantly associated with vaccine/IFNγ ELISpot–derived immunogenicity scores of MHC class I and II TSAs. Using gradient boosting with cross-validation, we developed an algorithm to predict MHC I TSA peptide immunogenicity based on peptide-intrinsic biochemical features. We modeled the immunogenicity of predicted neoantigens in the BBN963 basal-like bladder cancer model and demonstrated the capacity of epitopes with high predicted immunogenicity to control tumor growth significantly better than those with low predicted immunogenicity. This algorithm was additionally validated using graft-versus-leukemia (GvL) minor histocompatibility mismatch antigens (mHA) in the P815 mastocytoma allogeneic transplant model. Applying this algorithm to predicted MHC I neoantigens from The Cancer Genome Atlas (TCGA) pan-cancer data set, we observed significant positive association between highly immunogenic neoantigens (HIN; in the top 95th percentile of predicted immunogenicity score) and microsatellite instability (MSI) high–driven immune features in colon adenocarcinoma (COAD) and significant negative association between signatures of anti–PD-1 therapy responsiveness and HIN numbers in lung adenocarcinoma (LUAD) cancer types. Lastly, we provide evidence in favor of antitumor cytotoxic T-cell responses generated against a predicted out-of-frame neoantigen, suggesting a proportion of predicted out-of-frame SNV tumor antigens may be presented by the tumor to generate an immune response. Prediction of peptide immunogenicity on a framework of peptide/MHC binding should improve understanding of antitumor T-cell responses and neoantigen selection for therapeutic vaccine applications.

The B16F10 cell line was purchased from ATCC (CLR-6475) and cultured according to the ATCC protocol. The P815 cell line was purchased from ATCC (TIB-64), transduced with luciferase as previously described (29), and cultured according to the ATCC protocol. The BBN963, UPPL1541, and MB49 cell lines were obtained and passaged as previously described (15). The T11 model was obtained and passaged as previously described (30). All cells used in this study were derived from viably frozen stocks of the above cell lines, with aliquots derived within ≤5 passages of the original stock. No Mycoplasma testing was performed. No further authentication was performed on cell lines directly purchased from ATCC (B16F10, P815) or those received directly from the deriving lab (T11: Charles M. Perou, UNC Lineberger; BBN963, UPPL1541: William Y. Kim, UNC Lineberger). The MB49 cell line was authenticated through transcriptomic analysis, as previously described (15).

All experiments described in this study were approved by the UNC Institutional Animal Care and Use Committee (IACUC). Animals used in this study, their vendor source, and respective tumor cell lines included C57BL/6J (Jackson Laboratories; B16F10), C57BL/6 (Charles River Laboratories; BBN963, MB49, UPPL1541), DBA/2J (Jackson Laboratories; P815), and BALC/c (Jackson Laboratories; T11). Tumor injection routes and cell numbers for all models and experiments included B16F10: flank subcutaneous (s.c.), 10⁵ cells; BBN963: flank s.c., 10⁷ cells; MB49: flank s.c., 10⁵, UPPL1541: Flank s.c., 10⁶, T11: mammary fat pad intradermal, 10⁴, P815: tail-vein intravenous, 3 × 10⁵. Tissue collection and DNA/RNA isolation is described in the “Neoantigen and mHA prediction” section below. Graft-versus-host disease (GvHD) scoring was performed as previously described (31), with score defined as the sum of five components of posture, fur, activity, skin, and weight loss on a 0 to 2 scale.

All single-cell suspensions mentioned in the below methods sections were derived using the below listed protocol. Tissues were homogenized in cold PBS using the gentleMACS Dissociator and the samples were passed through a 70-μm cell strainer using a 5-mL syringe plunger. The samples were centrifuged for seven minutes at 290 RCF, 4°C, decanting the supernatant. The remaining pellet was resuspended into 1 mL of ACK lysis buffer (150 mmol/L NH4Cl, 10 mmol/L, KHCO3, 0.1 nmol/L Na2EDTA in DPBS, pH 7.3) for 2 minutes at room temperature before quenching with 10 mL of cold media. The samples were centrifuged for 7 minutes at 290 RCF, 4°C, resuspended in 10 mL of cold media, and passed through a 40-μm cell strainer.

Neoantigen prediction was performed as previously described (15). Briefly, mice were inoculated with tumor cells (Fig. 1A) in the route and counts listed above, and monitored until tumor size reached 100 mm³ by caliper measurement (⁠|$\frac{{l\ \times \ {w^2}}}{2}$|⁠, where w is the smaller of two perpendicular tumor axes), at which point mice were humanely sacrificed with CO₂ asphyxiation followed by cervical dislocation. RNA was extracted from single-cell suspensions of tumors using Qiagen RNeasy Mini kit (cat. #74104), and DNA was extracted from single-cell suspensions of tumors and matched-normal tail clippings or livers using Qiagen DNeasy kit (cat. #69504), all according to the manufacturer's protocol. Whole-exome and transcriptome library preparation was performed using Agilent SureSelect XT All Exon and Illumina TruSeq Stranded mRNA library preparation kits, respectively. In addition to tumor-derived samples, P815 tumor WES samples were additionally derived from cell line culture (10⁵ cells per sample). Libraries were sequenced via 2 × 100 runs on an Illumina HiSeq 2500 at the UNC High-Throughput Sequencing Facility (HTSF). Tumor mutations were called using UNCeqR (https://lbc.unc.edu/∼mwilkers/unceqr_dist/; ref. 32), filtering for SNV mutations with at least 5× coverage by RNA sequencing (RNA-seq). Translated 8–11 mer (class I) or 15 mer (class II) peptides were derived across all open reading frames, and then predicted for MHC binding affinity using NetMHCPan2.8 (MHC I; http://www.cbs.dtu.dk/services/NetMHCpan-2.8/) and NetMHCIIPan3.0 (MHC II; http://www.cbs.dtu.dk/services/NetMHCIIpan-3.0/). Class I minor mismatch antigens were predicted similarly in the P815 model against the BALB/c histocompatible donor. Predicted binders were filtered by binding affinity <500 nmol/L, generally accepted cutoff for immunogenicity as previously noted (6, 9, 33), with top epitopes screened for immunogenicity using a vaccine/ELISpot approach, as described below. Matched-normal samples for the T11 tumor model (BALB/c) have been previously deposited (https://www.ncbi.nlm.nih.gov/sra/?term=SRR3502701). All additional sequencing data can be accessed at the Gene Expression Omnibus (GEO; accession: GSE136619).

Predicted neoantigen peptides (MHC I: n = 210; MHC II: n = 68) were synthesized by New England Peptide, using custom peptide array technology (Supplementary Data File 1). Non–tumor-bearing wild-type animals were vaccinated with predicted neoantigen peptides of their respective predicted haplotype, given as a subcutaneous injection of a pool of 8 equimolar peptides (5 nmol total peptide) and 50 μg poly(I:C; Sigma; cat. #P1530) in PBS. A second identical injection was repeated 6 to 7 days after primary injection. Mice were humanely sacrificed with CO₂ asphyxiation followed by cervical dislocation 5 to 6 days after the second injection. Spleens were harvested and prepared into single-cell suspension, as described above. Splenocytes were plated in triplicate at 5 × 10⁵ cells per 100 μL media (RPMI-1640; Gibco cat. #11875-093) with 10% FBS (Gemini; cat. #900-108) onto an IFNγ capture antibody-coated ELISPOT plate (BD Biosciences; cat. #551083) according to the manufacturer protocol for 48 to 72 hours, along with 1 nmol of a single peptide against which the respective mouse was vaccinated. Immunogenicity was defined as the average number of spot-forming cells (SFC) identified using an ELISpot plate reader (AID Classic ERL07), with no-peptide background subtracted from each epitope.

Variables used in neoantigen immunogenicity regression and modeling were derived using the “aaComp” command of the R package “Peptides” (v2.4; https://cran.r-project.org/web/packages/Peptides/index.html). Using features derived from this command (Tiny, Small, Aliphatic, Aromatic, Non-polar, Polar, Charged, Basic, Acidic), variables were derived by the presence (1) or absence (0) of each feature at each absolute and relative position along each antigen, at the site of SNV mutation along each antigen, at the first or last 3 amino acid residues (beginning/end) or middle residues (middle) of each antigen, or difference (loss: −1, gain: 1, or no change: 0) of each feature in the mutated versus reference antigen along SNV mutation site.

For all generalized linear models (GLM) and predictive modeling analyses, low variance variables (defined by the “nearZeroVar” function of the “caret” package) were removed prior to further analysis. GLM using the R “glm” function were used for all nonmodeling univariable and multivariable linear regression analyses, with significance reported as false discovery rate (FDR)-adjusted P values (q-value) using the R “p.adjust” command. Backward stepwise regression for multivariable modeling was performed using the R “stepAIC” command of the “MASS” package (https://cran.r-project.org/web/packages/MASS/index.html), optimized on Akaike Information Criterion (AIC). Backward stepwise regression was performed by starting with all variable candidates and testing the deterioration of the model with removal of each variable.

For immunogenicity prediction modeling, analyses were performed using the R package “caret” as a wrapper for running each multivariable approach: GLM, elastic net, random forest, gradient boosting, and linear and radial support vector machine methods. For cross-validation, data were split into exploration (n = 141) and validation (n = 69) sets using the caret “createDataPartition” function, confirming statistically nonsignificant differences in measured immunogenicity between exploration and validation sets (Mann–Whitney P > 0.2). Model performance was derived from Pearson correlation coefficients between ELISpot immunogenicity and predicted immunogenicity scores, using a 10,000-fold cross-validation (2/3rd random resampling) approach within the exploration set, with the input predictor variables limited to those that demonstrated significant univariable correlation in >50% of 1,000-fold bootstrapping iterations (2/3rd resampling) within the exploration set. The final gradient boosting machine-learning algorithm immunogenicity prediction of MHC I epitopes can be accessed at https://github.com/vincentlaboratories/neoag.

To explore for computation evidence of out-of-frame transcripts, StringTie (https://ccb.jhu.edu/software/stringtie/; ref. 34) and Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki; ref. 35) were used for de novo assembly of transcripts from BBN963 RNA-seq data, according to standard workflow provided in the above links.

BBN963 basal-like bladder cancer model treatment began with pretumor vaccination with 30 μg of a single peptide (or no-peptide control) and 50 μg poly(I:C) adjuvant injected in 100 μL PBS intradermally in the flank of 8- to 10-week-old female C57BL/6 mice (Charles River). Twelve days after vaccination, 1 × 10⁷ BBN963 cells were injected in 100 μL PBS subcutaneously in the flank, ipsilateral to the vaccine site. On day 21 after primary vaccination, animals were given a vaccine booster with 30 μg of the initial respective peptide with no poly(I:C) adjuvant. This booster was delivered in 100 μL PBS intradermally in the skin directly adjacent to the tumor. Animals were monitored for tumor growth via caliper measurement and survival every 2 to 3 days for the remainder of the study, with UNC IACUC defined endpoints of area >200 mm² or ulceration >5 mm in the longest diameter.

For P815 treatment studies, 8- to 12-week-old male BALB/c donors (Jackson Laboratory) were vaccinated on days 0 and 7 with 100 μg total peptide (3–4 pooled equimolar peptides, or no-peptide control) and 50 μg poly(I:C) adjuvant in 100 μL PBS intradermally in the flank. DBA/2 recipients were treated with 800 rad total body irradiation on day 13. On day 14, splenic-derived T cells and bone marrow cells were isolated from donor BALB/c animals. T cells were isolated from single-cell splenocyte suspensions uisng the Miltenyi Pan T-Cell Isolation Kit II (cat. #130-095-130), according to the manufacturer's protocol. Bone marrow cells were isolated as previously described (36). Recipient DBA/2 animals were given tail-vein i.v. injections of 3 × 10⁶ T cells, 3 × 10⁶ bone marrow cells, and 3 × 10⁵ P815-luciferase tumor cells (or bone marrow–only control). DBA/2 recipients were given a booster vaccine on day 21 after primary vaccine (100 μg total peptide, 50 μg poly(I:C)), with animals monitored every 2 to 3 days for survival, with UNC IACUC defined endpoints of bilateral hind-limb paralysis. Luciferase imaging studies were performed on days 8, 13, 22, 26, and 35 after transplant, using an IVIS imaging system on animals given 3 mg intraperitoneal D-luciferin (PerkinElmer; cat. #122799) 10 minutes prior to imaging.

Peptide/MHC tetramer and cell-surface protein staining were performed as described previously (37). Briefly, viable, single-cell suspensions derived from tumors (approximately 10⁷ total cells) were treated with 50 nmol/L dasatinib (Sigma-Aldrich; cat. #CDS023389) for 30 minutes at 37°C, and then stained using approximately 10 μg/mL tetramer on ice for 30 minutes. Tetramers were generated using the MBL Quickswitch Quant H-2 Kb Tetramer Kit-PE (cat. #TB-7400-K1), using peptides VALLPSVMNL or SIINFEKL irrelevant control, according to the manufacturer's protocol. Cells were then washed and incubated on ice with biotin-conjugated anti-PE antibody (5 μg/mL; BioLegend; PE001) for 20 minutes, followed by 2 washes, then further incubation with streptavidin, R-PE conjugate (SAPE, 5 μg/mL) for 10 minutes on ice. Cells were then washed and stained for viability using BD fixable viability dye FVS620 according to the manufacturer's directions. Last, cells were Fc blocked (anti-mouse CD16/CD32; 2.4G2, BD Biosciences) for 10 minutes on ice, followed by surface staining for 20 minutes on ice with the following markers: CD45 (BV510; 30-F11), CD4 (FITC; RM4-5), CD8 (APC/Fire-750; 53-6.7; all antibodies purchased from BD Biosciences). Aquisition was performed using a BD LSRFortessa flow cytometer. FlowJo flow cytometry software version 10 was used for analyses of all flow-cytometric data. Cells were selected using gates defined by single color controls and FMO or irrelevant tetramer controls, with tetramer-positive/negative CD8⁺ T cells defined within live, singlet (by FSC-A versus FSC-H), CD45⁺, CD4⁻/CD8⁺ gates.

For tetramer-sorted T-cell isolation, CD8⁺ T cells were isolated from BBN963 tumor single-cell suspensions (as described above) using the Miltenyi Dead Cell Removal Kit (130-090-101) followed by the Miltenyi CD8a⁺ T-Cell Isolation Kit (130-104-075), both according to the manufacturer's protocol. CD8⁺ T cells stained with tetramer as described above and sorted on the BD FACSJazz. Tetramer-sorted T cells or column-sorted (Miltenyi CD8a⁺ T-Cell Isolation Kit) CD8⁺ T cells from OT1 (C57BL/6-Tg(TcraTcrb)1100Mjb/J, Jackson Laboratories; cat. #003831) splenocytes were cultured in complete RPMI media [RPMI-1640 (Gibco; cat. #11875-093) with 10% FBS (Gemini; cat. #900-108), 1% sodium pyruvate (100 nmol/L; Gibco cat. #11360-070), 1% nonessential amino acids (10 mmol/L; Gibco cat. #11140-050), 1% l-glutamine (Gemini cat. #400-106), 1% HEPES buffer (1 mL; Corning; cat. #25-060-Cl), and 0.1% 2-mercaptoethanol (55 nmol/L; Gibco cat. #2198502)] in the presence of recombinant murine IL7 (10 ng/mL, PeproTech cat. #217-17), IL15 (10 ng/mL, PeproTech cat. #210-15), and IL2 (100 IU/mL, PeproTech cat. #212-12) for 72 hours. Antigen-specific T-cell expansion was performed using a previously described protocol (38). Briefly, all sorted T cells were recovered at 10⁶ cells/mL for 72 hours in RPMI complete media in the presence of IL7 (10 ng/mL), IL15 (10 ng/mL), and IL2 (10 ng/mL) at 37°C and 5% CO₂. T cells were then cocultured in RPMI complete media and IL7, IL15, and IL2, alongside peptide-pulsed DCs (2.5 μg/mL peptide) which had been pulsed approximately 18 hours prior to coculture. Media and cytokines were changed every 2 to 3 days, letting cells expand for 7 to 10 days after coculture with peptide-pulsed dendritic cells before downstream assays.

For flow-cytometric–based cytotoxicity assays, target cells (BBN963 or an irrelevant splenocyte control) were prelabeled in 5 μmol/L CFSE for 15 minutes prior to coculture. Tetramer-sorted and antigen-expanded T cells (per above section) were cocultured alongside targets at a 1:1 ratio, with 1 × 10⁵ of each target and effector population. Cells were plated on a v-bottom 96-well polypropylene plate, centrifuged at 300 × g for 1 minute, and incubated at 37°C, 5% CO₂ for 4 hours. After incubation, cells were stained using FVS700 viability dye (BD Biosciences; cat. #564997) according to the manufacturer's directions. Aquisition was performed on a BD LSRFortessa flow cytometer. FlowJo flow cytometry software version 10 was used for analyses of all flow-cytometric data. Cells were identified as targets (CFSE⁺) or effectors (CFSE⁻), looking for percent viability among targets. Percent killing was reported as frequency of dead targets, background subtracted from no-effector control wells.

The cytotoxic activity of T cells was evaluated using a standard 4-hour ⁵¹Cr release assay (39). In brief, 5 × 10^{3 51}Cr-labeled (PerkinElmer cat. #NEZ030001MC) BBN963 target cells per well were plated in triplicate in a 96 well v-bottom plate with different ratios (10:1 and 5:1 effector:target) of effector cells and incubated for 4 hours at 37°C. The supernatant was collected and analyzed with a gamma-counter (PerkinElmer). Before labeling, target cells were incubated for 2 hours at 37°C with the specific peptides (100 nmol/L) and washed twice with complete medium. Target cells were incubated with medium alone or in 1% Triton X-100 (Sigma-Aldrich) to determine the spontaneous and maximum ⁵¹Cr release, respectively. The mean percentage of specific lysis of triplicate wells was calculated as follows: [(test counts − spontaneous counts)/(maximum counts − spontaneous counts)] × 100%.

MapSplice-aligned, RSEM-quantified RNA-seq expression matrices and survival data were downloaded from FireBrowse (http://firebrowse.org/). Expression matrices were merged between all cancer types, upper quartile normalized within each sample, and log₂ transformed. Immune gene signatures (IGS) were derived from previously described signatures (40–44), with expression calculated as the mean expression of each gene within the signature. TCGA LAML samples were omitted from analysis in order to prevent skewing of IGS patterns. Inclusion criteria were those defined by TCGA pan-immune working group, according to previous studies (17). MHC I neoantigen expression used for machine-learning algorithm immunogenicity prediction was obtained from publicly available data derived in previous studies (17). TCGA pan-cancer data set (n = 11,092; LUAD n = 515; COAD n = 283) analyses were performed according to the above “Computational analysis” methods section.

Differential gene-expression analysis was performed using DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html; ref. 45). Gene set enrichment analysis (http://software.broadinstitute.org/gsea/index.jsp; ref. 46), Ingenuity pathway analysis (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/; ref. 47), and DAVID gene ontology analysis (https://david.ncifcrf.gov/; ref. 48) were performed from respective web portals. T-cell receptor diversity analysis was performed from previously published MiXCR-derived TCR reads (17)—MiXCR is an analytic tool for TCR inference from whole-transcriptome RNA-seq data (49).

Statistical analyses for survival (displayed as Kaplan–Meier plots) were performed using log-rank test, with no statistical correction. Differences in cytotoxicity in tetramer-sorted cytotoxicity assays were determined via Welch t test, with no statistical correction. Differences in tetramer-positive T-cell populations were determined via Mann–Whitney U test, with no statistical correction. All above analyses were performed in GraphPad Prism 8. All other analyses and corresponding statistical tests are described in the “Computational analysis” section.

Neoantigens and mHA were predicted in six tumor models (B16F10, BBN963, MB49, UPPL1541, P815, and T11), allowing us to characterize neoantigens in the H2^b and H2^d haplotypes (Fig. 1B). We predicted a total of 210 MHC I epitopes and 68 MHC II epitopes and determined their immunogenicity using a vaccine/ELISpot screening approach. Distribution of epitope ELISpot scores (SFC) varied by model, with MB49, B16F10, and P815 tumors including nine of the 10 most immunogenic epitopes, and BBN963 including seven of the 10 least immunogenic epitopes (Supplementary Fig. S1). MHC II epitopes were not predicted for P815 GvL mHA. With the goal of identifying peptide-intrinsic features that associated with and predicted for immunogenicity, we derived a set of features for each peptide, including the amino acid sequence and characteristic at each (I) absolute position, (II) relative site, (III) site of mutation and changes in amino acid sequence and characteristic at mutational site, and (IV) presence of amino acid or characteristic at the beginning, middle, or end of each peptide (Fig. 1C).

Univariable regression considering intrinsic peptide features as the predictor variable and immunogenicity (measured as IFNγ ELISpot values of T cells derived from vaccinated mice) of class I antigens demonstrated 38 significant features (q-value < 0.05; Fig. 2A) associated with ELISpot immunogenicity. Among these features, the most significant positive associations were changes in the mutation position to a small amino acid (Mutated_position_change_of_Small_feature), valine at relative site 2 (“Relative_site_2_V”), and basic amino acids of the reference sequence at the mutated position (Reference_AA_at_mutated_position_Basic). In contrast, the most signficant negative associations were small amino acids of the reference sequence at the mutated position (Reference_AA_at_mutated_position_Small), changes in the mutation position to a basic amino acid (Mutated_position_change_of_Basic_feature), and polar amino acids at position 6 (Absolute_position_6_Polar). We additionally sought to determine the correlation among the 38 significant features, observing relatively low correlation (Fig. 2B). Features that demonstrated significant correlation were related, such as (i) the amino acid at the mutated position of the reference sequence with charged or basic features, (ii) valine or small amino acids at absolute position 11, and (iii) the presence of a valine or small amino acid at the last position and valine at relative site 8.

Next, we evaluated the independence of the variables identified using univariable analysis using multivariable regression. To increase confidence of our multivariable regression, we performed backward stepwise regression, optimizing on AIC, as described in the Materials and Methods. Variables whose loss resulted in an insignificant change to the model performance (as measured by the AIC) were removed from the set of variables until no further variables could be removed without a significant decrease in model fit. Sixteen significant features from this step were inputted into multivariable regression. The resulting model of 8 features indicated 33.4% variation (P < 0.0001) in immunogenicity was explained by the prediction (Supplementary Fig. S2A). Significant features of the multivariable model included valine at the last position (Last_position_V; P = 0.0001), tyrosine at position 3 (Absolute_position_3_Y; P = 0.003), changes in the mutated position to a small amino acid (Mutated_position_change_of_Small_feature; P = 0.007), cysteine at relative site 4 (Relative_site_4_C; P = 0.012), lysine at relative site 5 (Relative_site_5_K; P = 0.015), tiny amino acids at relative site 6 (Relative_site_6_Tiny; P = 0.016), basic amino acids of the reference sequence at the mutated position (Reference_AA_at_mutated_position_Basic; P = 0.027), and valine at relative site 2 (Relative_site_2_V; P = 0.041; Supplementary Fig. S2B; Supplementary Table S1). To ensure this model was accurately representing both H^b and H^d haplotypes, we tested the immunogenicity for each of these five significant features, split categorically, which demonstrated similar trends between both haplotypes (Supplementary Fig. S3). We did not observe significant differences in ELISpot immunogenicity among predicted in-frame (n = 131) and out-of-frame (n = 79) epitopes, emphasizing that peptide-intrinsic features were the primary driver for immunogenicity (P > 0.05; Supplementary Fig. S4).

Among class II epitopes, 15 peptide-intrinsic features were significantly associated with ELISpot immunogenicity (GLM q-value < 0.05; Fig. 2C). Among the most significant positively associated included changes in the mutation position to a nonpolar amino acid (Mutated_position_change_of_NonPolar_feature), valine at position 1, tyrosine at position 6, and basic amino acid at position 2. The more significant negatively correlated feature was a change in the mutation position into a small amino acid (Mutated_position_change_of_Small_feature), which was positively correlated in class I epitopes. Negatively correlated features also included changes in the mutation position to a polar amino acid (Mutated_position_change_of_Polar_feature), and small/tiny amino acids at the mutated site. Among the significant features (Fig. 2D), we observed one cluster of closely intercorrelated features (Sig 1; n = 7, right-hand side of dendrogram), as well as a second cluster of loosely intercorrelated features (Sig 2; n = 8, left-hand side of dendrogram). With each respective tumor model defined as a binary variable (1 = true, 0 = false), the mean expression of cluster 1 features was significantly correlated with the B16F10 model and inversely correlated with MB49, whereas the mean expression of cluster 2 was significantly correlated with MB49 tumors (Supplementary Fig. S5; Spearman q-value <0.05). This combined with the greater burden of immunogenic class II neoantigens identified in these two models suggested relatively greater contribution of these models (particularly MB49) in the regression results. Using the same backward AIC stepwise regression approach described above, multivariable GLM regression was performed on six features. The resulting model indicated that 50.7% variation (P < 0.0001) in immunogenicity was explained by the prediction (Supplementary Fig. S6; Supplementary Table S2), with three significant features [tyrosine at positive 6 (Absolute_position_6_Y), valine at positive 1 (Absolute_position_1_V), and changes in the mutated position to a small amino acid (Mutated_position_change_of_Small_feature)] primarily driving the fit (P = 0.024, 0.002, and 1.6 × 10⁻⁵, respectively). As with class I epitopes, we did not observe significant differences in immunogenicity among predicted in-frame (n = 59) and out-of-frame (n = 9) epitopes (P > 0.05; Supplementary Fig. S3).

Our analysis of class I epitopes using multivariable GLM suggested an optimized multivariable model may be able to discriminate between high and low immunogenicity peptides (Fig. 2E). With the goal of designing a predictive model for neoantigen and mHA immunogenicity, we split our class I epitope database into an exploration (2/3 of epitopes, n = 141) and validation (1/3 of epitopes, n = 69) set (Fig. 3A). Class II modeling was not attempted due to the low number of epitopes available within our database (n = 68). In order to reduce noise within our model, we collapsed ELISpot scores with absolute values less than or equal to the absolute value of the most negative count (–53 spots) to zero. This was performed because we were not focused on the ability of the model to characterize exact immunogenicity values within the low immunogenicity range. Within the exploration set, we used a 10,000-fold cross-validation (2/3rd random resampling) approach, which demonstrated that only gradient boosting consistently performed better than chance. Our final gradient boosting algorithm contained seven predictive features: valine at position 1 (Absolute_position_1_V), valine at the last position (Last_position_V), small amino acids at the last position (Last_position_Small), basic amino acids of the reference sequence at the mutated position (Reference_AA_at_mutated_position_Basic), changes in the mutated position to a small amino acid (Mutated_position_change_of_Small_feature), lysine at relative site 1 ("Relative_site_1_K"), and presence of valine within the first 3 positions (First_three_AA_V). The class I validation set was run through this final gradient boosting algorithm, demonstrating significantly accurate performance when comparing the linear fit between the actual immunogenicity by ELISpot and the predicted immunogenicity by modeling (P = 0.01893, coefficient = 0.30; Fig. 3B). This model provided a high negative predictive value (83.6% predicted values <53), ideal in the setting for filtering out a large pool of predicted tumor antigens in order to select epitopes for therapeutic targeting.

To test whether our final algorithm increased the likelihood of identifying therapeutically relevant, immunogenic epitopes for antitumor vaccine responses, we used two tumor models within our validation set: BBN963 basal-like bladder cancer neoantigens (epithelial tumor) and P815 mastocytoma GvL mHA (hematopietic tumor). Epitopes were binned into predicted high immunogenicity (top quartile) and predicted low immunogenicity (bottom quartile) groups for comparison of relative efficacy (Supplementary Table S3). In BBN963 tumors, three predicted high (B2: VALLPSVML; C2: VSLTLFSSWL; A5: SNVMQLLL) and two predicted low (B5: ETLLNSATI; B12: MISRNRHTL) immunogenicity neoantigens were identified. Animals were vaccinated with 30 μg of one peptide (or no-peptide control) alongside 50 μg poly(I:C) as adjuvant, challenged with tumor 12 days after vaccination, then given a 30 μg peptide boost on day 21 after the initial vaccination (Fig. 3C). Animals vaccinated with predicted a high immunogenicity peptide survived longer than those vaccinated with either predicted low immunogenicity peptide (P = 0.0006) or no-peptide control (P = 0.0031; Fig. 3D; Supplementary Fig. S7A). In contrast, no significant difference in survival was observed between predicted low immunogenicity peptide and no-peptide control groups (P = 0.9674). We additionally observed better control of tumor size in animals vaccinated with predicted high immunogenicity peptide (Supplementary Fig. S7B–S7D).

In P815 tumors, two predicted high (AFQRVTCTTL and QYSSANDWTV) and three predicted low (HYAANEWI, KFFPNCIFL, and LYISPNPEVL) immunogenicity GvL mHAs were identified. BALB/c donor animals were vaccinated with a pool of predicted high or low immunogenicity peptides (100 μg each peptide) or no-peptide control, with 50 μg poly(I:C) as adjuvant on days 0 and 7. DBA/2 recipient animals were lethally irradiated on day 13; transplanted with 3 × 10⁶ BALB/c T cells, 3 × 10⁶ BALC/c bone marrow cells, and 3 × 10⁵ P815 tumor cells on day 14; and finally given a third booster vaccine on day 21 (Fig. 3E). Animals given predicted high immunogenicity T cells survived longer than those given predicted low immunogenicity T cells (median survival 44.5 and 28 days, respectively), both of which survived longer than no-peptide control T cells (median survival 19 days; Fig. 3F). Additionally, we observed significantly lower tumor burden in high immunogenicity versus low immunogenicity peptide vaccinated animals by luciferase imaging on day 26 (P < 0.05, Mann–Whitney U test; Supplementary Figs. S8 and S9). All groups receiving donor T cells demonstrated measurable GvHD clinically after transplant, without significant differences in weight loss or GvHD clinical scores between groups up to day 30 (Supplementary Fig. S10). In summary, these experiments demonstrated the in vivo biological relevance of our immunogenicity prediction model, with significant differences observed between predicted high and low immunogenicity epitopes.

Although the immunogenicity prediction algorithm was designed and validated in mice, we hypothesized that similar rules of immunogenicity may exist among human neoantigens. To study this, we ran previously predicted MHC I neoantigens from TCGA through our machine-learning algorithm, generating immunogenicity scores for each epitope (17). As expected, we observed a correlation between the number of HINs identified by our model (>95th percentile) with number of total neoantigens (Pearson correlation P < 0.0001; Supplementary Fig. S11). Therefore, we performed regression studies between HIN count and immune features without controlling for total neoantigen burden. We observed significant association between HIN count and IFNγ, cytotoxicity, CD8⁺ T-cell and total T-cell, and B-cell IGS among the data set (not controlling for cancer type; Fig. 4A). Assessing these associations individually by tumor type, we observed that the most significant associations were encompassed by the colon (COAD) and lung (LUAD) adenocarcinoma cancer types (Fig. 4B). Within COAD, there was a positive association between HIN count and many T-cell and cytotoxicity signatures. To identify potential drivers of this pattern, we looked for the association between HIN count and MSI status, observing that MSI-high COAD tumors had significantly higher HIN counts (Fig. 4C; Supplementary Fig. S12).

In contrast, LUAD demonstrated a negative association with signatures of anti–PD-1 responsiveness and several immune cell signatures. Regression analysis between these negative IGS features and LUAD oncogene/tumor suppressor copy numbers demonstrated preferential association with MYC copy number (q-value < 0.05; Fig. 4D; Supplementary Fig. S13). To demonstrate that MYC amplification provided a protumorigenic signal in LUAD, we observed significantly greater expression of genes corresponding with cell-cycle gene patterns (Supplementary Fig. S14A, S14C, and S14D), as well as enrichment of downstream genes involved in the MYC pathway (Supplementary Fig. S14B) among MYC-amplified tumors. This increased proliferation pattern was additionally associated with decreased sharing of T-cell receptor sequences from tumor-infiltrating T cells in MYC-amplified tumors, suggesting a potentially altered antitumor immune response (Supplementary Fig. S14E). We did not find HIN count to correlate with MYC copy number (Pearson P > 0.5), suggesting tumor immunogenicity burden and MYC target expression may be independent predictors for immune exclusion and checkpoint inhibition resistance.

Through the design of our neoantigen prediction algorithm, we considered predicted tumor epitopes across all open reading frames. As such, subsets of our predicted neoantigens were frameshifted epitopes that contained a mutation. We hypothesized that through mechanisms such as novel splice variants and ribosomal dysfunction, a subset of these out-of-frame predicted antigens could arise in the tumor, allowing for a viable target with greater heterogeneity from self-antigen. Indeed, two of the predicted high immunogenicity neoantigens used in our BBN963 vaccine studies were predicted to be out-of-frame (B2 and C2), although still providing therapeutic efficacy over predicted low immunogenicity and no-peptide controls. One of these antigens (B2) demonstrated computational evidence of translation in the out-of-frame context using two de novo transcriptome assemblers, Trinity (35) and StringTie (34), whereby the presence of a 5′ start codon was identified with no intervening stop codon up to the B2 antigen site.

With therapeutic and computational evidence in favor of B2 antigen-mediated antitumor immune responses against BBN963 tumors, we next confirmed the presence of B2/MHC tetramer–specific CD8⁺ T cells infiltrating within the tumor of BBN963-bearing animals, suggesting an antigen-driven T-cell response (Fig. 5A; Supplementary Fig. S15). Tetramer sorting and peptide-pulsed dendritic cell cocultures of B2-specific T cells demonstrated approximately 40-fold expansion of T cells within 10 days (<5 × 10⁵ to >2 × 10⁶), with maintenance of a B2-enriched population (Supplementary Fig. S16). Using a flow-cytometric cytotoxicity assay, coculture of B2-specific T cells with the BBN963 cell line demonstrated >1.7× increase in killing of target cells (15.25%) compared with the OTI T-cell irrelevant control (8.85%). Neither B2-specific nor OTI T cells demonstrated killing of irrelevant splenocytes control cells (1.1% and 0.85%, respectively; Fig. 5B; Supplementary Fig. S17). B2-specific T-cell killing of BBN963 over that of OTI T-cell controls was additionally confirmed using a ⁵¹Cr-release cytotoxicity assay (Fig. 5C). Altogether, these results suggested the presence of a cytotoxic CD8⁺ T-cell response against the out-of-frame B2 neoantigen in BBN963.

The study presented here addressed two unanswered questions regarding tumor antigens: (i) what features of a tumor antigen sequence are associated with immunogenicity and (ii) can inclusion of noncanonical, out-of-frame epitopes provide viable targets for antitumor therapeutic vaccination? We demonstrated that peptide-intrinsic features of predicted tumor antigens could discriminate epitopes with therapeutic efficacy, and that inclusion of out-of-frame epitopes among this pool could provide antitumor immunity against these alternative antigens. We showed that reading frame was not a significant determinant for immunogenicity (at least among peptides with predicted binding affinity <500 nmol/L), and that exclusion of frame-filtering could identify out-of-frame epitopes with therapeutic antitumor, cytotoxic activity. Although the optimal rules for immunogenicity may differ between in-frame and out-of-frame tumor antigens, our relatively limited training set was underpowered to discriminate between these two classes. As such, future studies should be performed to address the biological differences between in-frame and out-of-frame tumor antigens, and what methods can most optimally identify clinically relevant epitopes in each class.

Our analysis of class I epitopes demonstrated similar trends in expression of features significantly associated with immunogenicity, as well as potential generalizability of our final gradient boosting algorithm for human MHC. We were unable to demonstrate here whether MHC haplotype may influence immunogenicity prediction, given our murine models were limited to two haplotypes. That said, it may be the case that certain features may significantly impact immunogenicity in a way that is conserved across haplotypes. A potential bias in our analysis is the variation in distribution of ELISpot scores among various models. This variation is likely a product of both our selection process (i.e., peptide selection based on a threshold predicted MHC binding affinity) as well as the number of predicted epitopes available for screening in each model. As methodology for antigen prediction and validation was conserved across all models, as well as biological validation performed across two independent tumor models, we do not believe there to be significant underlying biological differences among epitopes identified between different tumor models.

Despite the increased interest in neoantigen-based therapeutic tumor vaccine therapy, prediction algorithms capable of directly predicting for neoantigen immunogenicity are lacking (50), and no neoantigen immunogenicity predictor trained specifically on tumor antigen data exists. Current neoantigen immunogenicity predictors are instead trained on databases containing immunogenicity scores from all potential MHC binding epitopes, of which the biology may not closely match that of mutation-derived tumor antigens. An example of this biological disparity is observed in the vastly different immune response rates between neoantigens and tumor-specific endogenous retroviral epitopes (51), whereby the concept of a “self” and “non-self” antigen is not considered as a feature for immunogenicity prediction among current algorithms. Training our model specifically on “self” tumor antigenic sources instead allowed for greater specificity of selection for peptide-intrinsic features that correlated with ex vivo–validated IFNγ release scores. Our final predictive algorithm demonstrated capacity to select for therapeutically relevant epitopes, as observed in our treatment studies where predicted high immunogenicity peptides controlled tumor burden significantly better than predicted low immunogenicity peptides and no-peptide control groups. This model demonstrated a high true-negative rate, which is ideal in the setting of filtering out many weakly immunogenic epitopes to select for a small pool of targets for therapy.

Validation experiments in BBN963 and P815 models were performed as a combination of prophylactic and therapeutic vaccines, rather than strictly treating animals after tumor inoculation. This method was selected due to the intrinsically low efficacy of free-peptide vaccines, whereby differences in therapeutic efficacy may not be observed between predicted high and low immunogenicity antigens (52). As such, although these experiments provided evidence for the biological relevance of our computational model, development of more robust therapeutic vaccine platforms are still necessary for improving response rates to peptide-based, tumor-specific antigen vaccines. From these studies, we observed that predicted high immunogenicity peptides had greater benefit in the therapeutic vaccine setting than predicted low immunogenicity peptides. As such, we reasoned that although HIN count and total neoantigen burden were correlated, the most HINs were the key drivers of immunity. Thus, we performed regression studies between HIN count and immune features without controlling for total neoantigen burden. Analysis of human neoantigen data from TCGA demonstrated association between presence of HINs with features of immune response in colon and LUAD. Athough the association between HIN count and immune gene signature expression, as well as MSI status, in COAD agreed with the classic view of a tumor antigen–driven immune response, the negative association with immune features (including signatures of anti–PD-1 responsiveness) in LUAD is less clear. Jerby-Arnon and colleagues demonstrated an association between resistance to immune-checkpoint inhibition and MYC target expression (53). As such, we initially hypothesized that MYC target expression may be the common driver for immune exclusion, anti–PD-1 nonresponsiveness, as well as high HIN burden. However, MYC copy number did not correlate with HIN count, suggesting MYC expression and HIN burden are independent pathways of immune exclusion and checkpoint inhibitor resistance in LUAD. Further studies are necessary to more closely examine the relationship between tumor immunogenic neoantigen burden and immune features, elucidating why higher immunogenicity burden is unexpectedly negatively associated with IGS patterns.

Currently, it is not well understood what frequency of tumor antigens arises from conventional in-frame antigens versus nonconventional antigenic sources, such as retroviral/retrotransponson expression, intron expression, and out-of-frame translation. A study from Laumont and colleagues suggests that noncoding regions are the main source of tumor-specific antigens in acute lymphoblastic leukemia patient samples (20), providing evidence that current methods for neoantigen prediction may be limited by filtering for in-frame exon regions. Laumont and colleagues used an RNA-based screening approach whereby k-mers derived from tumor RNA-seq reads were directly screened against matched-normal RNA k-mers, keeping only tumor-specific regions. Compared with conventional exome-based TSA calling, this RNA-based screening approach allowed for identification of a broader repertoire of epitopes, consistent with the increased frequency of noncanonical TSAs identified by Laumont and colleagues compared with this current study. Although Laumont and colleagues used a mass spectrometry approach to confirm expression of out-of-frame epitopes, no computational methods have been used to identify such noncanonical epitopes. As such, our study relies upon a naïve, nonbiased approach for screening out-of-frame antigens, whereby we combined conventional exome-based SNP antigen calling with identification of potential epitopes across all open reading frames. We demonstrated that the frame of an epitope did not associate with immunogenicity, but inclusion of out-of-frame epitopes could provide therapeutic benefit. This analysis highlighted how some proportion of SNV neoantigens predicted to be out-of-frame may still maintain expression and capacity to trigger a cytotoxic T-cell response against the tumor. If such antigens are further confirmed in human cancers, there are important implications that will need to be addressed: (i) whether the biology and immunogenicity of these out-of-frame “SNV” antigens more closely mirrors that of classic SNV neoantigens or whether they are instead more similar to INDEL-derived neoantigens or alternative neoantigens, such as tumor-specific endogenous retroviral antigens; and (ii) if reading frame filters should be applied to current neoantigen calling algorithms in order to most effectively capture the targetable antigen landscape of a tumor.

J.S. Serody reports receiving commercial research grants from Merck and GlaxoSmithKline and is a consultant/advisory board member for PIQUE Therapeutics. B.G. Vincent is cofounder of and has ownership interest (including patents) in Select Immunogenomics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: C.C. Smith, S. Chai, B.G. Vincent

Development of methodology: C.C. Smith, S. Chai, J.S. Parker, B.G. Vincent

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C.C. Smith, A.R. Washington, E. Landoni, K. Field, J. Garness, B. Savoldo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C.C. Smith, A.R. Washington, S.J. Lee, K. Field, S.R. Selitsky, J.S. Parker, B.G. Vincent

Writing, review, and/or revision of the manuscript: C.C. Smith, S. Chai, K. Field, L.M. Bixby, S.R. Selitsky, J.S. Parker, B. Savoldo, J.S. Serody, B.G. Vincent

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C.C. Smith, A.R. Washington, K. Field, L.M. Bixby, B.G. Vincent

Study supervision: C.C. Smith, B.G. Vincent

This work was supported by NIH grants F30 CA225136 (C.C. Smith), U54 CA198999 (J.S. Serody), and P50 CA058223 (J.S. Serody), as well as the UNC University Cancer Research Fund (B.G. Vincent) and the Susan G. Komen Career Catalyst Research Grant (B.G. Vincent).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.