Current tumor neoantigen calling algorithms primarily rely on epitope/major histocompatibility complex (MHC) binding affinity predictions to rank and select for potential epitope targets. These algorithms do not predict for epitope immunogenicity using approaches modeled from tumor-specific antigen data. Here, we describe peptide-intrinsic biochemical features associated with neoantigen and minor histocompatibility mismatch antigen immunogenicity and present a gradient boosting algorithm for predicting tumor antigen immunogenicity. This algorithm was validated in two murine tumor models and demonstrated the capacity to select for therapeutically active antigens. Immune correlates of neoantigen immunogenicity were studied in a pan-cancer data set from The Cancer Genome Atlas and demonstrated an association between expression of immunogenic neoantigens and immunity in colon and lung adenocarcinomas. Lastly, we present evidence for expression of an out-of-frame neoantigen that was capable of driving antitumor cytotoxic T-cell responses. With the growing clinical importance of tumor vaccine therapies, our approach may allow for better selection of therapeutically relevant tumor-specific antigens, including nonclassic out-of-frame antigens capable of driving antitumor immunity.

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.

Cell lines

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).

Animal studies

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.), 105 cells; BBN963: flank s.c., 107 cells; MB49: flank s.c., 105, UPPL1541: Flank s.c., 106, T11: mammary fat pad intradermal, 104, P815: tail-vein intravenous, 3 × 105. 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.

Tissue dissociation

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 and mHA prediction

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 mm3 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 CO2 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 (105 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).

Figure 1.

Summary of tumor antigen prediction and identification of peptide-intrinsic features. A, Number of MHC class I and II neoantigens/mHA per tumor model contained within the study. B, Schematic of neoantigen/mHA prediction and ELISpot validation workflow. C, Summary of the major classes of peptide-intrinsic features identified for each antigen, including amino acid sequence and characteristics at (I) each absolute position, (II) each relative site, (III) the mutation position, and (IV) the start, middle, and end of each peptide. Red boxes around columns in (I) demonstrate each absolute position and in (IV) demonstrate the start, middle, and end distinctions. Red lettering in (III) provides an example for SNV mutation site between reference and antigenic sequences.

Figure 1.

Summary of tumor antigen prediction and identification of peptide-intrinsic features. A, Number of MHC class I and II neoantigens/mHA per tumor model contained within the study. B, Schematic of neoantigen/mHA prediction and ELISpot validation workflow. C, Summary of the major classes of peptide-intrinsic features identified for each antigen, including amino acid sequence and characteristics at (I) each absolute position, (II) each relative site, (III) the mutation position, and (IV) the start, middle, and end of each peptide. Red boxes around columns in (I) demonstrate each absolute position and in (IV) demonstrate the start, middle, and end distinctions. Red lettering in (III) provides an example for SNV mutation site between reference and antigenic sequences.

Close modal

Vaccine/ELISpot screening

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 CO2 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 × 105 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.

Computational analysis

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.

Peptide treatment studies

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 × 107 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 mm2 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 × 106 T cells, 3 × 106 bone marrow cells, and 3 × 105 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.

Tetramer studies

Peptide/MHC tetramer and cell-surface protein staining were performed as described previously (37). Briefly, viable, single-cell suspensions derived from tumors (approximately 107 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.

Ex vivo T-cell expansion and cytotoxicity assays

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 106 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% CO2. 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 × 105 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% CO2 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 51Cr release assay (39). In brief, 5 × 103 51Cr-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 51Cr release, respectively. The mean percentage of specific lysis of triplicate wells was calculated as follows: [(test counts − spontaneous counts)/(maximum counts − spontaneous counts)] × 100%.

TCGA data analyses

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 log2 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

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.

Correlates of immunogenicity in class I MHC epitopes

Neoantigens and mHA were predicted in six tumor models (B16F10, BBN963, MB49, UPPL1541, P815, and T11), allowing us to characterize neoantigens in the H2b and H2d 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.

Figure 2.

Linear regression analysis between peptide-intrinsic features and tumor antigen immunogenicity. A and C, Volcano plots representing the GLM coefficient (x-axis) and −log10(q-value) (y-axis) for each peptide-intrinsic feature as a predictor for immunogenicity in class I (A) and class II (C) neoantigens/mHA. Dashed line represents q-value = 0.05. Spot color represents −log10(q-value) magnitude, and size represents magnitude of the coefficient. B and D, Heat map representing Spearman correlations between each significantly correlated feature from A and C for class I (B) and class II (D) neoantigens/mHA, respectively. Colored cells represent significantly correlated features (q < 0.05), with magnitude of the correlation coefficient represented by color. E, ELISpot-derived immunogenicity scores for class I neoantigens/mHA classified as predicted high (>100) or low (<100) immunogenic by multivariable GLM regression. Data represent median (middle line), with boxes encompassing the 25th to 75th percentile, whiskers encompassing 1.5 × the interquartile range from the box, and independent values shown by dots. Statistics performed with Mann–Whitney U test, with significance defined as P < 0.05.

Figure 2.

Linear regression analysis between peptide-intrinsic features and tumor antigen immunogenicity. A and C, Volcano plots representing the GLM coefficient (x-axis) and −log10(q-value) (y-axis) for each peptide-intrinsic feature as a predictor for immunogenicity in class I (A) and class II (C) neoantigens/mHA. Dashed line represents q-value = 0.05. Spot color represents −log10(q-value) magnitude, and size represents magnitude of the coefficient. B and D, Heat map representing Spearman correlations between each significantly correlated feature from A and C for class I (B) and class II (D) neoantigens/mHA, respectively. Colored cells represent significantly correlated features (q < 0.05), with magnitude of the correlation coefficient represented by color. E, ELISpot-derived immunogenicity scores for class I neoantigens/mHA classified as predicted high (>100) or low (<100) immunogenic by multivariable GLM regression. Data represent median (middle line), with boxes encompassing the 25th to 75th percentile, whiskers encompassing 1.5 × the interquartile range from the box, and independent values shown by dots. Statistics performed with Mann–Whitney U test, with significance defined as P < 0.05.

Close modal

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 Hb and Hd 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).

Correlates of immunogenicity in class II MHC epitopes

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−5, 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).

Machine-learning algorithm for immunogenicity prediction in class I MHC epitopes

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.

Figure 3.

Performance and validation of the gradient boosting model (GBM) approach for predicting neoantigen/mHA immunogenicity. A, Schema of the cross-validation approach used for GBM building. B, Performance of the final GBM in the validation set, showing actual (x-axis) versus predicted (y-axis) immunogenicity scores. Size of each point represents number of antigens at each coordinate. Red line represents the line of best fit, with P value of fit shown above the graph. C and E, Schema for in vivo validation experiments, with tumor vaccine studies performed in BBN963 basal-like bladder cancer (C) and the P815 mastocytoma syngeneic transplant model (E). D and F, Kaplan–Meier survival curves for animals bearing BBN963 basal-like bladder cancer (D) and P815 mastocytoma syngeneic transplants (F). Animals treated with predicted high (red) or low (blue) immunogenicity antigens, no-peptide control (black), or bone marrow–only control (gray). Data in D represent two independent experiments. Data in F represent one independent experiment. Statistics performed with log-ranked testing (**, P < 0.01; ***, P < 0.001). n.s., not significant.

Figure 3.

Performance and validation of the gradient boosting model (GBM) approach for predicting neoantigen/mHA immunogenicity. A, Schema of the cross-validation approach used for GBM building. B, Performance of the final GBM in the validation set, showing actual (x-axis) versus predicted (y-axis) immunogenicity scores. Size of each point represents number of antigens at each coordinate. Red line represents the line of best fit, with P value of fit shown above the graph. C and E, Schema for in vivo validation experiments, with tumor vaccine studies performed in BBN963 basal-like bladder cancer (C) and the P815 mastocytoma syngeneic transplant model (E). D and F, Kaplan–Meier survival curves for animals bearing BBN963 basal-like bladder cancer (D) and P815 mastocytoma syngeneic transplants (F). Animals treated with predicted high (red) or low (blue) immunogenicity antigens, no-peptide control (black), or bone marrow–only control (gray). Data in D represent two independent experiments. Data in F represent one independent experiment. Statistics performed with log-ranked testing (**, P < 0.01; ***, P < 0.001). n.s., not significant.

Close modal

In vivo validation of the class I immunogenicity prediction model

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 × 106 BALB/c T cells, 3 × 106 BALC/c bone marrow cells, and 3 × 105 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.

Correlates of predicted immunogenicity in human class I 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).

Figure 4.

Correlative analysis of predicted neoantigen immunogenicity in TCGA human data sets. A, Volcano plot representing GLM coefficient (x-axis) and −log10(q-value) (y-axis) between numbers of HINs and IGS in TCGA pan-cancer data sets. B, Heat map representing GLM regression between numbers of HINs and IGS for each TCGA cancer subset. Color represents direction of coefficient (red: positive; blue: negative), and shade represents −log10(q-value) magnitude. C, Number of HINs (x-axis) versus MSI score (y-axis) for a TCGA colorectal carcinoma (COAD) data set. D, Volcano plot representing GLM coefficient (x-axis) and −log10(q-value) (y-axis) between average expression of IGS in B with significantly negative association with HIN burden and oncogene/tumor suppressor copy numbers in a TCGA LUAD data set. A and D, Dashed line represents q-value = 0.05.

Figure 4.

Correlative analysis of predicted neoantigen immunogenicity in TCGA human data sets. A, Volcano plot representing GLM coefficient (x-axis) and −log10(q-value) (y-axis) between numbers of HINs and IGS in TCGA pan-cancer data sets. B, Heat map representing GLM regression between numbers of HINs and IGS for each TCGA cancer subset. Color represents direction of coefficient (red: positive; blue: negative), and shade represents −log10(q-value) magnitude. C, Number of HINs (x-axis) versus MSI score (y-axis) for a TCGA colorectal carcinoma (COAD) data set. D, Volcano plot representing GLM coefficient (x-axis) and −log10(q-value) (y-axis) between average expression of IGS in B with significantly negative association with HIN burden and oncogene/tumor suppressor copy numbers in a TCGA LUAD data set. A and D, Dashed line represents q-value = 0.05.

Close modal

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.

Out-of-frame neoantigen epitopes promote antitumor immunity

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 × 105 to >2 × 106), 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 51Cr-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.

Figure 5.

Analysis of out-of-frame epitope B2-specific T cells. A, Percentage of B2 tetramer–positive (left) versus irrelevant SIINFEKL-tetramer control (right) BBN963 tumor-infiltrating CD8+ T cells. Statistics performed with Mann–Whitney U test (*, P < 0.05). B, Four-hour flow-cytometric cytotoxicity assay comparing percentage of specific killing in 1:1 cocultures of B2 tetramer–specific T cells or OTI irrelevant T-cell controls versus BBN963 target or irrelevant splenocyte target control. Percentage of killing represents spontaneous target death background subtracted values. Statistics not performed for B due to n = 2 sample size across all groups. C, Four-hour chromium-51 release assay comparing percentage of specific killing in cocultures of B2 tetramer–specific T cells or OTI irrelevant T-cell controls versus BBN963 targets at 10:1 and 5:1 effector-to-target ratios. A–C, Error bars represent mean ± standard deviation. Data from each graph represent one independent experiment, respectively. Statistics performed with Welch t test (**, P < 0.01; *, P < 0.05).

Figure 5.

Analysis of out-of-frame epitope B2-specific T cells. A, Percentage of B2 tetramer–positive (left) versus irrelevant SIINFEKL-tetramer control (right) BBN963 tumor-infiltrating CD8+ T cells. Statistics performed with Mann–Whitney U test (*, P < 0.05). B, Four-hour flow-cytometric cytotoxicity assay comparing percentage of specific killing in 1:1 cocultures of B2 tetramer–specific T cells or OTI irrelevant T-cell controls versus BBN963 target or irrelevant splenocyte target control. Percentage of killing represents spontaneous target death background subtracted values. Statistics not performed for B due to n = 2 sample size across all groups. C, Four-hour chromium-51 release assay comparing percentage of specific killing in cocultures of B2 tetramer–specific T cells or OTI irrelevant T-cell controls versus BBN963 targets at 10:1 and 5:1 effector-to-target ratios. A–C, Error bars represent mean ± standard deviation. Data from each graph represent one independent experiment, respectively. Statistics performed with Welch t test (**, P < 0.01; *, P < 0.05).

Close modal

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.

1.
Gubin
MM
,
Artyomov
MN
,
Mardis
ER
,
Schreiber
RD
. 
Tumor neoantigens: building a framework for personalized cancer immunotherapy
.
J Clin Invest
2015
;
125
:
3413
21
.
2.
Castle
JC
,
Kreiter
S
,
Diekmann
J
,
Löwer
M
,
Van De Roemer
N
,
De Graaf
J
, et al
Exploiting the mutanome for tumor vaccination
.
Cancer Res
2012
;
72
:
1081
91
.
3.
Kreiter
S
,
Vormehr
M
,
van de Roemer
N
,
Diken
M
,
Löwer
M
,
Diekmann
J
, et al
Mutant MHC class II epitopes drive therapeutic immune responses to cancer
.
Nature
2015
;
520
:
692
6
.
4.
Kranz
LM
,
Diken
M
,
Haas
H
,
Kreiter
S
,
Loquai
C
,
Reuter
KC
, et al
Systemic RNA delivery to dendritic cells exploits antiviral defence for cancer immunotherapy
.
Nature
2016
;
534
:
396
401
.
5.
Sahin
U
,
Derhovanessian
E
,
Miller
M
,
Kloke
B-P
,
Simon
P
,
Löwer
M
, et al
Personalized RNA mutanome vaccines mobilize poly-specific therapeutic immunity against cancer
.
Nature
2017
;
547
:
222
6
.
6.
Ott
PA
,
Hu
Z
,
Keskin
DB
,
Shukla
SA
,
Sun
J
,
Bozym
DJ
, et al
An immunogenic personal neoantigen vaccine for patients with melanoma
.
Nature
2017
;
547
:
217
21
.
7.
Hoof
I
,
Peters
B
,
Sidney
J
,
Pedersen
LE
,
Sette
A
,
Lund
O
, et al
NetMHCpan, a method for MHC class i binding prediction beyond humans
.
Immunogenetics
2009
;
61
:
1
13
.
8.
Nielsen
M
,
Andreatta
M
. 
NetMHCpan-3.0; improved prediction of binding to MHC class I molecules integrating information from multiple receptor and peptide length datasets
.
Genome Med
2016
;
8
:
33
.
9.
Jurtz
V
,
Paul
S
,
Andreatta
M
,
Marcatili
P
,
Peters
B
,
Nielsen
M
. 
NetMHCpan-4.0: improved peptide–MHC class I interaction predictions integrating eluted ligand and peptide binding affinity data
.
J Immunol
2017
;
199
:
3360
8
.
10.
Andreatta
M
,
Karosiene
E
,
Rasmussen
M
,
Stryhn
A
,
Buus
S
,
Nielsen
M
. 
Accurate pan-specific prediction of peptide-MHC class II binding affinity with improved binding core identification
.
Immunogenetics
2015
;
67
:
641
50
.
11.
Zhang
Q
,
Wang
P
,
Kim
Y
,
Haste-Andersen
P
,
Beaver
J
,
Bourne
PE
, et al
Immune epitope database analysis resource (IEDB-AR)
.
Nucleic Acids Res
2008
;
36
:
W513
8
.
12.
Kim
Y
,
Ponomarenko
J
,
Zhu
Z
,
Tamang
D
,
Wang
P
,
Greenbaum
J
, et al
Immune epitope database analysis resource
.
Nucleic Acids Res
2012
;
40
:
W525
30
.
13.
O'Donnell
TJ
,
Rubinsteyn
A
,
Bonsack
M
,
Riemer
AB
,
Laserson
U
,
Hammerbacher
J
. 
MHCflurry: open-source class I MHC binding affinity prediction
.
Cell Syst
2018
;
7
:
129
32
.
e4
.
14.
Karosiene
E
,
Lundegaard
C
,
Lund
O
,
Nielsen
M
. 
NetMHCcons: a consensus method for the major histocompatibility complex class I predictions
.
Immunogenetics
2012
;
64
:
177
86
.
15.
Saito
R
,
Smith
CC
,
Utsumi
T
,
Bixby
LM
,
Kardos
J
,
Wobker
SE
, et al
Molecular subtype-specific immunocompetent models of high-grade urothelial carcinoma reveal differential neoantigen expression and response to immunotherapy
.
Cancer Res
2018
;
78
:
3954
68
.
16.
Colli
LM
,
Machiela
MJ
,
Myers
TA
,
Jessop
L
,
Yu
K
,
Chanock
SJ
. 
Burden of nonsynonymous mutations among TCGA cancers and candidate immune checkpoint inhibitor responses
.
Cancer Res
2016
;
76
:
3767
72
.
17.
Thorsson
V
,
Gibbs
DL
,
Brown
SD
,
Wolf
D
,
Bortone
DS
,
Ou Yang
TH
, et al
The immune landscape of cancer
.
Immunity
2018
;
48
:
812
30
.
e14
.
18.
Jayasinghe
RG
,
Cao
S
,
Gao
Q
,
Wendl
MC
,
Vo
NS
,
Reynolds
SM
, et al
Systematic analysis of splice-site-creating mutations in cancer
.
Cell Rep
2018
;
23
:
270
81
.
e3
.
19.
Kahles
A
,
Van Lehmann
K
,
Toussaint
NC
,
Hüser
M
,
Stark
SG
,
Sachsenberg
T
, et al
Comprehensive analysis of alternative splicing across tumors from 8,705 patients
.
Cancer Cell
2018
;
34
:
211
24
.
e6
.
20.
Laumont
CM
,
Vincent
K
,
Hesnard
L
,
Audemard
É
,
Bonneil
É
,
Laverdure
J-P
, et al
Noncoding regions are the main source of targetable tumor-specific antigens
.
Sci Transl Med
2018
;
10
:
eaau5516
.
21.
Weiss
RB
,
Dunn
DM
,
Atkins
JF
,
Gesteland
RF
. 
Slippery runs, shifty stops, backward steps, and forward hops: -2, -1, +1, +2, +5, and +6 ribosomal frameshifting
.
Cold Spring Harb Symp Quant Biol
1987
;
52
:
687
93
.
22.
Saulquin
X
,
Scotet
E
,
Trautmann
L
,
Peyrat
MA
,
Halary
F
,
Bonneville
M
, et al
+1 Frameshifting as a novel mechanism to generate a cryptic cytotoxic T lymphocyte epitope derived from human interleukin 10
.
J Exp Med
2002
;
195
:
353
8
.
23.
Macejak
DG
,
Sarnow
P
. 
Internal initiation of translation mediated by the 5′ leader of a cellular mRNA
.
Nature
1991
;
353
:
90
4
.
24.
Bullock
TNJ
,
Patterson
AE
,
Franlin
LL
,
Notidis
E
,
Eisenlohr
LC
. 
Initiation codon scanthrough versus termination codon readthrough demonstrates strong potential for major histocompatibility complex class I–restricted cryptic epitope expression
.
J Exp Med
1997
;
186
:
1051
8
.
25.
Bullock
TN
. 
Ribosomal scanning past the primary initiation codon as a mechanism for expression of CTL epitopes encoded in alternative reading frames
.
J Exp Med
1996
;
184
:
1319
29
.
26.
Malarkannan
S
,
Horng
T
,
Shih
PP
,
Schwab
S
,
Shastri
N
. 
Presentation of out-of-frame peptide/MHC class I complexes by a novel translation initiation mechanism
.
Immunity
1999
;
10
:
681
90
.
27.
Van Den Eynde
BJ
,
Gaugler
B
,
Probst-Kepper
M
,
Michaux
L
,
Devuyst
O
,
Lorge
F
, et al
A new antigen recognized by cytolytic T lymphocytes on a human kidney tumor results from reverse strand transcription
.
J Exp Med
1999
;
190
:
1793
800
.
28.
Bruce
AG
,
Atkins
JF
,
Gesteland
RF
. 
tRNA anticodon replacement experiments show that ribosomal frameshifting can be caused by doublet decoding
.
Proc Natl Acad Sci U S A
1986
;
83
:
5062
6
.
29.
Bruce
DW
,
Stefanski
HE
,
Vincent
BG
,
Dant
TA
,
Reisdorf
S
,
Bommiasamy
H
, et al
Type 2 innate lymphoid cells treat and prevent acute gastrointestinal graft-versus-host disease
.
J Clin Invest
2017
;
127
:
1813
25
.
30.
Roberts
PJ
,
Usary
JE
,
Darr
DB
,
Dillon
PM
,
Pfefferle
AD
,
Whittle
MC
, et al
Combined PI3K/mTOR and MEK inhibition provides broad antitumor activity in faithful murine cancer models
.
Clin Cancer Res
2012
;
18
:
5290
303
.
31.
Cooke
KR
,
Kobzik
L
,
Martin
TR
,
Brewer
J
,
Delmonte
J
,
Crawford
JM
, et al
An experimental model of idiopathic pneumonia syndrome after bone marrow transplantation: I. The roles of minor H antigens and endotoxin
.
Blood
1996
;
88
:
3230
9
.
32.
Wilkerson
MD
,
Cabanski
CR
,
Sun
W
,
Hoadley
KA
,
Walter
V
,
Mose
LE
, et al
Integrated RNA and DNA sequencing improves mutation detection in low purity tumors
.
Nucleic Acids Res
2014
;
42
:
e107
.
33.
Rajasagi
M
,
Shukla
SA
,
Fritsch
EF
,
Keskin
DB
,
DeLuca
D
,
Carmona
E
, et al
Systematic identification of personal tumor-specific neoantigens in chronic lymphocytic leukemia
.
Blood
2014
;
124
:
453
62
.
34.
Pertea
M
,
Pertea
GM
,
Antonescu
CM
,
Chang
TC
,
Mendell
JT
,
Salzberg
SL
. 
StringTie enables improved reconstruction of a transcriptome from RNA-seq reads
.
Nat Biotechnol
2015
;
33
:
290
5
.
35.
Haas
BJ
,
Papanicolaou
A
,
Yassour
M
,
Grabherr
M
,
Blood
PD
,
Bowden
J
, et al
De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis
.
Nat Protoc
2013
;
8
:
1494
512
.
36.
Coghill
JM
,
Carlson
MJ
,
Panoskaltsis-Mortari
A
,
West
ML
,
Burgents
JE
,
Blazar
BR
, et al
Separation of graft-versus-host disease from graft-versus-leukemia responses by targeting CC-chemokine receptor 7 on donor T cells
.
Blood
2010
;
115
:
4914
22
.
37.
Dolton
G
,
Tungatt
K
,
Lloyd
A
,
Bianchi
V
,
Theaker
SM
,
Trimby
A
, et al
More tricks with tetramers: a practical guide to staining T cells with peptide-MHC multimers
.
Immunology
2015
;
146
:
11
22
.
38.
Wölfl
M
,
Greenberg
PD
. 
Antigen-specific activation and cytokine-facilitated expansion of naive, human CD8+T cells
.
Nat Protoc
2014
;
9
:
950
66
.
39.
Quintarelli
C
,
Dotti
G
,
De Angelis
B
,
Hoyos
V
,
Mims
M
,
Luciano
L
, et al
Cytotoxic T lymphocytes directed to the preferentially expressed antigen of melanoma (PRAME) target chronic myeloid leukemia
.
Blood
2008
;
112
:
1876
85
.
40.
Chan
KS
,
Espinosa
I
,
Chao
M
,
Wong
D
,
Ailles
L
,
Diehn
M
, et al
Identification, molecular characterization, clinical prognosis, and therapeutic targeting of human bladder tumor-initiating cells
.
Proc Natl Acad Sci U S A
2009
;
106
:
14016
21
.
41.
Prat
A
,
Parker
JS
,
Karginova
O
,
Fan
C
,
Livasy
C
,
Herschkowitz
JI
, et al
Phenotypic and molecular characterization of the claudin-low intrinsic subtype of breast cancer
.
Breast Cancer Res
2010
;
12
:
R68
.
42.
Iglesia
MD
,
Vincent
BG
,
Parker
JS
,
Hoadley
KA
,
Carey
LA
,
Perou
CM
, et al
Prognostic B-cell signatures using mRNA-seq in patients with subtype-specific breast and ovarian cancer
.
Clin Cancer Res
2014
;
20
:
3818
29
.
43.
Kardos
J
,
Chai
S
,
Mose
LE
,
Selitsky
SR
,
Krishnan
B
,
Saito
R
, et al
Claudin-low bladder tumors are immune infiltrated and actively immune suppressed
.
JCI Insight
2016
;
1
:
e85902
.
44.
Bindea
G
,
Mlecnik
B
,
Tosolini
M
,
Kirilovsky
A
,
Waldner
M
,
Obenauf
AC
, et al
Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer
.
Immunity
2013
;
39
:
782
95
.
45.
Love
MI
,
Huber
W
,
Anders
S
. 
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
46.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
Proc Natl Acad Sci U S A
2005
;
102
:
15545
50
.
47.
Ingenuity Systems
. 
Ingenuity Pathway Analysis
.
WwwIngenuityCom
. 
2013
;
5020
.
48.
Huang
DW
,
Sherman
BT
,
Tan
Q
,
Collins
JR
,
Alvord
WG
,
Roayaei
J
, et al
The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists
.
Genome Biol
2007
;
8
:
R183
.
49.
Bolotin
DA
,
Poslavsky
S
,
Mitrophanov
I
,
Shugay
M
,
Mamedov
IZ
,
Putintseva
E V
, et al
MiXCR: software for comprehensive adaptive immunity profiling
.
Nat Methods
2015
;
12
:
380
1
.
50.
Kim
S
,
Kim
HS
,
Kim
E
,
Lee
MG
,
Shin
EC
,
Paik
S
, et al
Neopepsee: accurate genome-level prediction of neoantigens by harnessing sequence and amino acid immunogenicity information
.
Ann Oncol
2018
;
29
:
1030
6
.
51.
Smith
CC
,
Beckermann
KE
,
Bortone
DS
,
Cubas
AA
,
Bixby
LM
,
Lee
SJ
, et al
Endogenous retroviral signatures predict immunotherapy response in clear cell renal cell carcinoma
.
J Clin Invest
2018
;
128
:
4804
20
.
52.
Irvine
DJ
,
Hanson
MC
,
Rakhra
K
,
Tokatlian
T
. 
Synthetic nanoparticles for vaccines and immunotherapy
.
Chem Rev
2015
;
115
:
11109
46
.
53.
Jerby-Arnon
L
,
Shah
P
,
Cuoco
MS
,
Rodman
C
,
Su
MJ
,
Melms
JC
, et al
A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade
.
Cell
2018
;
175
:
984
97
.
e24
.