Integrated Germline and Somatic Features Reveal Divergent Immune Pathways Driving Response to Immune Checkpoint Blockade

The development of immune checkpoint blockade (ICB) therapeutics has shifted the cancer treatment paradigm, offering unprecedented hope for patients who once faced limited therapeutic options (1, 2). The remarkable successes of ICB, leading to complete remissions in some patients with advanced cancers have propelled this approach to the forefront of modern oncology (3). ICB is now a standard treatment in some tumor types; however, a substantial proportion of patients still fail to benefit while experiencing the side effects and costs of the therapeutics (4–6). Despite several landmark studies on biomarkers for immunotherapy response (7–9), identifying those patients who would effectively respond to immunotherapy remains a challenge (10).

Thus far, biomarkers have focused on measured characteristics of the tumor or the tumor immune microenvironment (TIME). Current FDA-approved biomarkers include tumor mutation burden (TMB), microsatellite instability status, and IHC staining of the tumor microenvironment to quantify PDL1 positivity (11). However, these predictors of response are imperfect and their application in clinical settings is not straightforward (12). Many more sophisticated measures of ICB response have also been proposed, including the potential immunogenicity of somatic mutations in the tumor (13, 14), measures of immunoediting such as the ratio of nonsynonymous:synonymous mutations of the immunopeptidome (15), evidence of impaired antigen presentation quantified from somatic copy number loss and mutation of MHC genes (16–18), and tumor clone phylogeny estimates as a proxy for intratumoral heterogeneity (19). Anagnostou and colleagues (20) successfully integrated somatic features such as these to predict ICB response using machine learning models with superior accuracy, suggesting nonlinear predictive models may capture additional biological complexity. These approaches show promise to improve over current FDA-approved measures, although performance gains have generally been modest.

More recent work has uncovered a role for germline genetic variation in influencing the characteristics of the TIME and ICB response. Although whole exome sequencing (WES) methods require a matched normal tissue as a background panel for somatic mutation detection (21), patient germline variation has largely been ignored in the development of predictive ICB modeling, even though germline variation has a considerable effect on adaptive immune traits (22–24). We reasoned that although individual common variants often have only a weak influence on traits, the sum of these variations could have a large impact on the TIME, as suggested by a study in which common germline variants were found to predict ICB responses independent of somatic biomarkers (25). With the exception of some rare germline variants, cancer often arises from mutagenic processes independent of host germline genetics (26).

In this study, we developed a machine learning framework that integrates both somatic and germline features into a unified model that aims to maximize the identification of patients who may benefit from ICB therapy. We used XGBoost for the model architecture as it has shown strong performance on limited training data, allows for nonlinear interactions among features, and is interpretable in that individual feature contributions to predictive performance can be quantified (27–29). A composite model using all features demonstrated superior performance across multiple independent test sets relative to predictors trained on germline or somatic features alone. Analysis of the composite model revealed feature interactions that contributed to model performance, the strongest of which occurred between MHC class-I (MHC-I) damage and a germline variant associated with increased infiltration of T-follicular helper cells. Further investigation of this interaction suggested an MHC-I–independent mechanism of ICB response associated with the MHC class-II (MHC-II) CD4⁺ T-cell axis in some patients. Grouping ICB responders by response type showed more durable ICB responses in the MHC-II–driven response axis. For the 34% of patients with RNA expression data, we also investigated characteristics of the TIME such as checkpoint expression, T-cell infiltration, and tertiary lymphoid structure (TLS) signatures (30–32). Overall, our results support the notion that nonlinear models using somatic and germline features together to predict ICB outcomes allow for the formulation of new hypotheses about biological mechanisms underlying the diversity of clinical responses to ICB.

Raw FASTQ files were obtained using SRA toolkit v2.9.6-1-ubuntu64 for the following immune checkpoint trials: Hugo and colleagues [SRA accession: SRP090294, SRP067938; Cancer: melanoma; n = 38; RNA sequencing (RNAseq) = 27], Van Allen and colleagues (SRA accession: SRP011540, Cancer: melanoma; n = 110; RNAseq = 40), Miao and colleagues (SRA accession: SRP128156, Cancer: clear cell renal carcinoma; n = 69; RNAseq = 33), Riaz and colleagues (SRA accession: SRP095809, SRP094781; Cancer: melanoma; n = 57; RNAseq = 46), Rizvi and colleagues (SRA accession: SRP064805, Cancer: non–small cell lung cancer; n = 35), Snyder and colleagues (SRA accession: SRP072934, Cancer: melanoma; n = 64), Liu and colleagues (SRA accession: PRJNA82747 Cancer: melanoma; n = 122; RNAseq = 122), and Cristescu and colleagues (SRA accession: PRJNA449580, Cancer: Melanoma, HSNCC, Urothelial; n = 213). Only pretreatment samples were utilized in this study. Across cohorts, a total of 708 ICB-treated patients were evaluated in this study. The Cancer Genome Atlas (TCGA; ref. 33) samples from similar tissues to the ICB cohorts were obtained from the TCGA data portal https://portal.gdc.cancer.gov/. Only patients with available genetic SNP, somatic, and RNAseq data were extracted (n = 3,377) from the following TCGA cohorts (LUAD, n = 522; KIRC, n = 537; SKCM, n = 470; HNSC, n = 528; BLCA, n = 412; KICH, n = 113; KIRP, n = 291; LUSC, n = 504).

FASTQ files were processed via an identical bioinformatics pipeline. DNA: Genomic reads were aligned to UCSC hg19 coordinates using BWA v0.7.17-r1188. Reads were sorted by SAMTOOLS v0.1.19, marked for duplicates with Picard Tools v2.12.3 and recalibrated with GATK v3.8-1-0. Germline variants were called from sorted BAM (Binary Alignment Map) files using DeepVariant v0.10.0-gpu. Somatic variants were obtained through the following additional steps. Aligned tumor/normal BAM files were submitted to standard Mutect2 somatic variant calling using GATK-4.1.3.0. First, BAM file formats were standardized using GATK-4.1.3.0 AddorReplaceReadGroups, then GATK-4.1.3.0 Mutect2 was used to call somatic variants using default settings (including the presence of a matched normal), the gnomAD v3.1 raw sites background SNP panel, and the Twist Exome Target bed file to limit variant calling to exonic regions. Potential somatic variants were filtered using GATK-4.1.3.0 FilterMutectCalls and only mutations with a filter flag of “PASS” were kept for subsequent analysis. Somatic mutations were further filtered to retain only those with a DNA allelic fraction >5%. The resulting variant call format (VCF) files were annotated by variant effect prediction using cache version 102_GRCh37 and default settings. RNA: When available, RNA FASTQ/BAM files were downloaded for 33 RCC and 240 patients with melanoma. BAM files were converted to FASTQ using bam2fq. Unpaired reads were removed using fastq pair. Paired reads were aligned with STAR v2.4.1 d to GRCh37 reference alignment. RSEM v1.2.21 was used for transcript quantification. Raw transcript counts were corrected for cohort-specific batch effects using ComBat before being transformed into transcript per million (TPM) values.

A set of 1,084 TIME-associated SNPs was sourced from Pagadala and colleagues (25). These SNPs were demonstrated in the aforementioned study to have significant associations with immune-related functions in TCGA and were successfully used to develop an earlier germline ICB response prediction model (25). Next, we filtered for SNPs present with a mutant allele fraction >0.05 in all studies, leaving 598 SNPs to run METAL (34) analysis with ICB response in the three training cohorts. METAL analysis calculates a single P-value for each SNP across the three training cohorts (Hugo and colleagues, Riaz and colleagues, and Snyder and colleagues) and indicates the direction of effect for each cohort. SNPs with an FDR of <0.25 and showing full agreement with direction of impact were included, resulting in 229 SNPs with a nominal ICB association. TCGA and discovery genotype processing was performed by Pagadala and colleagues and is described in detail in their methods. For this study, we obtained preprocessed genotype matrices for each of the cohorts examined.

TMB was defined as the sum of all nonsynonymous somatic coding mutations in each patient’s VCF file, including “protein coding,” “frameshift variant,” and “stop lost” mutations. To adjust for cohort-specific effects, TMB was transformed by the intra-cohort z-score before being included in the machine learning model. A similar convention is described in Vokes and colleagues (35).

A comprehensive list of immune escape–related genes was obtained from Zapata and colleagues (15). Somatic mutations with variant effect prediction impact annotations of “MODERATE” or “HIGH” were tallied from per patient VCF files. The final immune escape mutation counts were divided by each patient’s total TMB to generate a score reflecting disproportionate immune evasion—otherwise, the score is highly correlated with TMB.

A list of key antigen presentation pathway–related genes was obtained from MSigDB M1062, Reactome Antigen Presentation Folding Assembly, and Peptide Loading of Class-I MHC. All HLA genes were removed from this list as they are accounted for with better accuracy by HLA-specific tools and summarized in other features. Somatic mutations with an impact of “MODERATE” or “HIGH” were tallied from per patient VCF files. The resulting scores were divided by each patient’s total TMB to generate a score reflecting disproportionate damage to the antigen presentation pathway.

Intratumoral heterogeneity (ITH) and fraction of TMB subclonal both rely on accurate subclonal estimates, which are derived as follows. First, copy number calling was performed using CNVkit v0.9.10. A background panel of normals was constructed for each cohort separately using CNVkit reference to protect against batch effects. CNVkit batch was used to call copy number changes with each respective cohort’s matched background panel. We next used PureCN v2.6.4 (run via singularity image) with CNVkit-derived .cnr and .seg files, and Mutect2-derived filtered VCF files to generate purity and ploidy metrics to be used in subsequent subclone estimation. PureCN was run with default settings, repeat regions censored, and a random seed set to 123. Next, PyClone-VI v0.13.1 was run on mutation-specific integer copy number estimates derived from the CNVkit call (https://cnvkit.readthedocs.io/en/stable/heterogeneity.html) to estimate the clonal structure of the tumor. ITH was defined as the total number of subclones with at least five mutations (total range 0–11 subclones). Fraction of TMB subclonal was calculated by taking the total number of mutations belonging to small subclones (<5 mutations per subclone) and dividing it by the total number of mutations for each tumor. This generates an inverse estimate of clonal heterogeneity from ITH.

Immunoediting evaluates the ratio of nonsynonymous to synonymous mutations (dN/dS) in a tumor as a measure of selection (36). Immune dN/dS was adapted by Zapata and colleagues (15) in their toolkit SOPRANO (https://github.com/luisgls/SOPRANO) to calculate the immunoediting score for each patient using an hg19 reference and default settings. Essentially, this score derives from calculating dN/dS across all regions of the proteome predicted to bind the set of patient-specific MHC alleles (i.e., displayed for immune surveillance) and ranges from 0 to ∼5 with a score above one indicating a higher amount of nonsynonymous mutations to synonymous ones.

MHC-I damage was defined as the union of POLYSOLVER (16) and LOHHLA (17) results. First, class-I HLA alleles were genotyped via POLYSOLVER [see Patient Harmonic-mean Best Rank (PHBR) pipeline methods]. Next, LOHHLA (https://github.com/mskcc/lohhla), originally published in McGranahan and colleagues (17), is used to identify copy number losses of HLA alleles. Copy number and purity data are provided to the program and summary statistics about HLA copy number losses are generated. A given HLA allele was marked as lost if the Pval_unique of its loss was ≤0.05. POLYSOLVER mutation calling [Shukla and colleagues (16)] was used to generate somatic mutation calls of each HLA allele. If an HLA allele was flagged by either of these tools, it was marked as damaged. Alleles were only counted as damaged once even if flagged by both tools. Both programs were provided identical HLA genotypes on a per-patient basis.

We built XGBoost classifiers for three predictive tasks: ICB response prediction from germline, somatic, and combined features, respectively. Models were fit in two stages: feature selection, followed by model training and evaluation. First, we conducted recursive feature elimination (RFE) on an initial array of features using the Cristescu and colleagues (37) cohort, then trained classifiers to predict ICB response using Hugo and colleagues (38). (n = 34), Riaz and colleagues (n = 61; ref. 39), and Snyder (n = 64; ref. 40) melanoma cohorts. The trained model was then evaluated on three test cohorts: Vanallen and colleagues (n = 110; ref. 41), Miao and colleagues (n = 70; ref. 42), and Rizvi and colleagues (n = 34; ref. 43). Biological implication validation was conducted with the Liu and colleagues (n = 122; ref. 44) cohort.

RFE was performed on three feature sets: 229 germline SNPs only, 16 somatic variables only, and both sets combined. The RFE model was trained on Cristescu melanoma (n = 89) and tested on Cristescu HNSCC (n = 107) and Cristescu urothelial (n = 17) samples to ensure this step prioritized broadly useful biological features to use in the model training step. The model used for RFE was an XGBoost random forest classifier (python package version 1.6.2) with 20 total estimators and a maximum depth of 8. We used a nonlinear model for feature selection to allow for feature interactions even during the feature selection stage. All possible feature combinations and total model sizes were tested, and the mean squared error (MSE) of each was recorded. The model with the lowest MSE was selected, and the features included in that model were used for training in stage two. For the 229 germline SNPs, a model with a combination of 54 SNPs yielded the lowest MSE in the RFE cohorts. These 54 SNPs were collapsed into continuous gene-level expression quantitative trait loci (eQTL) scores by measuring the direction of their effect on gene expression in TCGA (see dataset methods for more details) and orienting alleles such that all SNPs affected gene expression in the same direction (Supplementary Fig. S1). This resulted in 23 simplified, gene-level continuous scores reflecting the total magnitude of the expected change in gene expression (Supplementary Fig. S1). For the composite model, RFE was performed on the set of features prioritized by the initial RFE performed for each data type separately.

We trained three different classifiers to predict ICB response, one using only germline features, one using only somatic features, and one on the combined feature set (the composite model). Using features passing RFE analysis, XGBoost random forest classifiers were trained on Riaz and colleagues (39), Hugo and colleagues (38), and Snyder (40). data sets with 1,200 total estimators and a maximum depth of 8. The performance of these models was then evaluated separately on the Vanallen and colleagues (41), Rizvi and colleagues (43), and Miao and colleagues (42) datasets. Aside from feature curation, this process was identical for all models, and a standard random seed was set for all models to ensure reproducibility. For each patient, the XGBoost random forest classifier returns a class prediction probability ranging from 0 to 1, which we refer to as the immune checkpoint (IC) index. For visualization purposes, we used sklearn MinMaxScaler to scale these values from 0 to 10. This process preserves the distribution of scores and therefore does not affect statistical comparisons. For each model and cohort, IC index scores were compared between responders and nonresponders using Mann–Whitney U tests. Comparisons between the effect size of each IC index were made using the Cliff’s Delta value of each model’s effect size. Response to immunotherapy was defined using the RECIST criteria (45). ROC plots were constructed using the scaled continuous IC index results, in which the outcome label was the response phenotype, and the AUC was used to summarize overall performance. Test datasets were then pooled for survival analysis via multivariable Cox proportional hazard analysis, in which the association of IC index with progression-free survival (PFS) was measured alongside covariates of age, sex, and tumor type, using the R packages “survival” and “survminer” (46, 47). Kaplan–Meier curves were constructed using tertile splits of IC index scores and P-values of pairwise comparisons between tertiles were computed with log-rank tests. Finally, positive and negative predictive values (PPV and NPV, respectively) were computed and compared between each model type using the “DTComPair” package (48). State-of-the-art ICB response prediction projects from Litchfield and colleagues, Chowell and colleagues, and Auslander and colleagues (49–51) have demonstrated remarkable accuracy in validation sets when RECIST stable disease (SD) category patients are included as nonresponders or excluded entirely. These SD patients are particularly difficult to classify because of their ambiguous TIME and somatic biomarker landscape, but still benefit from increased overall survival (45) and were counted as responders in predictive modeling tasks in this study.

The composition of immune infiltrates in the TIME was evaluated by digital cytometry via CIBERSORTx using the LM22 signature matrix with batch correction. The T-cell infiltration score was constructed from the CIBERSORTx CD8 T Cells score. The general TIME score used in Kaplan–Meier plotting was calculated as the linear combination of the therapeutic target, T-cell response, and TLS formation. CIBERSORTx T-follicular helper cell estimates were reused for MHC reliance analyses to corroborate the effect of SNPs associated with higher T-follicular helper cell infiltration. The TLS gene expression signature was generated from a set of TLS-related genes reported by Cabrita and colleagues (52) and Sautès–Fridman and colleagues (CCL19, CCL21, CXCL13, CCR7, CXCR5, SELL, LAMP3, CETP, RBP5, AICDA, BCL6, CCR6, and CD79B; 53) using the method put forth in Cabrita and colleagues (52), in which mean gene expression of key genes upregulated in TLS was calculated. CD4⁺ and CD8⁺ T-cell infiltration estimates were calculated using CIBERSORTx, in which the CD4/CD8 ratio was defined using “T-cell CD4 memory activated” + “T-cell follicular helper” infiltration divided by “T-cell CD8” infiltration categories. Only patients in the top two tertiles of CD8⁺ T-cell infiltration were included in direct CD4/CD8 ratio comparison analysis to remove patients with zero or very low levels of immune infiltrates.

Feature importance and interaction within nonlinear models were calculated using the SHapley Additive exPlanations (SHAP) machine learning interpretability suite (https://shap.readthedocs.io/en/latest/). SHAP is a unified approach to explain the output of any machine learning model. It is based on cooperative game theory and the concept of Shapley values. SHAP values assign each feature an importance value for a particular prediction in the context of a specific model. These values allow for nonlinear interactions between features to be accounted for on a per-patient basis, and they also allow us to rank pairwise feature interaction by magnitude. Each model was run through the standard SHAP python pipeline and the feature importances were recorded. For the composite model, feature interaction analysis was performed as well using the shap_interaction_values function.

Originally developed by Marty and colleagues (54, 55), the PHBR score is a measure of how well a given neoantigen is presented by the MHC based on computationally derived binding affinities between all possible peptides harboring the mutation and a patient’s set of HLA alleles. A detailed description can be found in the original publication (54). For each patient, all single nucleotide variant mutations were given an MHC-I PHBR score and an MHC-II PHBR score representing presentation by class-I and class-II, respectively. A neoantigen was considered to be well presented by MHC-I with a PHBR score ≤2 and well presented by MHC-II with a PHBR score ≤10 (56). Class-I HLA alleles were called using POLYSOLVER (v1.0.0; ref. 16) with default parameters, and Class-II HLA alleles were called using HLA-HD (v1.4.0; ref. 57) with default parameters.

Patients were stratified by the ratio of the total number of neoantigens well presented by MHC-II divided by the total number of neoantigens well presented by MHC-I. A patient was only considered for analysis if they had at least three mutations well presented by both MHC-I and MHC-II. Neoantigens that were both well presented by both MHC-I and MHC-II were not considered in this ratio. These ratios were divided into tertiles and defined as follows: the lowest tertile was MHC-I-reliant, the middle tertile was balanced, and the highest tertile was MHC-II reliant. To select for patients with MHC II–based immune responses, MHC II–reliant patients with no evidence of MHC-I damage or loss of heterozygosity were excluded.

Tissue types matching those from our analysis (melanoma, renal cell carcinoma, non–small cell lung carcinoma, head and neck squamous cell carcinoma, and urothelial/bladder cancer) were pulled from TCGA: LUAD, KIRC, SKCM, HNSC, BLCA, KICH, KIRP, and LUSC (see ICB and TCGA datasets). Stage II–IV cancers were analyzed to better match our ICB cohorts. Poorly infiltrated tumors were dropped from the analysis to ensure that cancers analyzed from TCGA were at least somewhat infiltrated by lymphocytes. To achieve this, we calculated the ImmunoScore (58), for all patients, and the bottom tertile (most poorly infiltrated) patients were dropped from the analysis. CD4/CD8 T-cell ratios were calculated in an identical manner as the ICB cohorts. Similarly, MHC reliance groupings were generated identically as in ICB discovery and reliance cohorts.

We estimated tumor intrinsic MHC-II expression by adjusting HLA DRB1 expression in the same cancer types and stages as the above section (see “TCGA immune infiltration analysis”). HLA DRB1 expression levels were corrected for inter-patient variation in immune infiltrates by multiplying by tumor purity fraction and 1 minus the sum of relative infiltration of professional antigen-presenting cells and CD4⁺ T cells, as measured by CIBERSORTx. We adjusted for the following canonical MHC-II expressing cell types in the CIBERSORTx LM22 matrix: B cells, CD4⁺ T cells, macrophages, and dendritic cells.

Five FDA-unapproved IC genes with ongoing clinical trials were investigated for an association with a particular MHC reliance group: LAG3, TIM3, TIGIT, OX40, and IDO1. Univariable analysis revealed significant associations with LAG3 in both discovery and validation ICB cohorts, which were subjected to further multivariable analysis accounting for PDL1 and CTLA4 expression. A median expression cutoff was used to create binary high- and low-expressing groups for each of the checkpoint genes. Age, sex, and tumor type were accounted for during multivariable Cox proportional hazard analysis, as well as prior CTLA4 treatment in the validation cohort, because of a large proportion of patients in Liu and colleagues (52) having received such treatment. Kaplan–Meier curves were generated using these same binary cutoffs and P-values were calculated using the log-rank test.

Statistical software used in this manuscript were R version 4.2.1 and Python version 3.9.2. Unless otherwise indicated all P-value significance thresholds were set at <0.05. Where indicated, P-values were corrected using the Benjamini–Hochberg method.

This study relied entirely on de-identified publicly available datasets and does not necessitate IRB review. The studies relied upon in this manuscript were all conducted in accordance with recognized ethical guidelines and approved by an institutional review board. Code Availability Code to reproduce models, analyses, and figures can be found at the following Github repository: https://github.com/cartercompbio/MHC_reliance. The relevant Code Ocean capsule can be found here: https://codeocean.com/capsule/9714470/tree/v1.

Paired tumor/normal WES data were obtained for eight independent ICB studies encompassing a range of tissue types and treatments across a total of 708 patients (Supplementary Table S1; refs. 37–44). Seven of these were used for machine learning, including feature selection, model training, and independent validation (Fig. 1), and the eighth (Liu and colleagues) was added later to validate the translational potential of biological findings. We first assembled a set of germline and somatic features that could be extracted from WES data and that have previously been reported to predict ICB response (Supplementary Table S2). Germline SNPs associated with the TIME and ICB response from Pagadala and colleagues (Supplementary Table S3; ref. 25) were further harmonized and aggregated at the gene level into numerical scores for their respective gene, here termed eQTL-scores (Supplementary Fig. S1A; Supplementary Table S4). SNPs associated with immune infiltration levels were encoded at the single SNP level instead. Somatic features from several impactful ICB response prediction studies were generated for each cohort, including TMB (59), dN/dS of the immunopeptidome (immunoediting; ref. 15), damage of MHC-I alleles (16, 17), and somatic mutation of genes in the antigen presentation pathway (Supplementary Table S5; ref. 60). Clinical features available for all data sets included patient age and sex (61). To train models to predict ICB response, we used a two-stage machine learning approach entailing feature selection followed by model training (Fig. 1). We first reduced the number of features via RFE using the Cristescu and colleagues (37) cohort before training an XGBoost (62) classifier to predict ICB response as class labels. XGBoost is a tree-based ensemble method that generates a continuous probability score, here scaled to range between 0 and 10. We combined three similar anti-PD1/anti-PDL1/anti-CTLA4 treated melanoma cohorts [Hugo and colleagues (38), Riaz and colleagues (39), and Snyder and colleagues (40)] into a single training set, and evaluated the potential of the classifier to generalize by applying it separately to three heterogeneous independent test cohorts: Van Allen (anti-CTLA4 treated melanoma; ref. 41), Rizvi (anti-PD1 treated NSCLC; ref. 43), and Miao [anti-PD1 or anti-PDL1 (42) some also with anti-CTLA4 treated RCC]. We compared models that relied only on germline features, only on somatic features, or on a combination of both (referred to as the composite model). We termed the scores produced by these models the IC index.

After RFE, we retained 24 germline features to train the germline model (Supplementary Fig. S2A), including 23 germline eQTL-scores representing genes involved in antigen processing/presentation [ERAP2 (63), ERAP1 (64), VAMP8 (65)], immune signaling [FCGR2B (66), PDCD1 (67), CTSS, CTSW)], and DNA replication [DHFR (68), TREX1 (69)] and a SNP associated with infiltration of T-follicular helper cells (T_FHQTL; ref. 70), which was strongly and consistently associated with response across all cohorts (Fig. 2A). RFE for the somatic-only model selected 13 features derived from clinical and tumor genomic data (Supplementary Fig. S2B), including TMB (12), clonality-aware derivatives of TMB such as ITH, and fraction of TMB subclonal (71, 72), as well as DNA based T-cell infiltration estimates (73) and measures of immune evasion (immunoediting, immune escape, MHC-I damage, and antigen presentation pathway damage; refs. 15–17, 74). RFE for the composite model selected 24 features, 18 (75%) of which were germline eQTL-scores and six (25%) of which were somatic features (Supplementary Fig. S2). Considered independently, only a minority of these features showed a significant association with ICB response, and although the direction of effects generally agreed, there was variability across datasets (Fig. 2A). Feature associations with ICB response were more similar across melanoma cohorts than other tumor types (Fig. 2B). Although TMB and clonal TMB features passed RFE in the somatic-only model they were eliminated in the composite model, which instead utilized fraction of TMB subclonal and ITH—features that are anticorrelated and correlated with TMB, respectively (Supplementary Fig. S3; fraction of TMB subclonal: R = −0.22, P = 5.9e⁻⁰⁸; ITH: R = 0.2, P = 1.6e⁻⁰⁶). The composite IC (cIC) index produced by the trained model remained somewhat correlated with TMB (Supplementary Fig. S4A; R = 0.2; P = 0.0035) even though TMB was not directly incorporated as a feature. The somatic IC index had a high correlation with TMB (Supplementary Fig. S4B; R = 0.46; P = 7.5e⁻¹³), whereas the germline IC index was completely uncorrelated with TMB (Supplementary Fig. S4C; R = −0.059; P = 0.39). Finally, purity and ploidy were somewhat correlated in our model (Supplementary Fig. S4D; R = 0.18; P = 4e⁻⁰⁶).

After training XGBoost models on the selected features using the combined training set, we compared the performance of each model on the three independent test sets. Although all three models could distinguish between responders and nonresponders, the cIC index showed the best performance, resulting in the largest mean shift in score distributions between responders and nonresponders (Fig. 3A), the highest Cliff’s Delta between responders and nonresponders (Fig. 3B), and the highest area under the ROC AUC (Fig. 3C). Improvements in ROC AUC from approximately 0.7 to 0.8 were observed in the Van Allen (41) and Rizvi (43) studies, but more modest improvements were observed in Miao (42), possibly due to the vastly different TIME landscape of renal cell carcinomas compared with melanomas (75). PFS of the highest tertile of cIC index scores was significantly higher than the lowest tertile in Kaplan–Meier analysis (Fig. 3D; P < 0.0001) and the cIC index was more predictive of PFS in a Cox proportional hazard analysis using age, sex, and tumor type as covariates (see “Materials and Methods”), with a more extreme HR and more significant P-value relative to germline and somatic-only models (Fig. 3E). Compared with germline and somatic-only models, the cIC index resulted in an increased PPV (Fig. 3F; P = 0.0012; P = 2e⁻⁰⁴), whereas NPV was not significantly different (Fig. 3G). In addition, we found that the germline IC index and somatic IC index scores were completely uncorrelated with each other, suggesting that these sources of data capture orthogonal information (Fig. 3H; R = 0.042; P = 0.54), helping to explain the improved performance of the composite model. The cIC index also outperformed baseline ICB response predictors including TMB, age, gender, and checkpoint expression (Supplementary Fig. S4E–S4I). On a pan-cohort basis, the difference in the cIC index of all responders versus all nonresponders was very significant (Supplementary Fig. S4J; P = 7e⁻⁰⁹).

Next, we compared the cIC index to characteristics of the TIME that can be obtained from RNA sequencing data, which were available for (72/214) 34% of test set patients. Several such measures, including effector CD8⁺ T-cell infiltrates (73, 76), joint B and CD4⁺ T-cell levels potentially indicative of TLS formation (32, 52, 77), and target checkpoint expression (PDL1/CTLA4; refs. 78, 79), have been previously correlated with ICB response. We evaluated CD8⁺ T-cell infiltration levels with CIBERSORTx (80), a digital cytometry tool that estimates immune cell fractions. To model TLS, we used the gene signature developed by Cabrita and colleagues (52) as a proxy for TLS formation. We found that patients split by high versus low cIC index (≥5) generally had similar TIME infiltration levels in all three categories (Supplementary Fig. S5A–S5C). Conversely, the TIME was significantly different between true positives and false positives, in which patients who were predicted to respond (cIC index ≥5) failed to respond and often had an immune-cold TIME, characterized by lower overall levels of immune infiltrates (Fig. 4A; ref. 81). This relationship was strongest in the checkpoint therapy target (CTLA4 for Van Allen and colleagues, PDL1 for Miao and colleagues; P = 0.0081) and TLS formation TIME categories (TLS gene signature P = 0.017), with CD8⁺ T cells showing near significant association (P = 0.055). These results imply that high cIC index patients with favorable germline and somatic biomarkers can nonetheless fail to respond to ICB due to a poorly infiltrated TIME.

We also investigated whether an immune hot TIME could rescue patients with low somatic and germline potential for response. Using a Cox proportional hazard model adjusted for age, sex, and data set, we found that each of the TIME infiltration estimates (checkpoint target: P = 0.0035, CD8⁺ T cells: P = 0.019, TLS formation: P = 0.043) was significantly associated with improved overall survival in high cIC index patients only, whereas low IC index patients failed to significantly benefit from an immune hot TIME (Fig. 4B). These results were mirrored in Kaplan–Meier plots of high and low cIC index patients (Fig. 4C and D) stratified by level of TIME infiltration. High cIC index patients benefit from an above median TIME (P = 0.0097), whereas low cIC index patients do not (P = 0.852). Similarly, a high cIC index with an immune hot TIME had the highest rate of response to ICB (Fig. 4E). These findings are consistent with previous studies indicating that immunogenic tumors respond at greater rates when there is high CD8⁺ T-cell infiltration but that high CD8⁺ T-cell infiltration alone is not sufficient for high rates of ICB response (37). Furthermore, although high cIC index scores yielded the strongest relationship with higher immune infiltration, we found this enhancement was primarily driven by germline factors rather than somatic ones (Supplementary Fig. S5D and S5E). Our analyses suggest that cIC index scores may be useful as general estimates of immunogenicity and could be used as additional indicators of when a patient could benefit from ICB beyond TIME profiling.

In order to better understand how selected germline and somatic features contribute to model performance, we analyzed feature importance using SHAP values (82), a game theory approach to improve the interpretation of the machine learning model. We noted differences in feature rankings particularly for ERAP1, MHC-I damage, and immunoediting, between XGBoost and linear models suggesting the presence of interactivity effects (Supplementary Fig. S6). Thus, we evaluated both individual feature contributions and pairwise interactions between features. SHAP analysis revealed several key feature interactions (Fig. 5A), the strongest of which was between somatic MHC-I damage, i.e., the cumulative MHC-I damage from somatic mutation and loss of heterozygosity (see “Materials and Methods”), and T_FHQTL. We further examined this interaction in terms of ICB response rates between categories (Fig. 5B) and observed higher rates of response when the T_FHQTL was present (P = 1.7e⁻⁰⁵ T_FHQTL vs. class-I MHC damage; P = 0.0063 T_FHQTL vs. neither), even when the potentially negative effect of MHC-I damage was present (P = 1.0; T_FHQTL vs. both). Because rates of ICB response were unaffected by MHC-I damage in patients carrying the T_FHQTL (Fig. 5B), we hypothesized that this SNP may promote immune responses upon ICB treatment that do not rely on MHC-I-based antigen presentation. Instead, patients with the T_FHQTL could be predisposed toward an MHC-II–driven mechanism of response. As T_FH cells primarily assist B cells in producing antibodies through MHC-II interactions (83), this may suggest a role for humoral immunity that is independent of cellular immune responses mediated via MHC-I and cytotoxic T cells. NK cells are also known to modulate adaptive immune responses via CD27 in MHC-I–deficient tumors (84). Although the role of follicular T helper cells (Tfh) in this context is not well understood, one study in mice with MHC-I–deficient tumors found that NK cell–CD4⁺ T-cell interplay led to tumor rejection without any CD8⁺ T-cell activity (85).

To further investigate this idea, we grouped tumors in the dataset according to whether somatic mutations were more prevalently presented by MHC-I or MHC-II molecules, suggesting the potential for the reliance of immune responses on particular MHC pathways of neoantigen presentation. First, we calculated PHBR scores (54, 55) for each nonsynonymous mutation in all patients. PHBR scores are mutation-centric scores that seek to summarize whether any peptides overlapping the mutated site will be presented by any of an individual’s HLA alleles. Patients with at least three mutations passing PHBR thresholds for both MHC-I and MHC-II were then split into groups termed MHC I reliant, MHC II reliant, or balanced based on the ratio of these class-specific neoantigens (Fig. 5C; Supplementary Fig. S7A; Supplementary Table S6), with reliant referring to an immune response potentially dependent on MHC-I versus MHC-II presented neoantigens. Among MHC I–reliant patients, we noted a significantly higher level of MHC-I damage in nonresponders versus responders (P = 0.0092; Fig. 5D) reflecting the notion that an MHC I–reliant response depends on the integrity of the MHC-I and associated antigen presentation pathway. Although balanced patients demonstrated an intermediate disparity in MHC-I damage between nonresponders versus responders (P = 0.02), this was not the case in MHC II–reliant patients (P = 0.74). Overall ICB response rates between these two groups were not significantly different (Supplementary Fig. S7B).

Next, we sought to understand how MHC reliance could modify the potential to benefit from the T_FHQTL. We reasoned that the most extreme cases of MHC-II reliance would be those that also had defects in the MHC-I antigen presentation pathway. The distribution of defects to the MHC-I antigen presentation pathway was statistically similar between each of the MHC reliance groups (Supplementary Fig. S7C and S7D), although MHC II–reliant patients with defects to the MHC-I antigen presentation pathway showed significantly less immunoediting than those with defects (Supplementary Fig. S7E), suggesting that these defects may limit a patient’s ability to mount an MHC-I driven immune response. MHC II–reliant patients with defects to the MHC-I antigen presentation pathway comprised 83% of MHC II–reliant tumors (154/171), so we focused further analyses on this subpopulation. We found a significant difference in the frequency of the T_FHQTL between responders versus nonresponders in the MHC I–reliant and balanced categories (P = 0.0042; P = 0.003; Fig. 5E) but not in the solely MHC II–reliant category (P = 0.12). This is somewhat mirrored in the subset of patients with tumor immune infiltration estimates available, in which T_FH cell estimates were higher in MHC I–reliant responders versus nonresponders (P = 0.03; Supplementary Fig. S7F) but not in the balanced or MHC II–reliant responders versus nonresponders (P = 0.48, P = 0.5). It is possible that MHC I–reliant responders benefit from an increased infiltration by T_FH cells, TLS formation, and associated helper effects that are important to maintain the function and precursor frequency of CD8⁺ T cells (86–91). Indeed, TLSs have been shown to enhance ICB response in melanoma (52, 77). Conversely, MHC II–reliant patients may receive less benefit from additional T_FH-cell infiltration because their neoantigen landscape is already predisposed toward the formation of TLSs. Indeed, we found that MHC I–reliant responders had higher TLS gene signature expression than nonresponders (P = 0.036; Fig. 5F), yet this difference was not significant in MHC II–reliant patients (P = 0.12; Fig. 5F). MHC II–reliant patients in general had a higher level of TLS gene signature expression than MHC I–reliant patients (P = 0.0088; Fig. 5F), which is consistent with the fact that TLS formation is more closely associated with the MHC-II/CD4⁺ T-cell axis (53, 92, 93). These initial observations point to the possibility that mechanistically divergent immune responses yield ICB responses based on how effectively neoantigens engage each MHC pathway.

We next sought to understand the clinical implications of differential MHC reliance. To validate our findings, we performed identical analyses on an additional independent ICB-treated cohort (n = 77) with paired transcriptomic data [Liu and colleagues (44)] and compared the results with those from our original set of seven cohorts referred to as the discovery set. We first investigated the effects of MHC reliance grouping on the composition of the TIME. We observed that CD4⁺/CD8⁺ T-cell ratios mirrored MHC reliance in responders, with higher ratios being observed in MHC II–reliant tumors (Fig. 6A; P = 0.0057 discovery; P = 0.025 validation). However, no such difference was found in nonresponders. We applied an identical methodology to immune-infiltrated (58) ICB-naive, tissue-matched cancer samples from TCGA and found a protective effect of the CD4⁺/CD8⁺ T-cell ratio in TCGA MHC II–reliant patients (Supplementary Fig. S8A; HR = −0.76; P = 0.0069) but a significantly adverse effect of that same ratio in TCGA MHC I–reliant patients (Supplementary Fig. S8B; HR = 0.59; P = 0.0352). Additionally, we found that an estimate of tumor intrinsic MHC-II expression was protective in MHC II–reliant patients only (Supplementary Fig. S8C; HR = −0.64; P = 0.0094; Supplementary Fig. S8D; HR = 0.14; P = 0.65). These data support the idea that there is a benefit to having some level of concordance between CD4⁺/CD8⁺ T-cell infiltration, MHC-II expression, and MHC-II/MHC-I neoantigen ratios. To investigate differences in response dynamics between CD4⁺ and CD8⁺ T-cell mediated responses, we compared the survival of responders MHC II–versus MHC I–reliant groups. Despite nonsignificant differences in response rates, MHC II–reliant responders had a significantly longer overall survival in both discovery and validation cohorts (Fig. 6B and C, discovery P = 0.0073; validation P = 0.0398), consistent with reports that CD4⁺ T-cell–based immune responses are tumor autonomous and therefore more difficult to evade in the long term (55, 94, 95). This observation was not solely reliant on one cohort or cancer type, as MHC reliance groupings were found to be balanced across all cohorts except Van Allen (Supplementary Table S7), and our findings were unchanged upon removal of the Van Allen cohort from this analysis (Supplementary Fig. S9).

Finally, we wanted to know if differences in MHC reliance could translate to differences in pathways of immune evasion. ICs are commonly overexpressed to suppress an active immune response. Currently, of the many checkpoints identified in the tumor microenvironment only PDL1 positivity in tumor sections is approved as a biomarker of ICB response, albeit its predictive value is modest (78). To investigate whether differences might exist about which checkpoints correlate with a beneficial antitumor immune response under different MHC reliance conditions, we evaluated the relationship between the expression of individual checkpoint genes and PFS post-ICB treatment by univariable Cox PH analysis. We focused on checkpoint genes with antibody inhibitors undergoing clinical trials (PDL1, CTLA4, LAG3, TIGIT, TIM3, IDO1, and OX40; ref. 96). When split by MHC reliance grouping, higher LAG3 expression was associated with benefit from ICB in the MHC II–reliant group (Supplementary Fig. S10; Supplementary Table S8). To adjust for potentially confounding effects of the correlated expression of canonical (97, 98) immune checkpoint genes, we performed a multivariable analysis centered on LAG3, PDL1, and CTLA4. We found that high PDL1 expression was generally associated with longer survival post-ICB treatment in MHC I–reliant patients (Fig. 6D, discovery P = 0.026; validation P = 0.062), CTLA4 expression with longer survival in balanced patients (Fig. 6D, discovery P = 0.054; validation P = 0.006), and LAG3 with longer survival in MHC II–reliant patients (Fig. 6D, discovery P = 0.014; validation P = 0.002). There was no association of checkpoint gene expression with the MHC reliance category (Supplementary Fig. S11A and S11B), and we observed some autocorrelation between LAG3, PDL1, and CTLA4 (Supplementary Fig. S11C and S11D). Among MHC II–reliant patients, higher expression of LAG3 was associated with significantly longer overall survival in both discovery and validation cohorts (Fig. 6E and F, discovery P = 0.0018; validation P = 0.0345). LAG3 is thought to play a prominent role in CD4⁺ T-cell regulation and may be a primary marker of activation (99, 100). Our results may, therefore, reflect a key role for LAG3 as a mediator of CD4⁺ T-cell–based response to ICB therapy.

ICB has emerged as a potent anticancer therapy; however, the fraction of patients who benefit from treatment remains disappointingly low. To improve the success of ICB, it is of the utmost importance to understand which factors govern the potential to respond via the immune system. Here, we used a machine learning framework to study somatic and germline biomarkers of response to ICB in human cohorts. We were able to extract both feature types from paired tumor-normal WES data across eight ICB-treated human studies. Germline immune eQTL biomarkers, whereas relatively new, show promise to capture complementary information from somatic features, and XGBoost models trained to predict a cIC index using both feature types performed better at predicting ICB response across different tumor types. When we interrogated patients with additional available RNAseq data, we found that the survival benefit of an immune hot microenvironment was contingent upon having a high cIC index score, that there was no response in patients with a low cIC index score, and that this was driven by germline features. This supports the notion that heritable differences in immune-cell function determine the effectiveness of an immune response once immune cells have reached the tumor. Furthermore, patients with a high cIC index score who failed to respond often had a “cold” TIME. This suggests that transcriptomic profiling might be useful as a supplemental prognostic tool of ICB response in high cIC index patients and that the cIC index score serves as a general proxy for clinical response to the immune invigorating effect of ICB.

To gain further insight about how various biomarkers relate to ICB response potential, we used state-of-the-art techniques for interpreting machine learning models and studied important features and feature interactions that drove model predictions. The strongest interaction involved an interplay between an SNP associated with increased infiltration of T-follicular helper cells (T_FHQTL) and MHC-I damage. Specifically, we observed a beneficial effect of the T_FHQTL on rates of response, independent of the deleterious effect of MHC-I damage. T_FH cells are the specialized subset of CD4⁺ T cells that help B cells produce antibodies in germinal centers (53). T_FH cells are normally located in secondary lymphoid organs at a close distance from B cells (83). However, there is increasing evidence that T_FH cells are part of TLS, intra-tumor organized clusters of immune cells including B and T cells and dendritic cells mimicking germinal centers in secondary lymphoid organs (53, 101). TLS are an increasingly common finding in cancer, and are linked with better prognosis (102, 103); increased infiltration by T_FH cells and TLS formation are a source of helper factors beneficial to both CD8⁺ and CD4⁺ T cells. Indeed, the number of TLS distinguishes ICB responders from nonresponders (32, 77).

MHC-I damage on cancer cells inherently hampers the cytotoxic function of CD8⁺ T cells, yielding low response rates. We found that response rates were rescued when patients had both the T_FHQTL and MHC-I damage, suggesting that rescue mechanisms of ICB response may be shifted toward MHC-II mediated immunity (MHC-II reliance). Using individual-level information about the ratio of neoantigens with binding affinity for MHC-I and MHC-II, we were able to allocate patients to either an MHC I– or MHC II–reliant group. That these groupings may initiate and sustain differential immune mechanisms in response to ICB is strengthened by the observation that MHC-II reliance promotes higher infiltration of CD4⁺ T cells and more durable clinical responses to ICB, potentially reflecting a direct effect on long-term memory CD4⁺ T-cell responses. In contrast, MHC I–reliant responses, which are centered on CD8⁺ T cells, are possibly more transient in the absence of CD4⁺ T-cell help (87).

When we examined the association of pretreatment checkpoint gene expression levels with ICB response, which was predominantly anti-PD1/anti-PDL1 treatment in the cohorts studied, we found that PDL1 expression was associated with better ICB response in MHC I–reliant patients but not in MHC II–reliant patients, whereas the reverse was true for LAG3. In patients in whom immune evasion is mediated by overexpression of PD1/PDL1, anti-PD1/anti-PDL1 therapies can be remarkably effective (104). In contrast, LAG3 has MHC-II as its major ligand (100), and it is widely regarded as a negative regulator of CD4⁺ T-cell activation (105). Higher expression of LAG3 could therefore indicate an effective ongoing MHC II–reliant antitumor response pre-ICB treatment. In our analysis, LAG3⁺ patients had better survival in the MHC II–reliant group, suggesting that MHC-II-driven immunity can support an effective response to anti-PD1/anti-PDL1 and that this could potentially be further amplified by an anti-LAG3 therapy. However, the lack of association of PDL1 expression with response in the MHC II–reliant group seems to suggest a mechanism independent of alleviating PDL1-based repression of CD8⁺ T cells. A similar phenomenon has been observed in microsatellite instability colorectal cancers with B2M loss that paradoxically remain among the best responders to anti-PD1/anti-PDL1 therapy (106). It is intriguing to think that anti-PD1/anti-PDL1 can be beneficial even if PDL1 is not highly expressed or the MHC-I antigen presentation machinery is not functional. Recent data show that LAG3 also associates with the T-cell receptor (TCR)–CD3 complex in both CD4⁺ and CD8⁺ T cells in the absence of binding to MHC-II, causing the dissociation of the tyrosine kinase Lck from the CD4 or CD8 co-receptors and loss of co-receptor–TCR signaling during T-cell activation (107). Our finding that LAG3 facilitates the CD4⁺ T-cell responses during ICB treatment could be explained by the fact that both LAG3 and ICB target the proximal signaling of the TCR (108), even though the reasons this creates an advantage in MHC II–reliant patients remains unclear. Perhaps this reflects the fact that the adult peripheral repertoire is richer in CD4⁺ than in CD8⁺ T cells. This bias may also explain the observation that patients with cancer vaccinated with neoantigens have a propensity to generate CD4⁺ T-cell responses (109).

The other implication is that the utility of each of these checkpoint genes as biomarkers of ICB response may be highly context-dependent. PDL1 expression was not associated with ICB response in MHC II–reliant patients responding via a CD4⁺ T-cell axis of adaptive immunity. This could explain in part why PDL1 positivity is a surprisingly poor general predictor of response rates (110). Future efforts to refine biomarkers of ICB response could attempt to leverage widely available germline information as well as understand the context of a patient’s MHC reliance status.

Our study has some limitations. First, there are a limited number of publicly available ICB-treated cohorts with DNA sequencing data available and even fewer with RNA sequencing data available. Second, larger feature selection and training cohorts could further improve model performance. Future studies could incorporate additional biomarkers, for example, genotypes associated with adverse immune events, such as rs16906115 affecting IL7 (111), that could lead to early stopping of therapy, or copy number alterations affecting key immune loci (112, 113). Third, we limited our features to those extractable from paired tumor-normal WES as tumor DNA to mirror what is more commonly available in real-world settings. Additionally, although the germline-derived features in the composite model are straightforward to compute once the bioinformatic infrastructure is in place, the variety and complexity of the somatic features may be more challenging to implement in the clinic. MHC reliance groupings were based solely on single nucleotide variants. Future versions of our PHBR pipeline will include support for frameshift and stop-loss variants, which may be more impactful in an immunogenicity context. Another limitation is that most ICB response classification approaches eliminate difficult-to-classify SD patients from their studies—despite the fact that these patients benefit from increased survival from ICB treatment. We chose to include these patients as responders to maximize potential clinical benefit, at the cost of increasing the complexity of our classification task. Finally, although our classifier—which was trained on patients with melanoma—showed some ability to generalize to other tumor types, especially non–small cell lung cancer, it may ultimately be essential to train and study tumor-type specific models.

Investigation of the factors that determine ICB response in patients with cancer is providing key insights into mechanisms that drive superior response. This study provides further evidence that CD4⁺ T-cell responses engaged by MHC-II antigen presentation are a critical component of superior immune responses and points to an alignment of checkpoint-based evasion with the immune cell types dominating the response. This sets the stage for future strategies to optimize the selection of checkpoint therapies from characteristics of the patient’s tumor and immune system.

M.S. Pagadala reports a patent for SD2022-086 issued. A. Castro reports personal fees from Tempus Labs outside the submitted work. J. Kong reports grants from T32CA121938 outside the submitted work. H. Carter reports grants from The Mark Foundation for Cancer Research, the National Cancer Institute, and the National Institute of General Medical Sciences during the conduct of the study; in addition, H. Carter has a patent for immunotherapy response prediction tools pending. No disclosures were reported by the other authors.

T.J. Sears: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. M.S. Pagadala: Conceptualization, data curation, software, methodology, writing–review and editing. A. Castro: Data curation, software, methodology. K.-h. Lee: Conceptualization, visualization, writing–original draft, writing–review and editing. J. Kong: Writing–original draft, writing–review and editing. K. Tanaka: Formal analysis, visualization. S.M. Lippman: Supervision, writing–original draft, writing–review and editing. M. Zanetti: Conceptualization, formal analysis, supervision, methodology, writing–original draft, writing–review and editing. H. Carter: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.

This work was funded by Mark Foundation Emerging Leader Award #18-022-ELA, NCI grant R01CA269919, and support from NCI grant U24CA248138 to H. Carter. Computational resources were supported by infrastructure grant 2P41GM103504-11. The results shown here are in part based on data generated by the TCGA Research Network. ICB datasets: For Rizvi and colleagues non–small cell lung cancer immunotherapy analysis, we used dbGaP data from accession phs000980.v1.p1. We thank the members of the Thoracic Oncology Service and the Chan and Wolchok labs at Memorial Sloan Kettering Cancer Center (MSKCC) for helpful discussions. We thank the Immune Monitoring Core at MSKCC, including L. Caro, R. Ramsawak, and Z. Mu, for exceptional support with processing and banking peripheral blood lymphocytes. We thank P. Worrell and E. Brzostowski for their help in identifying tumor specimens for analysis. We thank A. Viale for superb technical assistance. We thank D. Philips, M. van Buuren, and M. Toebes for help performing the combinatorial coding screens. This work was supported by the Geoffrey Beene Cancer Research Center (MDH, NAR, TAC, JDW, and AS), the Society for Memorial Sloan Kettering Cancer Center (MDH), Lung Cancer Research Foundation (WL), Frederick Adler Chair Fund (TAC), The One Ball Matt Memorial Golf Tournament (EBG), Queen Wilhelmina Cancer Research Award (TNS), The STARR Foundation (TAC and JDW), the Ludwig Trust (JDW), and a Stand Up To Cancer-Cancer Research Institute Cancer Immunology Translational Cancer Research Grant (JDW, TNS, and TAC). Stand Up To Cancer is a program of the Entertainment Industry Foundation administered by the American Association for Cancer Research. For Snyder and colleagues melanoma immunotherapy analysis, we used dbGaP data from accession phs001041.v1.p1. We thank Martin Miller at MSKCC for his assistance with the NetMHC server, Agnes Viale and Kety Huberman at the MSKCC Genomics Core, Annamalai Selvakumar and Alice Yeh at the MSKCC HLA typing laboratory for their technical assistance, and John Khoury for assistance in chart review. For Miao and colleagues renal cell carcinoma immunotherapy analysis, we used dbGap data from accession phs001493.v2.p1. This study was supported by an American Association for Cancer Research KureIt grant. Hugo and colleagues melanoma samples were acquired from SRA using accession numbers SRP067938 and SRP090294. Riaz and colleagues melanoma samples were acquired from SRA using accession number SRP095809. For Van Allen and colleagues melanoma sample, data was acquired from dbgap accession phs000452.v2.p1. For Liu and colleagues melanoma validation cohort, data was acquired from dbgap accession phs000452.v3.p1 and supported by the National Human Genome Research Institute (NHGRI) Large Scale Sequencing Program, Grant U54 HG003067 to the Broad Institute (PI, Lander).