Predicting HLA class II antigen presentation through integrated deep learning

Abstract

Accurate prediction of antigen presentation by human leukocyte antigen (HLA) class II molecules would be valuable for vaccine development and cancer immunotherapies. Current computational methods trained on in vitro binding data are limited by insufficient training data and algorithmic constraints. Here we describe MARIA (major histocompatibility complex analysis with recurrent integrated architecture; https://maria.stanford.edu/ ), a multimodal recurrent neural network for predicting the likelihood of antigen presentation from a gene of interest in the context of specific HLA class II alleles. In addition to in vitro binding measurements, MARIA is trained on peptide HLA ligand sequences identified by mass spectrometry, expression levels of antigen genes and protease cleavage signatures. Because it leverages these diverse training data and our improved machine learning framework, MARIA (area under the curve = 0.89-0.92) outperformed existing methods in validation datasets. Across independent cancer neoantigen studies, peptides with high MARIA scores are more likely to elicit strong CD4⁺ T cell responses. MARIA allows identification of immunogenic epitopes in diverse cancers and autoimmune disease.

Fig. 1 |. Rationale and framework for the development of a new method for prediction of HLA-II ligands.

a, Comparison of the number of unique HLA-DR ligands identified within two antigen presentation profiling studies^,, each exceeding all HLA-DR in vitro quantitative binding measurements from 239 previous studies within the IEDB (as of December 2018). b, Performance of NetMHCIIpan for discrimination of decoys from bona fide HLA-II ligands recovered by antigen presentation profiling. For each patient, NetMHCIIpan-predicted affinities and ranks were separately evaluated (x axis), and performance measured by ROC-AUC (y axis, dotted lines represent the median). NetMHCIIpan ranks (mean AUC = 0.68) slightly outperformed NetMHCIIpan binding affinities (mean AUC = 0.65, n = 18; two-tailed paired t test, P = 0.003; Supplementary Table 1). c, Limited sensitivity of NetMHCIIpan for classification of HLA-DR ligands. Depicted is the distribution of NetMHCIIpan ranks for all 6,063 peptides identified from the JeKo-1 cell line, where 22% of HLA-II ligands identified by MS had predicted values worse than the recommended NetMHCIIpan rank cut-off for binding (10%). d, In vitro binding assay results for HLA-II peptide ligands identified by MS but predicted by NetMHCIIpan not to bind HLA-II. Among ten such peptides predicted by NetMHCIIpan not to bind, nine were nevertheless confirmed to bind cognate HLA-DR alleles (04:03 and/or 07:01) by two independent flow cytometry experiments. Scatter plots depict binding of two exemplar FITC-conjugated peptides (x axis) to APC-conjugated HLA-DR proteins (y axis); remaining peptides are separately depicted in Supplementary Fig. 1. e, Training and evaluation scheme of MARIA, as a new machine learning framework for more accurate prediction of HLA-II ligands. Positive examples are HLA-II ligand peptide sequences directly identified by antigen presentation profiling of human cells and tissues by immunoprecipitation (i.p.) and MS, and negative examples are length-matched random human peptides (decoys). The model separately considers binding affinities estimated using in vitro binding data. Patient HLA-II allele or genotype and gene expression information are obtained from next-generation sequencing. A RNN integrates information and produces a predictor for HLA-II ligand presentation by minimizing training errors. Independent test sets determine the final performance of the model. See Supplementary Fig. 2 for detailed machine learning schemes.

Fig. 2 |. Features, model architecture and validation performance of MARIA.

a, Comparison of gene expression levels of HLA-DR ligands and non-ligands. Gene expression was estimated by RNA-seq for HLA-DR-presented genes, all protein-coding genes and non-presented protein-coding genes, respectively. HLA-DR ligand genes have significantly higher gene expression levels than the set of all protein-coding genes (n = 34,049, 23,165 and 19,464, respectively; **P < 1 × 10⁻⁵, Mann–Whitney U test). Some HLA-DR ligands (8.4%) had undetectable levels of expression; those in this set were enriched for extracellular protein (GO enrichment; Fisher’s exact test, P < 1 × 10⁻¹⁷). Violin curves represent the probability distribution function of gene expression, black boxes represent middle two quartiles and white dots represent the median. See Supplementary Fig. 3 for detailed analysis on HLA-I ligands and the predictive power of gene expression levels. b, Cleavage signature analysis for HLA-DR ligands. Frequencies of 20 amino acids at 6 positions upstream (−6 to −1) and downstream (+1 to +6) of HLA-DR ligands (n = 12,150) are compared to the background distribution (n = 23,218) to determine amino acid enrichment and depletion surrounding HLA-DR ligands. Colors of the heat map and sizes of the logo plot letters indicate fold change. The logo plot only includes statistically significant enrichment (P < 0.001, two-tailed independent t test by IceLogo). The minus symbol in the top row of the heat map indicates presented peptides that are located at the beginning or end of source protein sequences. See Supplementary Fig. 4 for cleavage signatures across different cell types. c, Workflow of MARIA for predicting HLA-DR ligand presentation score. Two separate models first calculate HLA-DR peptide binding scores and peptide cleavage scores. The neural network further integrates peptide sequence and estimated gene expression level with two scores, via a recurrent layer and merge layers, to generate a presentation score indicating likelihood of HLA-II presentation. d, Performance of MARIA and four alternative predictors on 10% of the held-out validation set (true MCL HLA-II ligands, n = 3,300; random human decoy peptides, n = 10,000; the same sample set is used in e and f). MARIA scores incorporating gene expression levels, peptide sequence, binding scores and cleavage scores outperformed methods using each of these features individually (DeLong test, P < 1 × 10⁻⁵; AUC = 0.92). See Supplementary Fig. 5 for detailed training data source and cross-validation performance. e, Comparison of model precision and specificity across a range of presented MCL HLA-DR peptide prevalences. Sensitivity for each model was controlled at 30% for all calculations, with corresponding specificity denoted adjacent to inset legend. The shaded areas represent the 95% confidence interval around the mean value, on the basis of tenfold cross-validation. f, Comparison of precision and recall for different models for predicting HLA-DR ligands using various types of training data. Precision was calculated assuming 1% prevalence of presented HLA-DR ligands. The shaded areas represent 95% confidence interval around the mean value (line), based on tenfold cross-validation.

Fig. 3 |. Benchmarking MARIA performance against existing binding-based methods with independent HLA-DR test sets.

a, Overlap and sequence motifs of two HLA-DR ligand sets identified from two monoallelic K562 cell lines. A proportion (31%) of peptides appeared in both the HLA-DRB1*01:01 (n = 2,430) and HLA-DRB1*04:04 (n = 2,072) ligand sets when considering substring matches. The sequence motifs with highest statistical significance (P < 1 × 10⁻⁷, multiple hypergeometric test implemented by MEME) are shown. For full potential motifs, see Supplementary Table 4. b, Performance of MARIA and six alternative methods when differentiating 1,361 K562 HLA-DRB1*01:01 ligands from 1,361 human decoys. MARIA outperformed the second-best method (SMM Align; DeLong test, P < 1 × 10⁻⁵). Limited by the IEDB Concensus3 package, only ligand sequences ≥15 amino acids are included in this comparison. c, Performance of MARIA and four alternative methods differentiating 2,032 K562 DRB1*04:04 ligands from 2,032 human decoys. MARIA achieved an AUC of 0.89 AUC as compared to an AUC of 0.56 for NetMHCIIpan. RNN and SNN trained on MCL ligands obtained AUC values of 0.83 and 0.78, respectively.

Fig. 4 |. MARIA trained on human HLA-DQ ligand peptides identified celiac-related gluten antigens.

a, Overlap and sequence motifs of two HLA-DQ ligand sets. A majority (65%) of peptides were present in both HLA-DQ2.2 (n = 7,374) and HLA-DQ2.5 (n = 4,249) ligand sets when considering substring matches. The sequence motifs with highest statistical significance (P < 1 × 10⁻⁷, multiple hypergeometric test implemented by MEME) are shown. Full potential motifs are included in Supplementary Table 4. b, Performance of MARIA trained on HLA-DQ2.2 ligand sequences and tested on a held-out human HLA-DQ2.2 peptide set (n = 650). MARIA was trained on 90% of the HLA-DQ2.2-associated peptide sequences shown in a. MARIA achieves an AUC of 0.89 when differentiating DQ2.2 ligands from length-matched decoys. By comparison, NetMHCIIpan percentiles obtained an AUC of 0.68. Dashed red lines indicate the 90th percentile, the default cut-off for NetMHCIIpan. See Supplementary Fig. 7 for detailed training schemes of MARIA for HLA-DQ2.2. c, Performance of MARIA and NetMHCIIpan when identifying immunogenic gluten peptide fragments (n = 69). MARIA trained on human DQ2.2 ligands identified 49% of HLA-DQ2.2-binding gluten peptides with 92% specificity. By comparison, NetMHCIIpan had 6% sensitivity and 88% specificity. Dashed red lines indicate the 90th percentile, the default cut-off for NetMHCIIpan. The x axes in b and c capture the percentiles for depicted bins, where higher percentiles reflect higher likelihood of presentation, by defining the percentiles as 100% minus the absolute rank reported by each method.

Fig. 5 |. MARIA identifies lymphoma immunoglobulin HLA-DR presentation hotspots in patients with MCL.

a, Correlation of MARIA-predicted and experimentally identified HLA-DR-presented immunoglobulin antigens. Eighteen MCL immunoglobulin sequences were analyzed by a version of MARIA trained on non-immunoglobulin HLA-DR ligands to determine the presentation hotspots (left, blue). The same 18 MCL samples were profiled with LC–MS/MS to determine the regions of immunoglobulin presented by HLA-DR (right, orange). Predicted and observed presentation hotspots were significantly correlated on both heavy chains and light chains (Spearman’s ρ of 0.63 and 0.55, P = 1 × 10⁻⁶⁵ and 7.5 × 10⁻¹⁹; n = 1,015 and 311, respectively). MARIA-predicted ligand numbers were normalized with the MS-identified maximum ligand numbers for visualization purposes. See Supplementary Fig. 8 for the presentation heat map predicted by NetMHCIIpan. b, MARIA identified HLA-DR presentation hotspots in the immunoglobulin heavy chain variable region (IGHV). MARIA-predicted HLA-DR-presented peptides from IGHV FR3 regions more than the other six regions across patients (P < 1 × 10⁻⁵, Mann–Whitney U test), consistent with MS findings (P < 1 × 10⁻⁵, Mann–Whitney U test). Each dot represents predicted or experimentally identified ligand coverage in a 15-amino-acid sliding window on the aligned IGHV sequence (n = 38 for the FR3 region and n = 87 for the non-FR3 regions).

Fig. 6 |. MARIA scores predict melanoma HLA-II-presented antigens and are associated with post-vaccine CD4⁺ T cell responses.

a, Performance of MARIA on an independent melanoma HLA-II ligand set3. MARIA trained on MCL ligands achieved an AUC of 0.89 when differentiating patient melanoma HLA-II peptides from length-matched decoys, as compared to NetMHCIIpan with an AUC of 0.64. Shuffling correct training labels diminished the prediction performance of MARIA, reducing its AUC to 0.53. b, Neoantigen gene expression in patients with melanoma is not associated with postvaccination CD4⁺ T cell responses. Personalized gene expression values were obtained from tumor RNA-seq in two personalized melanoma vaccine trials^,. In both trials, there is no difference in gene expression values between positive and negative vaccine candidates for their ex vivo CD4 cytokine release tests (n = 127 and 97; P = 0.49 and 0.50, two-tailed unpaired t test). NS, not significant. c,d, Post-vaccination CD4⁺ T cell responses are associated with MARIA scores. Peptide sequences from the same two clinical trials were scored with MARIA (c, n = 127 for Sahin et al., and d, n = 97 for Ott et al.). Each candidate was stratified into three categories on the basis of the highest MARIA percentile scores among 15-amino-acid oligomer sliding windows: low (<95th), medium (95-99.5th) and high (>99.5th). Dashed red lines indicate average response rates of the whole cohort. c, A majority (73%) of peptides in the MARIA high category elicited positive CD4⁺ T cell responses after vaccinations as compared to 26% in the low category and 47% in the medium category (χ² test, 2 degrees of freedom, P = 0.019). d, A higher proportion (38%) of peptides in the MARIA high category elicited positive CD4⁺ T cell responses after vaccinations as compared to the low category (8.1%) and the medium category (23%) (χ² test, 2 degrees of freedom, P = 0.023). See Supplementary Fig. 10 for detailed analysis on melanoma cancer vaccines. e, Relationship between MARIA percentile scores and CD4⁺ T cell responses to tumor-associated antigens across cancer types and studies. When considering seven different studies of CD4⁺ T cell immune responses to cancer-associated antigens (including this one)^,,– (rows), we identified immunogenic (positive; n = 27, rows 1-7) and non-immunogenic (negative; n = 494, row 8) peptides across diverse tumor types. Each of these 521 peptides (dots) were then tested by MARIA, allowing comparison of percentile scores (x axis, right) with immunogenicity (blue, immunogenic; green, non-immunogenic). As depicted by the summarized inset table, 74% of immunogenic peptides (20 of 27, blue) scored above the 95th MARIA percentile threshold. T_eff, effector T cells.