Genomic correlates of response to immune checkpoint blockade in microsatellite-stable solid tumors.

Tumor mutational burden correlates with response to immune checkpoint blockade in multiple solid tumors, although in microsatellite-stable tumors this association is of uncertain clinical utility. Here we uniformly analyzed whole-exome sequencing (WES) of 249 tumors and matched normal tissue from patients with clinically annotated outcomes to immune checkpoint therapy, including radiographic response, across multiple cancer types to examine additional tumor genomic features that contribute to selective response. Our analyses identified genomic correlates of response beyond mutational burden, including somatic events in individual driver genes, certain global mutational signatures, and specific HLA-restricted neoantigens. However, these features were often interrelated, highlighting the complexity of identifying genetic driver events that generate an immunoresponsive tumor environment. This study lays a path forward in analyzing large clinical cohorts in an integrated and multifaceted manner to enhance the ability to discover clinically meaningful predictive features of response to immune checkpoint blockade.

Introduction

Immune checkpoint inhibitors, including monoclonal antibodies targeting the immune inhibitory proteins programmed cell death protein 1 (PD-1) and cytotoxic T lymphocyte associated antigen 4 (CTLA-4), significantly extend patient survival across many cancer types[1]. However, pre-treatment predictive biomarkers of response remain unclear. Immunohistochemical (IHC) assessment of tumor immune activity, including staining for the ligand to PD-1 (PD-L1) and the presence of CD8+ T cells, may help predict response to anti-PD-1 therapies[1,2,3], but these histological snapshots may not fully represent a dynamic tumor-immune microenvironment[4]. Further, they are limited by variable sensitivity and specificity by tumor type[5] and minimal predictive value for anti-CTLA-4 therapies.

Molecular analysis of patient tumors has revealed that a high burden of tumor-specific mutant peptides (neoantigens) generated from somatic nonsynonymous coding mutations may increase tumor immunogenicity and the likelihood of patient benefit from immune checkpoint therapy in both microsatellite stable and unstable tumor types[6,7,8,9,10]. However, not all studies of immunotherapy-treated cohorts have found a strong relationship between mutational burden and response[11,12]. Even in studies that identified a statistically significant association, the distributions of mutational burden between responders and non-responders overlapped substantially[10], and the opportunities for therapeutic targeting of this feature to enhance response remain uncertain. Even the most optimal cutoff yields weak sensitivity (74%) and specificity (59.3%) to discriminate potential clinical benefit[13], limiting the utility of tumor mutational burden as a clinical biomarker for individual patients.

Given the major limitations of mutational burden as a predictive biomarker, additional genomic studies have suggested that clonal mutations and neoantigens[14], mutations and copy number alterations affecting particular genes and signaling pathways[15,16,17,18,19], and overall tumor aneuploidy[12,20] may have additional predictive value for response or resistance to immune checkpoint therapies. However, existing studies have focused on cancer types individually without identifying pan-cancer relationships; have been of limited sample size; and have had inconsistent computational methods. These limitations have restricted the power of these studies to identify meaningful associations and have likely contributed to difficulties in validating prior findings in independent prospective patient cohorts[21].

We hypothesized that an expanded and uniformly analyzed cohort of clinically-annotated patient samples would provide greater power to detect significant associations between pre-treatment tumor characteristics and response to immune checkpoint therapies. Thus, we gathered raw tumor and germline pre-treatment whole exome sequencing data from tumors from immune checkpoint-treated patients from seven published studies (N=171)[7,8,10,22,23,24,25] and combined it with data from 78 newly-sequenced pre-treatment tumors. By harmonizing clinical annotations and whole exome analyses across 249 samples and multiple cancer types, we aimed to 1) assess the generalizability of prior hypotheses regarding immune checkpoint blockade response to other histological or drug settings; 2) apply new computational techniques for inference of tumor biology and immunogenicity; and 3) determine whether our cohort was of sufficient size to reveal statistically robust associations about specific genetic mediators of selective response.

Results

Consolidation of a clinically annotated cohort of whole exome sequencing of tumors from patients treated with immune checkpoint blockade

Whole exome sequencing (WES) from clinically annotated tumor samples with matched germline blood or adjacent normal tissue was available for 314 patients (Supplementary Table 1, Fig. 1a). Standard quality control measures were taken to ensure adequate power to detect tumor-specific mutations[26] (Fig. 1a, Supplementary Fig. 1a–c). Our final analysis cohort included 249 patient tumors across six cancer types: melanoma (N=151), non-small cell lung cancer (NSCLC) (N=57), bladder cancer (N=27), head and neck squamous cell carcinoma (HNSCC) (N=12), anal cancer (N=1), and sarcoma (N=1) (Fig. 1a). These patients were treated with anti-PD-1 (N=74), anti-PD-L1 (N=20), anti-CTLA-4 (N=145), or a combination of anti-CTLA-4 and anti-PD-1/L1 therapies (N=10). A small minority of patients received anti-PD-1, anti-PD-L1, or anti-CTLA-4 therapy in combination with another immunotherapy, targeted therapy, or cytotoxic chemotherapy (N=7) (Supplementary Table 2). Across these 249 samples, average mean target sequencing coverage was 150-fold for tumor tissue and 119-fold for matched germline tissue. Mean estimated tumor purity was 58% (range 10%–97%) (Supplementary Table 1).

Fig. 1

Clinical cohort consolidation, response stratification, and mutational load investigation

(a) Data quality control for 249 samples included in final analysis. (b) Comparison of three published response metrics for determining OR vs. NR. (c) Comparison of tumor mutational burden between CR/PR vs. PD (For ‘All mutations’, p=0.0005 for CR/PR vs PD, p=0.0054 for CR/PR vs SD, and p=0.434 for SD vs PD; for ‘Nonsynonymous mutations’, p=0.0003 for CR/PR vs PD, p=0.0063 for CR/PR vs SD, and p=0.3769 for SD vs PD; for ‘Clonal nonsynonymous mutations’, p=0.00005 for CR/PR vs PD, p=0.011 for CR/PR vs SD, and p=0.151 for SD vs PD). Outlying points from patients with mutations/Mb > 101 are not shown (2 CR/PR, 1 SD, 3 PD). (d) Intratumoral heterogeneity across response groups (n=249 biologically independent samples, p=0.001 for CR/PR vs. PD, p=0.5122 for CR/PR vs SD). (e) Clinical response to immune checkpoint therapy broken down by intratumoral heterogeneity. For (c, d), p-values calculated by two-sided Mann-Whitney U *p<0.05, **p<0.005, ns = not significant. Boxplots show the median, first and third quartiles, whiskers extend to 1.5 × the interquartile range, and outlying points are plotted individually.

In selecting a framework to assess clinical response in this study, we encountered multiple patient stratification methods from prior studies of response predictors to immune checkpoint therapy[10,11,12,27]. These varied mostly in their treatment of patients with stable disease (SD) by RECIST 1.1[27], who have minimal to no change in tumor burden following therapy. Applying three existing response definitions to our cohort, we observed substantial differences in patient classification into responder and nonresponder groups (Fig. 1b, Supplementary Table 3). Given the evolving viewpoints on classifying response to immune checkpoint blockade[28], we adopted a conservative method of defining objective response (OR) as complete response (CR) or partial response (PR) by RECIST (major decrease in tumor burden following treatment) and no response (NR) as progressive disease (PD) by RECIST (major increase in tumor burden following treatment) for the main analyses. Patients with SD were considered separately, and analyses using two other response definitions[10,12] that stratify SD patients into OR vs. NR by duration of overall survival (OS) or progression-free survival (PFS) are available in the Supplementary Materials.

Mutational burden and response to immune checkpoint therapy

In examining whole exome genetic features in this cohort, we began with tumor mutational burden, as this has been the most widely reproduced association with response to immune checkpoint therapy. We found that in this combined cohort, patients with CR/PR had significantly higher tumor mutational burdens compared to patients with PD (p<0.05 for all; Mann-Whitney U) (Fig. 1c). This finding persisted within cancer types (Supplementary Fig. 2), and was particularly prominent in patients treated with anti-PD-1/PD-L1 inhibitors (Supplementary Fig. 3a). Patients with SD tended to have mutational burdens intermediate between those with PD or CR/PR, with higher mutational loads in patients with SD with long OS compared to short OS (Supplementary Fig. 3b).

While these findings are consistent with the growing body of literature supporting the association between mutational burden and immune checkpoint therapy response, we noted that the ranges of mutational burdens between response groups overlapped considerably (Fig. 1c), and we found that tumor mutational burden had poor predictive power to differentiate CR/PR vs. PD as a single variable in this cohort (AUC = 0.66) (Supplementary Fig. 4). To build upon the utility of tumor mutational burden as a predictive variable for response to immune checkpoint therapy, past studies have determined that clonal (found in every cancer cell) mutations rather than subclonal (found in a subset of cancer cells) mutations are more strongly associated with response in lung adenocarcinomas and some melanomas, potentially due to stronger T cell responses to neoantigens generated from clonal vs. subclonal mutations[14].

We queried our cohort for this association by using ABSOLUTE to infer mutational clonality[26], and demonstrated that clonal nonsynonymous mutational burden strongly predicted CR/PR vs. PD across cancer types and response categorizations (Fig. 1c–d, Supplementary Fig. 2). Patients with a large proportion of subclonal mutations (>50%) – which we term high intratumoral heterogeneity – were substantially more likely to have PD than CR/PR across all tumors described herein (Fig. 1e) (p=0.0014; Fisher’s exact). Thus, while mutational burden begins to explain the variance in patient response to immune checkpoint therapy, intratumoral heterogeneity contributes additional biological insight.

Mutations in specific genes associated with response or resistance to immune checkpoint therapy

Given the complexity we observed in the association between mutational load and response to immune checkpoint therapy, we pursued additional analyses of exome-level features that could provide further nuance to this association. We next investigated whether somatic mutations in specific genes were associated with response to immune checkpoint therapy, hypothesizing that this analysis would identify genes driving biological processes generating large numbers of mutations or those creating an immunoresponsive phenotype independent of mutational burden.

Given that response rates were similar across all cancer types (Supplementary Table 3), we first compared nonsynonymous mutations in CR/PR vs. PD across all genes; however, these associations were too weak to pass multiple hypothesis test correction (Supplementary Fig. 5). Thus, we limited our analysis to known hotspot mutations in cancer driver genes and to loss-of-function alterations in known tumor suppressors[26] (Supplementary Table 4), as these events are more likely to have a significant impact on tumor biology. Clonal driver alterations in PIK3CA, KRAS, and PBRM1 were enriched in CR/PR, while clonal driver mutations in EGFR were enriched in PD (p<0.05 for all, Benjamini-Hochberg false discovery rate (FDR) q=0.18 for KRAS and PIK3CA) (Fig. 2a). After correcting for tumor mutational burden, KRAS and PIK3CA remained associated with CR/PR (p<0.05; logistic regression) (Fig. 2b), although these observations did not pass FDR correction and may be confounded by additional aspects of a tumor’s genetic profile.

Fig. 2

Mutations in specific genes associated with response to immune checkpoint therapy

(a) Response-associated mutations in CR/PR vs. PD (Two-failed Fisher’s exact test, n = 70 biologically independent samples with CR/PR, n = 123 with PD). Dashed red line indicates Fisher’s exact p=0.05, and dashed dark red line indicates false discovery rate (FDR) q=0.05. (b) Response-associated mutations corrected for mutational burden. Sample size and dashed lines for p- and q-value cutoffs same as in (a). (c–d) Tile plot showing known hotspot and non-hotspot mutations in (c) PIK3CA and (d) KRAS by response group (bottom), with cancer type indicated by capital letters. Top four rows of (c) represent mutations arising in APOBEC-associated mutational contexts (See Supplementary Fig. 6). (e) Simulated statistical power calculation for detection of response-associated genes. Significance of association between response and presence of mutation in gene (Two-tailed Fisher’s exact test) is shown on the y-axis, for varying samples sizes (x-axis). Colors represent frequency of mutations and specificity of alteration to OR vs. NR. Dashed horizontal line represents Bonferroni-correct p=0.05 after correcting for multiple comparisons over the 116 cancer driver genes assessed in this study. Simulated cohort contains 40% CR/PR and 60% PD patients.

Driver mutations within a given gene may occur in different functional domains and have different phenotypic effects, often dependent on the cancer context. Thus, we next examined these cohort-wide gene associations for trends within cancer types. Of the patients with clonal hotspot mutations in PIK3CA, those with CR/PR had melanoma, HNSCC, anal cancer, or bladder cancer, while the majority of those with SD or PD had lung cancer (Fig. 2c, Supplementary Table 5, Supplementary Fig. 6). Hotspot mutations in KRAS occurred predominantly in CR/PR across multiple cancer types (Fig. 2d). Genetic events in PIK3CA and KRAS were too infrequent to fully clarify their relationships to response in this study, but these results demonstrate that single-gene associations with response to immune checkpoint therapy can provide additional information beyond mutational burden, in a manner that may be dependent on or independent of cancer type.

Despite combining data from multiple studies and cancer types, our analyses were still statistically underpowered to detect important relationships; thus, we sought to estimate the sample sizes needed for discovery of single-gene correlates of response with appropriate correction for multiple hypothesis testing[26]. We modeled statistical significance values for common or rare variants associated with CR/PR vs. PD at various sample sizes (Fig. 2e). In the best-case scenario, where a variant is both relatively common (~10% prevalence) and specific to responders, sample sizes of around 300 would be adequate to detect significant associations. Meanwhile, detection of rare response-associated variants (~1% frequency), even if highly specific, would necessitate sample sizes in the thousands. Thus, continuing to combine independent clinically-annotated cohorts will aid in increasing statistical power to detect common response-associated variants, but for detection of rarer events, applying insights from experimental studies for hypothesis-driven validation will be a crucial adjunct.

Integrated analysis of response- and resistance-associated mutations with mutational signatures

As a further step towards understanding the mutational processes that generate an immunoresponsive tumor environment, we next investigated whole exome signatures of mutagenic biological processes. We used a previously described non-negative matrix factorization technique to identify known mutagenic processes in lung cancer, melanoma, bladder cancer, and HNSCC[26,29,30]. Using this technique, the somatic mutations within a tumor are probabilistically assigned to underlying mutational signatures, which are patterns of somatic mutations thought to arise from carcinogenic processes (e.g. a predominance of C-to-A transversions in tobacco-smoking-associated cancers or C-to-T transitions in ultraviolet-light-associated tumors) (Supplementary Table 6, Supplementary Table 7). As has been seen previously in analyses of non-clinically-annotated NSCLC exomes[31], pre-treatment NSCLC tumors in this cohort with a high proportion of smoking-associated mutations tended to have low intratumoral heterogeneity, high mutational burden, and frequent KRAS mutations (Fig. 3a–c). The proportion of smoking-associated mutations was also higher in CR/PR patients compared with PD (Supplementary Fig. 7a). Meanwhile, tumors with EGFR hotspot mutations trended towards enrichment in subclonal mutations (p=0.035), had low mutational burdens, were over-represented in never-smokers (Two-tailed Fisher’s exact test,p = 0.00017), and had generally poor responses to immune checkpoint therapy (Fig. 3a, d), which is also consistent with previously published results[32,33]. Thus, the association between KRAS and EGFR mutation status and response to immunotherapy (Fig. 2a–b) may be related not only to mutational burden but also to mutational signature and mutational clonality. After controlling for smoking history, mutational burden remained a significant predictor of response (Fig. 3e, Supplementary Fig. 7b). Concurrent consideration of mutational signatures, clonal architecture, and hotspot mutations in NSCLC enhances understanding of somatic mediators of immunotherapy response and resistance.

Fig. 3

Integrated analysis of EGFR mutational status, intratumoral heterogeneity, and mutational signatures in lung cancer

(a) Stacked plot showing mutational burden (histogram, top), estimated tumor purity (tile plot, top), mutations in EGFR (tile plot, middle), mutational signatures (filled histogram, middle), and clinical response and clinical covariates (bottom). (b) Interaction between smoking-related mutational signatures and mutational burden. For aging/unknown vs smoking signatures p=0.0001; for aging/unknown vs APOBEC, p=0.427; for APOBEC vs smoking, p=0.158. (c) Proportion of patients with a given dominant mutational signature by clonal mutation composition. (d) Proportion of subclonal mutations in EGFR-mutant vs. EGFR-wildtype tumors (n = 57 biologically independent samples, p=0.035, unpaired two sample t-test). (e) Relationship between mutational burden and response by dominant mutational signature (for Aging/unknown dominant samples, CR/PR vs PD p=0.0044, CR/PR vs SD p=0.07, SD vs PD p=0.0597; for APOBEC/smoking dominant samples, CR/PR vs PD p=0.0023, CR/PR vs SD p=0.08, SD vs PD p=0.252). Two-sided Mann-Whitney U for (b, e), *p<0.05, **p<0.005, ns = not significant. Boxplots (b, d, e) show the median, first and third quartiles, whiskers extend to 1.5 × the interquartile range, and outlying points are plotted individually.

We also examined mutational signatures in melanoma tumors, which are dominated by exposure to ultraviolet (UV) light (S7), prior chemotherapeutic treatment with alkylating agents (S11), and other signatures not clearly associated with specific environmental exposures (S1, S5). Dominant mutational signature explained a large proportion of the variance in mutational burden and was highly correlated with intratumoral heterogeneity (Fig. 4a, Supplementary Fig. 8a). After stratifying by dominant mutational signature, no significant difference in mutational burden was observed between CR/PR vs. PD (p>0.05 for all) (Fig. 4b). The non-UV/non-alkylating group had a higher proportion of patients with PD compared to the other two groups (Fig. 4c), and was composed largely of mucosal, uveal, and acral lentiginous melanomas, though dominant mutational signature and histology did not overlap perfectly (Supplementary Fig. 8b). The observation that mutational load is not a significant predictor of response after correcting for dominant mutational signature in melanoma raises the possibility that in this cancer type, mutational burden itself may not directly mediate response but rather may serve as a proxy for an underlying biological process that both increases tumor immunogenicity and promotes accumulation of somatic mutations.

Fig. 4

Integrated analysis of mutational burden, intratumoral heterogeneity, and mutational signatures in melanoma, HNSCC, and bladder cancer

(a) Nonsynonymous mutational burden (alkylating- vs. UV-dominant, p=0.0005; UV- vs. other-dominant, p=6.079e-14) and clonal mutational burden (alkylating- vs UV-dominant, p=0.938, UV- vs. other-dominant, p=6.404e-15) stratified by dominant mutational signature in melanoma patients only. b) Association between mutational burden and response within dominant mutational signature in melanoma (CR/PR vs PD in alkylating-dominant, p=0.2; CR/PR vs PD in UV-dominant, p=0.549; CR/PR vs PD in ‘Other’ dominant, p=0.689). (c) Likelihood of response by dominant mutational signature group in melanoma. (d–e) Stacked plots of mutational burden (histogram, top), tumor purity (tile plot, top), mutational signatures (filled histogram, middle), and histology and clinical response to immune checkpoint therapy (tile plots, bottom) for HNSCC (d) and bladder cancer (e). (f) Proportion of mutations attributable to APOBEC mutational signatures (S2/13) vs. nonsynonymous mutational burden (n = 39 biologically independent samples, p=2.66e-08 for slope; p=0.848 for intercept). Symbols indicate cancer type. (g) Proportion of mutations probabilistically attributable to the APOBEC mutational signature in CR/PR vs. PD in HNSCC and bladder cancer (p=0.019, n = 39). Two-sided Mann-Whitney U for all, *p<0.05, **p<0.005, ns = not significant. Boxplots (a, b, g) show the median, first and third quartiles, whiskers extend to 1.5 × the interquartile range, and outlying points are plotted individually.

In bladder cancer and HNSCC, similar analyses demonstrated association of APOBEC-associated signatures (S2 and S13) with higher mutational burdens (p=0.002 for bladder, p=0.03 for HNSCC, two-sided Mann Whitney U) and greater likelihood of CR/PR (p=0.019) (Fig. 4d–g). APOBEC signatures have been linked with tobacco exposure[34], but viral infection and PIK3CA hotspot mutations are other potential etiologies[35]. These correlations may contribute to and/or confound the observed association between PIK3CA hotspot mutations and CR/PR described above (Fig. 2c, Supplementary Fig. 6). APOBEC mRNA expression has also been previously associated with increased PD-L1 IHC staining[36] and high tumor mutational burden[37] in urothelial carcinoma. Overall, genetic features – including mutational burden, intratumoral heterogeneity, tumor driver mutations, and mutational signatures – appear to have inter-related associations with response in this cohort. Determining which of these features ultimately drives response to immune checkpoint therapy in these patients will require further clinical and experimental study.

Copy number alterations associated with response or resistance to immune checkpoint therapy

In addition to somatic mutations, copy number alterations (CNAs) may also contribute to selective response, and CNAs affecting the interferon-γ pathway have been implicated in intrinsic resistance to immune checkpoint therapies in melanoma[17]. In this cohort, we assessed CNAs expected to interfere with interferon-γ signaling after correcting for tumor purity[26], which can strongly influence the number of called CNAs in cohorts with heterogeneous tumor purity and ploidy (Supplementary Fig. 9a–c). After correcting for tumor purity, interferon-γ-related CNA events were more infrequent than previously described, but consistent with prior studies, these events were enriched in PD vs. CR/PR (19/123 vs. 3/70; p=0.019, Fisher’s exact) (Fig. 5a). This relationship persisted within cancer types and therapy classes, though this study was insufficiently powered to detect a significant association in many of the subgroup analyses (Fig. 5b–c, Supplementary Fig. 10a–b, Supplementary Fig. 11a–b).

Fig. 5

Tumor copy number alterations associated with response to immune checkpoint therapy

(a) Amplifications and deletions of genes in the interferon-γ signaling pathway in CR/PR vs. PD across all samples (Two-tailed Fisher’s exact p=0.019, 95% CI 0.05 to 0.885, n=193), by (b) drug class (n=120 anti-CTLA4-treated, n = 65 anti-PD1/PD-L1-treated) and (c) cancer type (n = 23 bladder, 7 HSNCC, and 125 melanoma). Error bars denote proportion of CR/PR or PD with CNA +/− standard error. *p<0.05. (d) Difference in proportion of CR/PR vs. PD harboring focal cancer driver CNAs. Negative log10(p-value) for a two-tailed Fisher’s exact test for enrichment of a gene-level CNA in CR/PR vs. PD is shown on the y-axis. Genes more commonly affected by CNAs in CR/PR are shown on the right, while those more commonly deleted or amplified in PD are shown on the left. Dashed red line indicates p=0.05. (n=193 biologically independent samples) (e) Truncating mutations in PTEN by response group. (f) Prevalence and response association of truncating alterations in genes encoding SWI/SNF subunits[62]. Dark blue tiles indicate membership in either PBAF or BAF, which are complexes within the SWI/SNF family. Only genes encoding SWI/SNF subunits harboring truncating mutations in at least two patients are shown.

Next, we analyzed focal CNAs affecting 63 commonly amplified and deleted tumor suppressors and oncogenes[26] (Supplementary Tables 8–9) to nominate additional mediators of selective response or resistance. While these events were rare and no gene achieved statistical significance alone, amplifications of PAK1, YAP1, and CCND1 on chromosome 11q and amplifications of MDM2 and CDK4 on chromosome 12q were seen predominantly in PD patients (Fig. 5d), with the latter of note as CDK4/6 inhibition was recently associated with increased tumor immunogenicity[38]. Additionally, homozygous PTEN deletion occurred exclusively in patients with intrinsic resistance to immune checkpoint therapy (N=4 PD; p=0.30) (Fig. 5d), re-capitulating the observation of biallelic PTEN loss in resistance to immunotherapies in prior clinical and experimental studies [23,39]. Clonal biallelic loss of PTEN via truncating mutation was not as clearly associated with PD, although many of these were splice site mutations of uncertain biological significance (Fig. 5e). Truncating PTEN mutations were also seen in resistant tumors from patients with tumor shrinkage following anti-CTLA-4 therapy at other sites, although these patients were excluded from the main analysis due to their mixed clinical response (Supplementary Fig. 12a–b, Supplementary Table 1). Thus, both mutations and CNAs can influence response to immune checkpoint therapy by activating or suppressing pathways that interact with tumor-immune signaling, though our analyses indicate that these features are elaborately interrelated, and require significantly increased cohort size and mechanistic validation for robust interpretation.

Validation of previously described response predictors for immune checkpoint therapies

Prior studies have identified additional pathways and genes associated with response or resistance to immune checkpoint therapy, which were not identified via the unbiased analyses presented thus far, potentially due to insufficient power after correcting for multiple hypothesis testing. We attempted a focused validation of these previous findings in this cohort, beginning with an analysis of loss-of-function alterations in PBRM1, a member of the PBAF form of the SWI/SNF complex. Prior work in clear-cell renal cell carcinoma (ccRCC), a cancer type not represented in this study, demonstrated that biallelic PBRM1 loss correlates with response to anti-PD-1/L1 therapy, while mutational load does not[19]. Additional functional data in melanoma have supported this association, and implicated other related genes in the SWI/SNF family of chromatin remodelers[40]. In this cohort, truncating mutations in PBRM1 and other related epigenetic regulators in the SWI/SNF complex were also over-represented in responders, with biallelic loss events occurring almost exclusively in the tumors of CR/PR patients (Fig. 2a–b, Fig. 5f, Supplementary Fig. 13). Other mutations that have previously been described in association with acquired resistance to immune checkpoint therapies include JAK/STAT pathway mutations and alterations in antigen-presentation machinery (e.g. beta-2-microglobulin, tapasin)[41,42]. Loss-of-function mutations and deletions in these pathways were rare in this cohort (Supplementary Fig. 14a–c), perhaps reflecting differences in genetic mechanisms of intrinsic vs. acquired resistance, and further work is required to determine the functional significance of these variants.

Neoantigens associated with response to immune checkpoint therapy

Given the complexities of computationally assessing the impact of oncogenic pathway alterations and tumor mutational burden on response to immune checkpoint therapy, we lastly examined mutations that could more directly influence tumor-immune interactions by generating tumor-specific neoantigens that induce a T cell-mediated anti-tumor response[43]. We inferred neoantigens in silico[26] (Supplementary Table 10) and found that, on average, each nonsynonymous mutation generated 2.24 predicted neoantigens, with extremely high correlation (R2=0.99; p<0.00001) (Fig. 6a), making it difficult to disentangle the effects of mutational vs. neoantigen burden on response using in silico methods alone.

Fig. 6

Response-associated in silico predicted neoantigens

(a) Relationship between predicted neoantigen burden (y-axis) and nonsynonymous mutational burden (x-axis). Linear regression excludes one outlier (Pat110) (n = 249 biologically independent samples). (b) Prioritization of clinically actionable predicted neoantigens by presence in cancer driver genes and exclusive presence in CR/PR. (c) Response-associated predicted neoantigens generated by cancer driver mutations.

Recent studies have demonstrated that personalized cancer vaccines targeting neoantigens specific to an individual’s tumor can lead to durable clinical benefit alone or in concert with immune checkpoint blockade therapies[44,45]. More than 99% of predicted neoantigens in this study arose from passenger mutations, which occur throughout the exome, are frequently found in subclonal tumor populations, and are largely unique to each patient’s tumor. However, 871 predicted neoantigens were generated by driver mutations, and eight of these “driver” neoantigens occurred recurrently in CR/PR but not in PD patients, in an HLA-dependent manner (Fig. 6b–c). Additionally, as expected given the known oncogenic effects of the mutations yielding these neoantigens, these eight neoantigens were clonal in all samples, suggesting that a T cell-mediated response, if present, would target all cancer cells[14]. Thus, driver alterations can generate tumor neoantigens and may contribute to provoking an effective immune response to checkpoint blockade therapy in HLA-matched patients, though further experimental study is required to clarify the biological significance of these putative neoantigens.

Discussion

Whole exome sequencing and analysis of 249 tumors from patients treated with anti-PD-1/L1 or anti-CTLA-4 therapies suggest that genomic features beyond mutational burden, including genetic driver events, tumor heterogeneity, and mutational signatures, may affect response to immune checkpoint blockade. In this work, we combined data from multiple institutions using a standardized computational pipeline and applied a uniform and well-accepted definition of radiographic response to cancer therapy to more robustly assess genetic predictors of response to immune checkpoint therapy.

In so doing, we validated past findings and expanded their generalizability to new cancer contexts, discovered new correlative biomarkers of response using a larger sample size with more statistical power, and investigated the relationships between predictive biomarkers in enhanced detail. For instance, biallelic PTEN loss was first clinically described in acquired resistance to anti-PD-1 therapy in metastatic uterine leiomyosarcoma, and these clinical results and prior pre-clinical findings suggest it may be relevant to intrinsic resistance and mixed response in metastatic melanoma as well[23,39]. Copy number alterations leading to loss of intact interferon-γ signaling were previously noted in metastatic melanoma treated with anti-CTLA-4 agents[41], but this mechanism may also have relevance to anti-PD-1/L1-treated patients and in additional cancer types. PBRM1 loss was first implicated in increasing tumor intrinsic responsiveness to immune checkpoint therapy in ccRCC[19], but shared biology of loss of PBRM1 or ARID2 – which both encode proteins within the PBAF form of the SWI/SNF complex – in melanoma, lung cancer, bladder cancer, and HNSCC may underlie the similar response association observed here. Notably, experimental descriptions of the mechanistic underpinnings of these genetic variants in influencing response or resistance to T cell-mediated killing have been instrumental in supporting computational results[17,39,40], emphasizing the importance of cross-validation between clinically-relevant tumor variants and mechanistically-driven investigations.

Such findings show that comprehensive consideration of multiple genomic features may help place existing associations such as mutational burden in a broader biological context. In the melanoma tumors from this study, mutational burden was no longer a significant predictor of response after correcting for dominant mutational signature, a finding that warrants further experimental and translational inquiry. KRAS and EGFR mutations in lung cancer have previously been described as response- or resistance-associated, respectively, and this study demonstrates a relationship between these driver mutations and carcinogenic exposures, intratumoral heterogeneity, and mutational burden[7,33]. The global approach outlined herein will be essential in future investigations of predictors of immune checkpoint therapy response.

Power calculations suggest that combining hundreds or even thousands of clinically annotated patient samples will be necessary to reliably detect specific predictors of response to immune checkpoint therapies. While we preliminarily assessed response predictors specific to anti-CTLA-4 vs. anti-PD-1/L1 therapies, or within a cancer type, these findings are biased by available data from clinically-annotated cohorts; anti-CTLA-4 therapies were used predominantly in melanoma and anti-PD-1/L1 therapies dominated most other cancer types. Further studies more directly comparing therapy classes within the same tumor histology and vice versa will be necessary, as will consideration of response predictors for combinations of checkpoint inhibitors with or without targeted or cytotoxic chemotherapies. Exceptions to single-feature genomic associations between this study and prior works likely can be explained by complex context-dependent effects, and emphasize the need for caution to avoid over-interpretation of results.

While sample size and cohort heterogeneity remain major limitations of this work, this study describes a path forward for gathering insights from multiple clinically annotated patient cohorts. Our work advances hypotheses of biological mechanisms, suggests clinically relevant biomarkers, and highlights the importance of further, larger studies to reliably and robustly identify biomarkers of response and intrinsic resistance to immune checkpoint blockade.

Online Methods

Clinical cohort consolidation

Patients from the Dana-Farber Cancer Institute with metastatic bladder cancer, HNSCC, lung cancer, and melanoma treated with anti-PD-1, anti-PD-L1, anti-CTLA-4, or a combination of these therapies were identified, and pre-treatment tumor tissue and matched germline blood were obtained for genetic sequencing. These studies were conducted in accordance with the Declaration of Helsinki and approved by the Dana-Farber Cancer Institute Institutional Review Board (Protocols 11-104, 05-042, 02-180, 09-472, 02-021, 15-330). Electronic medical charts were reviewed to assess best response by RECIST (v1.1), duration of PFS, duration of OS, patient demographic characteristics, and other relevant clinical details (e.g. smoking history). “Current/former” smokers were those who reported >5 pack-year (packs per day * years smoking) of tobacco use. “Never” smokers were those who reported ≤5 pack-years[46]. Clinical information from studies conducted outside the Dana-Farber Cancer Institute were downloaded from online supplemental materials[7,8]. Where OS was not reported in these studies, OS was censored at PFS.

Response stratification

Patients were divided into objective responder (OR) and non-responder (NR) groups according to three published response metrics utilizing best response by RECIST, duration of PFS, and duration of OS[10,11,12]. For patients in the MSKCC melanoma cohort, all non-responder patients were presumed to have PD as their best response by RECIST. Tumor samples from progressing lesions from patients who had clinical benefit from immune checkpoint therapy were excluded from the main analysis. Tumors from patients with PFS or OS below 30 days were excluded from analysis, as these patients may have had disease that was too advanced to experience clinical benefit from immune checkpoint therapy.

DNA extraction and sequencing

For samples newly sequenced from the Dana-Farber Cancer institute, DNA extraction from formalin-fixed, paraffin-embedded (FFPE) tumor blocks was performed as previously described[47]. Exome sequencing and data processing to produce a BAM file was performed using established analytical pipelines at the Broad Institute[23].

DNA sequencing quality control

These data were combined with raw sequencing data (BAM files) from previously published cohorts of tumor/normal sequencing from patients with metastatic melanoma[8,10], lung cancer[7,25], anal cancer[24], and sarcoma[23]. All 314 samples with tumor and germline sequencing data and clinical annotations were processed through standard quality control pipelines. Samples with poor sequencing coverage (tumor mean target coverage < 25x, normal mean target coverage < 15x) or high sample contamination[48] in tumor or normal tissue were excluded (Fig. 1a, Supplementary Table 1). Additionally, samples with germline sequencing from adjacent normal tissue were assessed for tumor-in-normal contamination using deTiN (https://github.com/broadinstitute/deTiN), and excluded if normal tissue contained ≥1% tumor nuclei (Fig. 1a, Supplementary Table 1). After mutation calling and somatic copy number alteration assessment, tumors with estimated purity below 10% were also excluded[49] (Supplementary Fig. 1a–b). These quality control measures were taken because high sample contamination and low tumor purity can lead to systematic under-calling of somatic SNPs and CNAs and interference with accurate assessment of tumor mutational burden and identification of response-associated molecular features[50].

Whole exome analysis

Somatic nucleotide polymorphisms (SNPs) were identified by MuTect[50], with computational filtering of artifacts introduced by DNA oxidation during sequencing[51] or FFPE-based DNA extraction using a filter-based method. Strelka[52] was applied to detect small insertions and deletions (indels). Annotation of identified variants was done using Oncotator (http://www.broadinstitute.org/cancer/cga/oncotator). Mutational clonality was estimated by ABSOLUTE, which uses allelic fraction of called mutations and allelic copy number information to determine mutational clonality and overall tumor purity and ploidy[49]. Clonal mutations were defined as those with estimated cancer cell fraction (CCF) of 1 or those whose probability of being clonal exceeded the probability of being subclonal. Nonsynonymous mutational burden was normalized by megabases covered at adequate depth to detect variants with 80% power using MuTect given estimated tumor purity by ABSOLUTE (Supplementary Table 1). Number of bases covered at a given depth threshold in the tumor sample was determined using the GATK DepthOfCoverage module. Putative driver mutations were collected using cBioPortal standards[53,54], including both 3D hotspots in tumor suppressors and oncogenes[55] and any loss-of-function variant in a tumor suppressor.

For copy number analysis, copy ratios were calculated for each captured target by dividing the tumor coverage by the median coverage obtained in a set of reference normal samples. The resulting copy ratios were segmented using the circular binary segmentation algorithm[56]. Segments were considered amplified or deleted if the |log2(copy ratio)| exceeded 0.5 (Supplementary Fig. 9a). For samples with uniform sample purity and ploidy, this definition of amplifications and deletions is adequate to detect copy number alterations. However, in this sample of tumors, called copy number alterations were heavily influenced by sample purity (Supplementary Fig. 9a). The |log2(copy ratio)| > 0.5 definition was insufficiently sensitive in low-purity samples, which have artificially depressed copy ratios due to a high proportion of sequencing reads from normal tissue. Conversely, this definition may be excessively noisy in high-purity tumors.

Thus, to correct segment copy ratios for sample purity, segment copy ratios were re-scaled by sample purity and ploidy with values derived from ABSOLUTE. Segments were considered amplified or deleted if the |log2(purity-corrected copy ratio)| exceeded 0.5. Specifically, the purity-corrected copy ratio was derived by dividing the purity-corrected (or rescaled) total copy number (rCN) for a given segment by the sample ploidy. This procedure yielded improved false-negative rates in low-purity tumors and false-positive rates in high-purity tumors, such that the proportion of a tumor genome considered amplified or deleted was less closely associated with tumor purity (Supplementary Fig. 9b).

However, while using a |log2(purity-corrected copy ratio)| > 0.5 definition of deletions and amplifications was effective for detecting large CNAs with high sensitivity, it did not provide adequate specificity for detecting focal events that would be more likely to be genomic driver CNAs for a tumor. Thus, we applied a previously described concept called focality[57] to identify CNAs that were large outliers in copy ratio, representing either homozygous deletions expected to completely eliminate tumor expression of a given gene or amplifications expected to greatly over-express a gene. In this process, the rCN from ABSOLUTE was used as input. For each segment in a tumor genome, the focality was calculated by considering the fraction of a sample’s genome with lower rCN than that segment (for amplified regions) or higher rCN (for deleted regions). Segments were considered deleted if their rCN was < 0.25 and their focality was > 0.995. Segments were considered amplified if focality exceeded 0.98 − 0.2 × log2 (rCN/5), and highly amplified if focality exceeded 0.98 − (1/7) × log2 (rCN/7). The results from this analysis are displayed in Supplementary Fig. 9c, where far fewer segments meet criteria for being called as amplifications or deletions under this focality-based definition. This focality-based definition of CNA was applied in Fig. 5 as well as Supplementary Figs. 10–12 to demonstrate putative driver CNAs in genes in the interferon-gamma signaling pathway, PTEN, and SWI/SNF family of chromatin regulators. Genes were considered amplified or deleted if all or part of the gene was in a segment with a called copy number alteration using this focality-based definition. For the interferon-gamma analysis, samples were defined as having a copy number alteration affecting interferon signaling if any of the regular interferon genes were deleted or any of the four interferon pathway inhibitors (SOCS1, SOCS3, PIAS1, and PIAS4) harbored high amplifications, as previously described[17]. All gene-level CNAs of interest were manually reviewed.

A similar focality procedure was applied to allelic copy number calls from ABSOLUTE to determine heterozygous deletions and amplifications to identify loss-of-heterozygosity events for Supplementary Figs. 12–13.

Mutational signature deconvolution was conducted using a non-negative matrix factorization technique as previously described[30]. Mutational signatures were chosen from those previously described in COSMIC (http://cancer.sanger.ac.uk/cosmic/signatures)(Supplementary Table 6). The vectors for the commonly observed mutational signatures for each cancer type were used as input for inference of their contribution to observed mutations. Thus, for example, the signatures selected for melanoma pertained to UV exposure, prior alkylating agents, and other exposures, while those used for mutational signature deconvolution in lung cancer included tobacco exposure.

For neoantigen prediction, the 4-digit HLA type for each sample was inferred using Polysolver[58]. Putative neoantigens were predicted for each patient by defining all novel amino acid 9-mers and 10-mers resulting from each somatic nonsynonymous point mutation and determining whether the predicted binding rank—a proxy for predicted binding affinity to the patient’s germline HLA alleles—was < 2%. Strong binders had rank < 0.5%, while weak binders had rank between 0.5% and 2% using NetMHCpan (v3.0)[59,60,61].

Statistical analysis

Assessment of enrichment of binary molecular features (e.g. wildtype or mutant gene; CNA present or absent) with response (CR/PR vs. PD) was done with Fisher’s exact tests. Assessment of difference in means or medians for a continuous variable between two response groups was done with the non-parametric Mann-Whitney U test unless otherwise specified. Correction for multiple hypothesis testing was done controlling for false discovery rate (FDR) by the Benjamini-Hochberg method, unless otherwise noted. Receiver operating characteristic (ROC) analyses were done using the pROC and Epi packages in R.

For the power calculation analysis (Fig. 2e), correction for multiple hypothesis testing was modeled with a Bonferroni correction over the 116 genes with known cancer driver status assessed previously. Response rates were set at 40% CR/PR vs. 60% PD, which is a generous estimate for response rate in an unselected population. P-values were calculated using Fisher’s exact tests comparing the prevalence of mutations in a given gene in CR/PR vs. PD.

Alpha level for all comparisons was 0.05 unless indicated otherwise. All statistical analyses were done in R (v.3.3.2).