Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma

Liu, David; Schilling, Bastian; Liu, Derek; Sucker, Antje; Livingstone, Elisabeth; Jerby-Arnon, Livnat; Zimmer, Lisa; Gutzmer, Ralf; Satzger, Imke; Loquai, Carmen; Grabbe, Stephan; Vokes, Natalie; Margolis, Claire A.; Conway, Jake; He, Meng Xiao; Elmarakeby, Haitham; Dietlein, Felix; Miao, Diana; Tracy, Adam; Gogas, Helen; Goldinger, Simone M.; Utikal, Jochen; Blank, Christian U.; Rauschenberg, Ricarda; von Bubnoff, Dagmar; Krackhardt, Angela; Weide, Benjamin; Haferkamp, Sebastian; Kiecker, Felix; Izar, Ben; Garraway, Levi; Regev, Aviv; Flaherty, Keith; Paschen, Annette; Van Allen, Eliezer M.; Schadendorf, Dirk

doi:10.1038/s41591-019-0654-5

Download PDF

Article
Open access
Published: 02 December 2019

Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma

David Liu^1,2^na1,
Bastian Schilling ORCID: orcid.org/0000-0001-8859-4103^3,4,5^na1,
Derek Liu^1,2,6,
Antje Sucker^4,5,
Elisabeth Livingstone^4,5,
Livnat Jerby-Arnon²,
Lisa Zimmer^4,5,
Ralf Gutzmer⁷,
Imke Satzger⁷,
Carmen Loquai⁸,
Stephan Grabbe⁸,
Natalie Vokes^1,2,
Claire A. Margolis^1,2,
Jake Conway^2,6,
Meng Xiao He^2,6,9,
Haitham Elmarakeby^1,2,
Felix Dietlein^1,2,
Diana Miao ORCID: orcid.org/0000-0003-1754-2319^1,2,6,
Adam Tracy²,
Helen Gogas¹⁰,
Simone M. Goldinger¹¹,
Jochen Utikal ORCID: orcid.org/0000-0001-5316-0241^12,13,
Christian U. Blank ORCID: orcid.org/0000-0002-7945-5846¹⁴,
Ricarda Rauschenberg^15,16,17,
Dagmar von Bubnoff¹⁸,
Angela Krackhardt^5,19,
Benjamin Weide²⁰,
Sebastian Haferkamp²¹,
Felix Kiecker²²,
Ben Izar^1,2,
Levi Garraway²³,
Aviv Regev ORCID: orcid.org/0000-0003-3293-3158²,
Keith Flaherty ORCID: orcid.org/0000-0002-3402-0478²⁴,
Annette Paschen^4,5,
Eliezer M. Van Allen ORCID: orcid.org/0000-0002-0201-4444^1,2 &
…
Dirk Schadendorf ORCID: orcid.org/0000-0003-3524-7858^4,5

Nature Medicine volume 25, pages 1916–1927 (2019)Cite this article

51k Accesses
466 Citations
135 Altmetric
Metrics details

Subjects

An Author Correction to this article was published on 19 June 2020

This article has been updated

Abstract

Immune-checkpoint blockade (ICB) has demonstrated efficacy in many tumor types, but predictors of responsiveness to anti-PD1 ICB are incompletely characterized. In this study, we analyzed a clinically annotated cohort of patients with melanoma (n = 144) treated with anti-PD1 ICB, with whole-exome and whole-transcriptome sequencing of pre-treatment tumors. We found that tumor mutational burden as a predictor of response was confounded by melanoma subtype, whereas multiple novel genomic and transcriptomic features predicted selective response, including features associated with MHC-I and MHC-II antigen presentation. Furthermore, previous anti-CTLA4 ICB exposure was associated with different predictors of response compared to tumors that were naive to ICB, suggesting selective immune effects of previous exposure to anti-CTLA4 ICB. Finally, we developed parsimonious models integrating clinical, genomic and transcriptomic features to predict intrinsic resistance to anti-PD1 ICB in individual tumors, with validation in smaller independent cohorts limited by the availability of comprehensive data. Broadly, we present a framework to discover predictive features and build models of ICB therapeutic response.

Mutated processes predict immune checkpoint inhibitor therapy benefit in metastatic melanoma

Article Open access 19 September 2022

Andrew Patterson & Noam Auslander

Patients deriving long-term benefit from immune checkpoint inhibitors demonstrate conserved patterns of site-specific mutations

Article Open access 07 July 2022

Daniel R. Principe

Pan-cancer analysis of longitudinal metastatic tumors reveals genomic alterations and immune landscape dynamics associated with pembrolizumab sensitivity

Article Open access 26 August 2021

S. Y. Cindy Yang, Scott C. Lien, … Trevor J. Pugh

Main

While ICB has resulted in durable clinical response in multiple tumor types^{1,2,3,4,5,6,7}, only a subset of patients respond, and predictors of response are not fully characterized. Both tumor-intrinsic and tumor-extrinsic biomarkers of response and resistance to ICB in melanoma have been proposed, including tumor mutational burden (TMB) and neoantigen load^8,9,10,11, immunohistological detection of PD-L1 and CD8¹² and genetic alterations affecting antigen presentation^13,14, interferon (IFN)-γ signaling pathways¹⁵, alternative survival and proliferation pathways^13,16,17 and aneuploidy^18,19. Gene expression signatures expressed in tumors²⁰ and the tumor immune microenvironment²¹ have also been implicated. However, these observations have often been made in preclinical models or in small clinical cohorts without validation in larger, independent cohorts of patients with melanoma. Furthermore, whether these observations are exclusive to a specific ICB regimen (that is, anti-PD1, anti-CTLA4 or a combination of these) is incompletely characterized. Broadly, the expanding suite of pathways that has been invoked to mediate selective ICB response in melanoma indicates that integrated systems biology models to predict response and survival are necessary, but these have yet to be well developed.

Clinically, the optimal role of anti-CTLA4 in conjunction³ or sequentially²² with anti-PD1 ICB is unclear. Understanding the differential biology underlying the response to anti-PD1 ICB in tumors with and without previous anti-CTLA4 therapy may inform the rational design of combination therapies and optimize therapy selection for individual patients.

Thus, we performed an integrative study employing genomic, transcriptomic and clinical data from a comprehensively clinically annotated and sequenced cohort of 144 patients with advanced melanoma undergoing anti-PD1 ICB with and without previous anti-CTLA4 ICB to discover biomarkers of response and resistance, and develop clinically applicable parsimonious predictive models.

Results

Genomic and clinical cohort characteristics and melanoma subtypes

We identified 206 patients diagnosed with advanced melanoma and treated with anti-PD1 ICB, and performed whole-exome sequencing (WES) on matched pretreatment tumor samples and normal tissue²³, and whole-transcriptome sequencing (RNA-seq) on available pretreatment tumor tissue. After quality control (Methods), WES data from 144 patients and RNA-seq data from 121 patients were available for final evaluation (Extended Data Fig. 1). Best objective response (BOR) to anti-PD1 ICB using RECIST (v.1.1) criteria (Methods) included 45% with progressive disease (PD), 14% with stable disease (SD), 3% with mixed response (MR), 26% with partial response (PR) and 12% with complete response (CR; Fig. 1a), for an overall response rate of 38%. Overall, 73% were cutaneous melanomas, 13% were of occult origin, 7% were mucosal and 7% were acral in origin. A total of 44% (n = 64) of patients had previous treatment with ipilimumab, whereas 56% (n = 80) were naive to ipilimumab. The median follow-up for survival was 29.9 months. Other clinical characteristics are detailed in Table 1.

**Fig. 1: Cohort genomic and clinical characteristics and association of TMB with response.**

Table 1 Cohort clinical characteristics

Full size table

Overall the median nonsynonymous TMB was 6.5 mutations per Mb (250.5 mutations per exome), with an interquartile range (IQR) of 2.0–14.4 mutations per Mb (77.75–578.5 mutations per exome). Overall, 39% of tumors had BRAF mutations, 30% had NRAS mutations and 17% had NF1 mutations (Fig. 1a). The median tumor purity (the proportion of sample DNA from tumor cells) was 0.67 (IQR 0.46–0.83) and the median tumor heterogeneity (the proportion of subclonal mutations) was 0.17 (IQR 0.12–0.25). The median purity-corrected tumor ploidy (Methods) was 2.15 (IQR 2.01–3.12), with 38% of tumors inferred to have genome doubling, consistent with previous reports²⁴. The predominant mutational signature in most tumors was related to ultraviolet (UV) exposure²⁵ (69% related to UV, 3% related to alkylating chemotherapy²⁵ and 28% related to another predominant mutational signature, mostly associated with aging²⁵; Fig. 1a). Individual tumor characteristics are detailed in Supplementary Table 1.

To discover differential features associated with response, we compared clinical responders (n = 55) to progressors (n = 65), excluding patients with SD (n = 20) and MR (n = 4) as the BOR. Overall survival (OS) and progression-free survival (PFS) were significantly different between these groups (log-rank P < 0.00001 for both comparisons; Extended Data Fig. 2a,b).

TMB was higher in responders than in progressors (Mann–Whitney–Wilcoxon (MWW), P = 0.026; Fig. 1b), but there was substantial overlap between responders and progressors. We hypothesized that the relationship between response and TMB might further be confounded by melanoma subtype. TMB was significantly different between different melanoma subtypes (Kruskal–Wallis, P = 2.4 × 10⁻⁵; Fig. 1c), with cutaneous and occult melanomas having similar and higher TMB than acral and mucosal melanomas²⁶ (median of 297.5 versus 58, nonsynonymous mutations; MWW, P = 1.1 × 10⁻⁶), with a higher response rate (~40% versus ~20%; Fisher’s exact test, P = 0.06; Fig. 1d). When stratified by melanoma subtype, responders did not have significantly higher TMB than nonresponders (Fig. 1e), and, in multivariate logistic regression adjusting for melanoma subtype, TMB was not a significant predictor (P = 0.24). Strikingly, responders with mucosal or acral melanoma had a lower TMB than progressors with cutaneous or occult melanoma (MWW, P = 0.03; Fig. 1f), suggesting that disease subtype confounds the association between TMB and response to anti-PD1 therapy.

Genomic and transcriptomic features associated with response

Higher tumor purity and heterogeneity were associated with progression (MWW, P = 0.04 and P = 0.02, respectively; Fig. 2a,c), whereas ploidy was lower in progressors (MWW, P = 0.04; Fig. 2b). The proportion of the tumor genome with copy number alterations (CNAs) trended toward being higher in patients with PD (MWW, P = 0.09; Extended Data Fig. 2c).

Given these observations, we performed an unbiased analysis for single-gene predictors of response to anti-PD1 ICB across all mutated genes detected in this cohort. After multiple-hypothesis test correction, no genes were significantly associated with response or resistance to therapy (Extended Data Fig. 3a and Supplementary Table 2), highlighting the large sample sizes needed for adequate power to detect these associations²⁷. We observed only rare mutations in major histocompatibility complex class I (MHC-I) antigen-presentation genes (TAP1, TAP2, B2M, HLA-A, HLA-B and HLA-C)^14,28. Mutations in SERPINB3 or SERPINB4²⁹ were not associated with response (Fisher’s exact test, P = 0.51 and P = 1.0, respectively). Loss of heterozygosity (LOH) in B2M¹⁴ was found in 9 of 55 (16%) responders and 16 of 65 (25%) progressors (odds ratio (OR) = 0.6) but was not significantly associated with resistance (Fisher’s exact test, P = 0.37). LOH in HLA-A, HLA-B or HLA-C²⁸ was not associated with response or resistance to therapy (Fisher’s exact test, P = 0.52, P = 0.57 and P = 0.84, respectively). Confirming previous findings, LOH of JAK1^13,30,31 was associated with resistance (OR = 0.33; Fisher’s exact test, P = 0.02). Biallelic CDKN2A alteration²⁷ was found in 15 of 55 (27%) responders and 25 of 65 (38%) progressors (OR = 0.6), but was not significant (Fisher’s exact test, P = 0.24).

We also performed an unbiased analysis of the association between focal gene amplifications and response to therapy. While no gene amplification was significant after multiple-hypothesis test correction (Extended Data Fig. 3b and Supplementary Table 3), amplification of TAP2, an integral part of the MHC-I antigen-loading pathway, was found exclusively in responders to therapy (n = 6; Fisher’s exact test, P = 0.008; Fig. 2d). TAP2 is located at 6p21 in a region encoding both MHC-I and MHC-II human leukocyte antigen (HLA) loci, and four out of six TAP2 amplifications were associated with larger amplifications across the region, while two out of six amplifications were more focal (Extended Data Fig. 4). Notably, tumors with amplifications in this region encompassing the MHC-I-related HLA-A, HLA-B and HLA-C genes (a region of approximately 1.5 Mb; n = 6) were also associated with response to therapy (Fisher’s exact test, P = 0.008; Fig. 2e), with four out of six also having amplifications in TAP2. Altogether, eight patients had amplification of either TAP2 or HLA-A, HLA-B or HLA-C amplification (Fig. 2f and Extended Data Fig. 4), and were exclusively responders (Fisher’s exact test, P = 0.001).

We then examined the expression of antigen-presentation molecules and their association with response. Interestingly, expression of all 13 MHC-II-associated HLA genes was higher in responders (collective two-sided binomial test, P = 0.0002; Fig. 2g), with four genes (HLA-DMA, HLA-DMB, HLA-DOB and HLA-DOB) individually passing a statistical significance threshold (MWW, P < 0.05; Supplementary Table 4). MHC-I antigen-presentation genes all trended toward having higher expression in responders (collective two-sided binomial test, P = 0.02; Fig. 2g and Supplementary Table 4), but none passed the statistical significance threshold.

To examine pathways differentially enriched in responders versus progressors, we next performed unbiased gene set enrichment analysis (GSEA³²; Methods) using the Hallmark gene sets³³. A total of 24 pathways were enriched (false discovery rate (FDR), q < 0.1) in responders, and 5 of the top 6 enriched pathways were immune related, including IFN-γ response, genes involved in allograft rejection, complement, the inflammatory response and interleukin (IL)-6–JAK–STAT3 signaling (Supplementary Table 5). No pathways were significantly enriched in progressors.

We further evaluated various transcriptomic signatures^{18,19,21,34,35,36,37,38,39,40} that had been proposed and demonstrated in various settings to be associated with response to immunotherapy (Methods), but we found no significant differences (P < 0.05) in these signatures between responders and progressors within our cohort (Supplementary Table 6).

Immune infiltrate has been associated with response to immunotherapy across multiple cancer types and immune therapies^12,34. We inferred the absolute level of immune infiltrate within each sample using an immune deconvolution algorithm (CIBERSORT⁴¹ using the LM22 signature matrix). We found no significant difference in the total immune infiltrate or abundance of individual immune cell subsets in responders versus progressors (Supplementary Table 7). We also generated profiles of expression of immune cell subset signatures derived from single-cell analyses⁴² and found that the expression of multiple signatures was significantly higher (unadjusted MWW, P < 0.05) in responders versus progressors, including for signatures of overall immune infiltrate, T cells, B cells, macrophages, CD8⁺ cytotoxic exhausted T cells and CD4⁺ exhausted T cells (Supplementary Table 7). Although the strength of association differed by deconvolution method, both approaches generally agreed on the direction of association, providing evidence for a moderate association of immune infiltrate with the response to anti-PD1 ICB.

Correlations between molecular features

To understand the relationship between predictors of response, we performed hierarchical clustering of the correlation coefficients between associated predictors (Fig. 2h). Clustered features are correlated and may reflect the same underlying biology, whereas separate clusters may reflect independent feature categories. A large immune-related cluster with MHC-II- and MHC-I-related subclusters was observed, with a separate independent cluster of mutation- and neoantigen-load-related features, suggesting independent feature categories. Tumor purity was negatively correlated with the immune cluster, suggesting that low tumor purity may be a proxy for higher immune infiltrates within the tumor sample. Other features, including tumor heterogeneity and ploidy, were independent from these two clusters. Extending this analysis to previously hypothesized signatures and Hallmark gene set signatures (Extended Data Fig. 5), immune activity signatures, including signatures of cytolytic²¹ and cytotoxic¹⁹ activity, IFN-γ and T effector cells^18,35, immune chemokines³⁸ and single-cell-derived immune cell signatures⁴², clustered together. A tumor-intrinsic resistance program signature⁴², signatures of T cell dysfunction and exclusion³⁷ and comparative immune-checkpoint gene expression³⁴ were distinct from the immune cluster. Overall, these findings suggested that multiple previously hypothesized predictors of ICB response reflect the same underlying biological state and additional independent classes of predictors exist that may provide additional predictive power.

Previous exposure to anti-CTLA4 ICB

Our cohort contained patients with previous exposure to ipilimumab in anti-CTLA4 ICB (n = 60) and patients who were naive to ipilimumab (n = 84; Fig. 3a). Despite the groups having similar response rates (Fig. 3b), we hypothesized that these two groups might have differential predictors of response and resistance to anti-PD1 ICB. We performed a focused analysis of patients who were treated with ipilimumab and biopsied after treatment (n = 44 with WES and n = 34 with RNA-seq) versus patients who were naive to ipilimumab (n = 84 with WES and n = 71 with RNA-seq). A composite, rank-based score of MHC-II HLA expression (single-sample GSEA (ssGSEA)⁴³; Methods) was higher in responders than in progressors in the overall cohort and the subgroup treated with ipilimumab, but was not significantly different in the subgroup that was naive to ipilimumab (Fig. 3c–e; MWW, P = 0.03, P = 0.03 and P = 0.31, respectively). We found very similar results in the largest available independent validation cohort⁴⁴ with information on previous ipilimumab treatment (Fig. 3f–h), although the difference was not significant in this smaller cohort (n = 32 patients, 15 of whom were treated with ipilimumab and 17 of whom were naive for ipilimumab).

**Fig. 3: Differential predictors of response and progression in ipilimumab-treated tumors versus ipilimumab-naive tumors.**

We next examined the association of TMB, purity, ploidy and heterogeneity with response stratified by previous ipilimumab therapy (Extended Data Fig. 6a). Unlike previous studies^44,45, we found no specific association of a higher TMB with response in the ipilimumab-naive versus ipilimumab-treated subgroup (MWW, P = 0.15, both). However, higher heterogeneity and lower ploidy were associated with progressors only in the ipilimumab-naive subgroup (MWW, P = 0.06 and P = 0.004, respectively).

We analyzed the differential expression of specific immune-related genes in responders versus progressors in ipilimumab-treated and ipilimumab-naive subgroups and found that higher expression of various immune-related pathways distinguished responders from progressors in ipilimumab-treated but not ipilimumab-naive subgroups (all P values are unadjusted). Examples included the leukocyte chemoattractants CXCL9 and CXCL10 and their receptor CXCR3 (MWW, P = 0.05, P = 0.08 and P = 0.02, respectively, in the ipilimumab-treated subgroup), CD3D (MWW, P = 0.02), B cell markers CD19 (MWW, P = 0.04) and CD20 (MS4A1; MWW, P = 0.002) and macrophage marker CD163 (MWW, P = 0.03). Interestingly, CD4, FOXP3 and CTLA4 also followed this pattern of higher expression in responders in the ipilimumab-treated subgroup (MWW, P = 0.06, P = 0.06 and P = 0.008, respectively), but CD8A and CD8B had less evidence of association with response in either ipilimumab-treated (MWW, P = 0.17 and P = 0.27, respectively) or ipilimumab-naive (MWW, P = 0.93 and P = 0.49, respectively) subgroups. Expression of TAP2 was higher (MWW, P = 0.02) in responders than in progressors in the ipilimumab-treated subgroup but not in the ipilimumab-naive subgroup (MWW, P = 0.98). In contrast, TGFB2 expression was higher in progressors in the ipilimumab-naive subgroup (MWW, P = 0.002) but not in the ipilimumab-treated subgroup (MWW, P = 0.43). The complete set of gene expression comparisons in the overall, ipilimumab-treated and ipilimumab-naive cohorts is in Supplementary Table 4.

We then repeated GSEA to examine the pathways differentially enriched in responders versus progressors between ipilimumab-treated and ipilimumab-naive subgroups. The most differentially enriched pathways were related to immune response: the IFN-γ and IFN-α responses were significantly enriched in responders in the ipilimumab-treated subgroup (FDR, q < 0.0001, both), but not in the ipilimumab-naive subgroup (FDR, q = 0.13 and q = 0.996, respectively; Fig. 3i). Using permutation testing (Methods), we found a nonsignificant empiric P value of 0.183 and 0.18, respectively, for this difference in enriched pathways in these subgroups in our discovery cohort. However, we repeated the analysis in an independent validation cohort⁴⁶ and found similar results (Fig. 3j). Complete GSEA results are provided in Supplementary Table 5.

To further dissect the impact of MHC-II expression on patient response, we stratified the cohort into patients with high and low MHC-II expression (ssGSEA, median split). In the overall cohort, low MHC-II expression was associated with primary resistance (Fig. 4a; Fisher’s exact test, P = 0.01; OR = 2.9, 95% confidence interval (CI) 1.3–6.5), but this association was largely driven by the ipilimumab-treated subgroup (Fig. 4b; Fisher’s exact test, P = 0.02; OR = 9.9, 95% CI 1.5–63.7), with a nonsignificant association in the ipilimumab-naive subgroup (Fig. 4c; Fisher’s exact test, P = 0.32; OR = 1.9, 95% CI 0.7–5.1). A formal interaction test of previous ipilimumab treatment status with MHC-II expression for predicting response was consistent with a subgroup-specific effect of low MHC-II expression on the ipilimumab-experienced subgroup (OR = 0.20, 95% CI 0.02–1.64), although this was nonsignificant (P = 0.13) in this small cohort.

**Fig. 4: Progression versus response by immune infiltrate and MHC-II HLA expression stratified by ipilimumab treatment.**

A similar analysis of estimated total immune infiltrate levels⁴¹ (with tumors split by the median into high and low groups) showed that low immune infiltrate was significantly associated with intrinsic resistance in the ipilimumab-treated subgroup (Fig. 4e; Fisher’s exact test, P = 0.02; OR = 12.9, 95% CI 1.3–125.8), but not in the ipilimumab-naive subgroup (Fig. 4f; Fisher’s exact test, P = 1.00; OR = 0.99, 95% CI 0.4–2.6) or overall cohort (Fig. 4d; Fisher’s exact test, P = 0.17; OR = 1.8, 95% CI 0.8–3.9), indicating a subgroup-specific association in the ipilimumab-experienced subgroup of low immune infiltrate with resistance to therapy (OR = 0.08, 95% CI 0.007–0.97; P = 0.047).

Taken together, these findings suggest that evidence of immune response in the tumor microenvironment at the time of progression following anti-CTLA4 ICB is a necessary but not sufficient marker for response to anti-PD1 ICB therapy; patients without immune response to anti-CTLA4 ICB are very likely to also be intrinsically resistant to anti-PD1 ICB, highlighting a high-risk and poor-prognosis subgroup of patients (Fig. 4g).

Integrative predictive modeling of primary resistance

Patients with primary resistance to anti-PD1 ICB have poor survival (Extended Data Fig. 7), and the ability to predict these patients would enable individualized management regimens (for example combination ICB) to improve outcomes. Thus, we set out to develop parsimonious predictive models integrating clinical, genomic and transcriptomic features to predict PD (primary resistance) versus non-PD (CR, PR, SD and MR) and developed separate predictive models in ipilimumab-treated and ipilimumab-naive subgroups.

In the ipilimumab-treated group (n = 34 with WES and RNA-seq), there were 16 patients who had PD and 18 who had non-PD. Using a forward-selection approach to choose the features of a parsimonious predictive model (Methods), low MHC-II HLA expression was most strongly predictive of PD (Fig. 5a) and was correlated with MHC-I HLA, IFN-α and IFN-γ response pathway scores. In the final multivariate model, high MHC-II expression, low lactate dehydrogenase (LDH; below the median of 247 U l⁻¹) and the presence of lymph node metastases were independent predictors of non-PD (P = 0.03, P = 0.02 and P = 0.04, respectively, Supplementary Table 8 and Extended Data Fig. 8a,c), and the model had an area under the curve (AUC) of 0.90 in our discovery cohort (fivefold cross-validation mean AUC = 0.83; empiric P < 0.001; Fig. 5b, Extended Data Fig. 8f and Methods). Notably, TMB did not significantly improve model fit (log-likelihood ratio, P = 0.10), was not an independently predictive feature (P = 0.18) and did not meet Bayesian information criteria (BIC; Extended Data Fig. 8c) when added to the model.

**Fig. 5: Integrative predictive modeling of intrinsic resistance to anti-PD1 ICB.**

In the ipilimumab-naive group, there were 34 patients with PD and 41 with non-PD. In the multivariate predictive model, higher heterogeneity, lower ploidy and higher purity were independently predictive of PD (P = 0.025, P = 0.014 and P = 0.046, respectively; Fig. 5c, Supplementary Table 9 and Extended Data Fig. 8b,d), with an AUC of 0.77 (tenfold cross-validation mean AUC of 0.73; empiric P = 0.036; Fig. 5d, Extended Data Fig. 8e and Methods). TMB did not significantly improve the model fit (log-likelihood ratio test, P = 0.63), was not an independently predictive feature (P = 0.63) and did not meet BIC criteria (Extended Data Fig. 8d) when added to this multivariate model. Further, each model’s performance was specific to its subgroup (P = 0.004 and P = 0.018 for the interaction between ipilimumab-experienced and ipilimumab-naive model scores, respectively, and previous ipilimumab therapy). Each model had poor performance when applied to the opposing subgroup (AUC = 0.49 and AUC = 0.54, respectively; Extended Data Fig. 9a,b), suggesting that previous ipilimumab treatment status may stratify the appropriate predictors and predictive models to be applied in these subgroups.

We attempted to validate our models in independent cohorts but were limited by a lack of publicly available cohorts with all the molecular features and relevant clinical data on previous ipilimumab therapy and biopsy timing used in our integrated model. In a limited validation, we tested predictive models incorporating individual features where data were available in an independent validation cohort⁴⁴ (Methods), and found concordant predictions of primary resistance with low MHC-II HLA expression in ipilimumab-treated tumors and higher heterogeneity in ipilimumab-naive tumors (Extended Data Fig. 10), although neither of these predictors was significant in this small cohort (empiric P = 0.21 and P = 0.066, respectively; Methods).

In an exploratory analysis, we predicted PD in response to anti-PD1 ICB in tumors where the tumor biopsy was taken before ipilimumab treatment (n = 15; Fig. 3a) using our ipilimumab-naive predictive model. Of the eight tumors with PD in response to anti-PD1 ICB, five had the highest model scores (Extended Data Fig. 9c), with an overall AUC of 0.71. Interestingly, of the three poorly discriminated tumors with low model scores but PD responses, one (patient 82) was from a brain metastasis and one (patient 80) was an acral melanoma.

Finally, splitting the cohort into subsets with high and low model scores (split by the median), we found large differences in OS and PFS in both the ipilimumab-treated subgroup (Fig. 5e; median PFS: 38.1 months versus 2.8 months, log-rank P = 0.0001; median OS: unreached versus 9.0 months, log-rank P = 0.001) and ipilimumab-naive subgroup (Fig. 5f; median PFS: 24.7 months versus 3.1 months, log-rank P = 0.005; median OS: unreached versus 15.0 months, log-rank P < 0.0001).

Discussion

In this study, we analyzed a uniformly clinically annotated cohort of patients with advanced melanoma treated with anti-PD1 ICB monotherapy for whom WES and RNA-seq data were available. While we observed an association between response and TMB, this observation was confounded by disease subtype, strongly suggesting that TMB cannot be applied generically across melanoma subtypes as a predictive biomarker for anti-PD1 ICB.

Beyond TMB, we found that MHC-II expression, tumor heterogeneity, purity and ploidy were associated with ICB response. In two previous studies^47,48, MHC-II expression on tumor cells by immunohistochemistry was found to be predictive for response to anti-PD1 ICB and was hypothesized to represent a subset of tumors that could stimulate CD4⁺ helper T cell or cytotoxic activity. Consistent with this hypothesis, we found that MHC-II transcriptomic expression was correlated with expression of CD4 and the cytolytic molecules PRF1 and GZMA in our cohort. However, whether MHC-II expression represents expression on tumor cells or antigen-presenting cells within the tumor microenvironment cannot be determined from our bulk transcriptome data, and whether the association of PD1 ICB response with MHC-II expression is limited to tumor-cell-specific MHC-II expression is unclear. Notably, CD8⁺ T cell markers were not higher in responders versus progressors in our cohort, and although MHC-II, MHC-I, IFN-γ and IFN-α response pathway expression was correlated, MHC-II expression was the best predictor of response in our cohort. Further, we found evidence for the involvement of other immune compartments (for example, B cell markers enriched in ipilimumab-experienced responders) in ICB response, consistent with data from a recent trial demonstrating higher B cell infiltrate in responders to neoadjuvant immunotherapy⁴⁶, although the specific cell types, functional states and tumor immune interactions are not yet well characterized.

Tumor heterogeneity (the proportion of subclonal mutations) has previously been associated with poor prognosis across multiple tumor types and therapies^49,50,51. High heterogeneity suggests a highly mutagenic disease and a high degree of subclonality, with a higher likelihood of preexisting or rapidly evolving resistant clones. Interestingly, in our cohort, four patients had a very high TMB with an alkylating chemotherapy mutational signature (and known previous alkylating chemotherapy); low tumor heterogeneity distinguished the two responders from the two nonresponders (with SD and PD as the best response), who had high tumor heterogeneity with a majority of subclonal mutations. This association has also been observed in dacarbazine-experienced patients with melanoma treated with anti-CTLA4 ICB⁴⁹, suggesting that tumor heterogeneity may be significantly correlated with ICB resistance.

In our cohort, higher ploidy and lower purity were associated with ICB response, but the biological basis of these relationships is unclear. Genome doubling events are common in cancer and may accelerate genome evolution by increasing the tolerance of genome instability⁵², and higher aneuploidy has been associated with a worse response to ICB^18,19. However, whether genome doubling (and higher ploidy) is also associated with increased immunogenicity is unclear. Purity is negatively correlated with expression of markers of immune response and may be a proxy for the level of immune response in the tumor microenvironment in this setting rather than an artifact of tumor sample processing. However, tumor purity may also reflect differences in tumor biology leading to intrinsic resistance.

Notably, we found that previous exposure to anti-CTLA4 ICB affected the predictors of response to anti-PD1 ICB, although patients with and without exposure had similar response rates to anti-PD1 ICB. Immune-related markers are strongly enriched in responders compared to progressors with previous ipilimumab exposure, but this relationship is less clear in ipilimumab-naive patients. Specifically, post-ipilimumab tumors with poor immune response at progression were resistant to further anti-PD1 ICB, whereas ipilimumab-naive immune-poor and immune-infiltrated tumors were similarly likely to respond to anti-PD1 ICB. Whether anti-CTLA4 ICB induces or, alternatively, reveals an immune-resistant state in a subset of melanomas is an important question that deserves further evaluation. Further, cross-resistance to sequential ICB may also predict resistance to simultaneous combination ICB; this hypothesis should be evaluated in an appropriate cohort.

Building on these findings, we constructed predictive models integrating clinical, genomic and transcriptomic characteristics to identify patients with melanoma with intrinsic resistance to anti-PD1 ICB. Integrating multiple clinical and molecular features resulted in superior discrimination compared to models with any single feature or modality. In patients treated with ipilimumab, low MHC-II expression and high LDH predicted intrinsic resistance, whereas lymph node metastasis predicted improved response to therapy. MHC-II^47,48 and LDH⁵³ have previously been implicated in predicting anti-PD1 responsiveness. Lymph node metastases might provide a reservoir of tumor-specific immune cells, facilitating their activation by physiologic lymph node function; recent experimental data in a murine model suggests that lymph node metastases are necessary for PD1 response⁵⁴, and recent clinical data showing greater tumor-resident T cell clone response to neoadjuvant compared to adjuvant immunotherapy further supports this hypothesis⁵⁵. Similarly, integrating tumor heterogeneity, ploidy and purity for ipilimumab-naive disease resulted in a higher AUC than was obtained with any single-feature model alone. Beyond predicting response, these parsimonious models strongly stratified patients by PFS and OS, suggesting potential clinical applicability in identifying patients at high and low risk.

These findings will require validation in independent and larger cohorts; at the time of our study, limited data were publicly available where molecularly sequenced tumors with previous treatment data and all relevant clinical parameters were available for validation. Further, heterogeneity in sequencing approaches and data normalization between cohorts hindered our ability to develop standardized features to create and validate models. However, our results highlight the value of integrating rich clinical data with molecular tumor characterization and the need to generate such multimodal data.

Methods

Patient cohort and clinical end points

Patients were identified in databases of participating sites. For enrollment, patients were required to have advanced melanoma and to have received PD1 blockade as a palliative treatment. Tissue obtained before PD1 blockade was required for enrollment and was collected during routine medical care. Clinicopathological and demographic data were collected from patient records locally and are shown in Table 1. Age, stage and ECOG performance status were documented before the first application of anti-PD1 ICB. LDH was measured within 28 d of the first application of nivolumab or pembrolizumab. OS was defined as the time between the first application of anti-PD1 ICB and the date of death (any cause). For subjects without documentation of death, OS was censored on the last date the patient was known to be alive. BOR to anti-PD1 ICB was assessed according to RECIST criteria v.1.1 by the participating sites. Patients achieving CR or PR as BOR were grouped as responders, whereas patients showing PD as the best response were referred to as progressors. Patients were classified as MR when achieving unequivocal responses in individual existing lesions but also progression in others or new lesions. PFS was defined as the time between the first application of anti-PD1 ICB and the date of documented disease progression. For patients without documentation of progression, PFS was censored on the last date the patient was known to be without progression.

This retrospective study and associated informed consent procedures were approved by the central Ethics Committee (EC) of the University Hospital Essen (12-5152-BO and 11-4715). Approval by the local EC was obtained by investigators if required by local regulations.

Samples

Samples were collected retrospectively and obtained by excision or biopsy of melanoma tissue, collected locally at the participating sites and provided formalin-fixed and paraffin-embedded (FFPE). Samples were collected between January 2013 and June 2016. The median time from biopsy to initiation of anti-PD1 blockade was 2.1 months with 90% of samples being collected 6 months before the first application of nivolumab or pembrolizumab. All biopsies were from metastatic sites, with the exception of eight biopsies; seven were from a primary lesion and one was from a recurrence at a primary site, representing less than 6% of the overall cohort.

Whole-exome and whole-transcriptome sequencing

DNA extraction, whole exome library preparation and sequencing were performed for samples as previously described^10,23. Slides were cut from FFPE blocks and macrodissected for tumor-enriched tissue. Paraffin was removed from FFPE sections and cores using CitriSolv (Fisher Scientific), followed by ethanol washes and tissue lysis overnight at 56 °C. Samples were then incubated at 90 °C to remove DNA cross links. Extraction of DNA, and, when possible, RNA was performed using the QIAGEN AllPrep DNA/RNA mini kit (51306). Germline DNA was obtained from peripheral blood mononuclear cells and adjacent normal tissue.

Whole-exome capture libraries were constructed from 100 ng of DNA from tumor and normal tissue after sample shearing, end repair and phosphorylation and ligation to barcoded sequencing adaptors. Ligated DNA was size selected for lengths of 200–350 bp and subjected to exonic hybrid capture using Illumina library preps. The sample was multiplexed and sequenced using Illumina HiSeq technology. The Illumina exome sequencing approach uses Illumina’s in-solution DNA-probe-based hybrid selection method that applies principles similar to those of Broad Institute–Agilent Technologies’ in-solution RNA-probe-based hybrid selection method^57,58 to generate Illumina exome sequencing libraries.

Total RNA was assessed for quality using the Caliper LabChip GX2. The percentage of fragments with a size greater than 200 nucleotides (DV200) was calculated using software. An aliquot of 200 ng of RNA was used as the input for first-strand cDNA synthesis using Illumina’s TruSeq RNA Access Library Prep kit. Synthesis of the second strand of cDNA was followed by indexed adaptor ligation. Subsequent PCR amplification enriched for adaptor-ligated fragments. The amplified libraries were quantified using an automated PicoGreen assay.

A total of 200 ng of each cDNA library, not including controls, was combined into four-plex pools. Capture probes that target the exome were added and hybridized for the recommended time. Following hybridization, streptavidin magnetic beads were used to capture the library-bound probes from the previous step. Two wash steps effectively removed any nonspecifically bound products. These same hybridization, capture and wash steps were repeated to assure high specificity. A second round of amplification enriched the captured libraries. After enrichment the libraries were quantified with qPCR using the KAPA Library Quantification kit for Illumina sequencing platforms and were then pooled at an equimolar ratio. The entire process was in a 96-well format and all pipetting was performed using Agilent Bravo or Hamilton Starlet.

Pooled libraries were normalized to 2 nM and denatured using 0.2 N NaOH before sequencing. Flowcell cluster amplification and sequencing were performed according to the manufacturer’s protocols using either the HiSeq 2000 v.3 or HiSeq 2500. Each run generated 76-bp paired-end reads with a dual eight-base index barcode. Data were analyzed using the Broad Picard Pipeline, which includes demultiplexing and data aggregation.

Quality control and variant calling

Initial exome sequence data processing and analysis were performed using pipelines at the Broad Institute. After alignment from the Broad Picard Pipeline, BAM files were uploaded into the Firehose infrastructure (https://software.broadinstitute.org/cancer/cga/Firehose), which managed intermediate analysis files executed by analysis pipelines. Sequencing data were incorporated into quality-control modules in Firehose to compare the tumor and normal genotypes and ensure concordance between samples. Quality-control cutoffs were as follows: mean target coverage >50× (tumor) and >20× (matched normal), cross-contamination of samples estimation (ContEst⁵⁹) <5% and tumor purity ≥10% (Extended Data Fig. 1).

Power calculation quality control

To limit our analysis to samples where we had adequate power to call somatic variants, we performed a downstream per-sample power calculation. For each sample, we performed a Monte Carlo simulation of 1,000 true clonal mutations using the following procedures:

Sample the number of reads from the sample-specific coverage distribution
Draw the number of tumor reads from a binomial distribution using the estimated tumor purity
Draw the number of mutation reads from a binomial distribution given the assumption of a heterozygous mutation and no copy number variation
Characterize the mutation as detected or not on the basis of a log odds threshold of 6.3 (consistent with the MuTect⁶⁰ implementation).

The estimated power to detect clonal mutations is the proportion of simulated mutations detected (for example, 800 detected out of 1,000 simulated clonal mutations is 80% power), which is a function of both sample-specific sequencing depth of coverage and tumor purity. Three tumors were excluded using this threshold (Extended Data Fig. 1).

Variant calling

The MuTect algorithm⁶⁰ was applied to identify somatic single-nucleotide variants in targeted exons, with computational filtering of artifacts introduced by DNA oxidation during sequencing⁶¹ or FFPE-based DNA extraction using a filter-based method. Strelka⁶² was applied to identify small insertions or deletions. Identified alterations were annotated using Oncotator⁶³.

TMB and neoantigen load

For purposes of analysis, the TMB was calculated as the log of the number of nonsynonymous mutations detected from WES. Mutations per Mb was calculated by dividing the total number of nonsynonymous mutations by the number of bases with sufficient coverage in the tumor and normal samples (≥14× and ≥8×, respectively) to call mutations (https://software.broadinstitute.org/cancer/cga/mutect). Neoantigen prediction was performed as previously described¹⁰. Briefly, HLA type was inferred using POLYSOLVER⁶⁴, which uses a Bayesian classifier to determine the genotype of each patient. Neoantigens were predicted for each patient by defining all novel nine and ten amino-acid sequences resulting from mutations and determining whether the predicted binding affinity to the patient’s germline HLA alleles was <500 nM using NetMHCpan (v.2.4)⁶⁵.

Copy number variants

The total number of copy number alterations for individual tumors was inferred using adaptations of a binary segmentation algorithm⁶⁶ (CapSeg) comparing fractional exon coverage for tumor segments to a panel of normal samples, generating exomic segments and segment copy number. Copy number data were inspected visually and manually for focal amplifications and deletions and genes were annotated with Oncotator⁶³. For allelic copy number alterations, heterozygous single-nucleotide polymorphisms were identified and integrated into the binary segmentation algorithm (Allelic CapSeg) and allelic segments were further adjusted for tumor purity and ploidy using estimates derived from ABSOLUTE⁶⁷. We then called allelic amplifications and deletions, following previously described methodology and criteria⁶⁸ integrating segment focality and the purity- and ploidy-corrected allelic copy number.

Purity and ploidy

Purity and ploidy were estimated using the ABSOLUTE algorithm⁶⁷, which integrates variant allele frequency distributions and copy number variants to estimate absolute tumor purity and ploidy and infer the cancer cell fraction, the proportion of cancer cells in the sample that contain each mutation. Allelic segments following purity and ploidy correction were used to estimate allelic copy number.

Heterogeneity and aneuploidy

Heterogeneity was estimated as the proportion of mutations in each sample that were inferred to be subclonal. Clonal mutations were defined as having a cancer cell fraction ≥0.8, while other mutations were defined as subclonal; we chose this definition as a simple conservative approach with high specificity.

To estimate aneuploidy, we used the proportion of the genome inferred to have an allelic amplification or deletion, using the allelic segmentation described above.

Mutational signatures

De novo mutational signatures were generated in this cohort using an adaptation of non-negative matrix factorization⁶⁹ via the Brunet update method⁷⁰, as previously described in detail⁵¹, with the R package SomaticSignatures⁷¹ and non-negative matrix factorization⁷². Cosine similarity was used to compare the discovered signatures to the 30 existing discovered and validated signatures in COSMIC^25,73, with a threshold of 0.85, and we also manually visualized and inspected similarities in mutational motifs between our signatures and COSMIC signatures.

Transcriptomic analysis

Whole-transcriptome sequencing data from FFPE tissues were aligned using STAR⁷⁴ and quantified with RSEM⁷⁵ to yield gene-level expression in transcripts per million (TPM). For RNA-seq quality control, sequencing- and alignment-specific metrics were considered for each sample. The following alignment metrics (output by the STAR alignment method) were considered: percentage of uniquely mapped reads, average mapped read length, number of splices, mismatch rate per base, percentage of multimapped reads, percentage of reads mapped to too many locations, percentage of unmapped reads due to too many mismatches, percentage of unmapped reads due to reads being too short and percentage of unmapped reads due to other reasons. Additionally, we considered the raw number of reads, the average read length, the read duplication rate and DV200 for each sample. Samples were clustered across quality-control metrics using principal-component analysis, and outlier samples were manually evaluated and discarded. Three samples were removed owing to poor quality: patient 143 was excluded owing to an abnormally low absolute number of reads (number of reads <1 million); patient 107 was excluded owing to an abnormally high percentage of reads mapped to too many locations (>10% of reads), likely indicating high numbers of short or degraded reads; and patient 61 was excluded owing to multiple aberrant quality-control metrics resulting in overall poor quality when considering all metrics in aggregate as well as an aberrant expression profile compared to all other samples. Only transcriptomes from tumors whose WES also passed quality control were included; the final patient cohort for RNA-seq analysis included n = 121 transcriptomes.

We excluded certain classes of noncoding genes that constituted a large (>10%) proportion of TPM in the majority of samples. Specifically we excluded genes characterized as snoRNA, using biomaRt⁷⁶ to download Ensembl biotype annotations (using the dataset ‘hsapiens_gene_ensembl’ and version ‘GRCh38.p10’) and excluding genes whose biotype was ‘snoRNA’ (n = 380 genes). We then regenerated a new TPM metric for each sample to normalize the total transcriptome sum to 1 million.

For analysis, only genes with TPM > 0 in 25% or more of the samples were included. This excluded 6,158 genes, with 20,848 genes passing this threshold.

GSEA and ssGSEA

GSEA³² was performed using the Cancer Hallmarks gene sets³³ from MSigDB at https://cloud.genepattern.org/. We used default settings with 10,000 gene set permutations to generate P and q values, and we compared progressors and responders in the overall cohort, the ipilimumab-treated subgroup and the ipilimumab-naive subgroup separately.

To generate nonparametric gene set scores in individual samples, we generated ssGSEA projections⁴³ for gene sets using rank normalization, including the MHC-II HLA genes (HLA-DMA, HLA-DMB, HLA-DOA, HLA-DOB, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1 and HLA-DRB5) and MHC-I HLA genes (HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, TAP1, TAP2 and B2M). ssGSEA scores were also generated for Cancer Hallmarks gene sets.

Absolute immune infiltrate and immune subsets

Estimation of the total immune infiltrate in each sample and immune cell subsets was performed using CIBERSORT with the LM22 gene set⁴¹ on the CIBERSORT website (http://cibersort.stanford.edu). The absolute mode was enabled and quantile normalization was disabled using the RNA-seq TPM matrix for the cohort.

A separate immune infiltrate and immune cell subsets score analysis was performed using single-cell-derived signatures, and the methodology used for generating normalized signature scores was as described by Jerby-Arnon and colleagues⁴².

As there is no gold standard for inferring immune infiltrate from bulk RNA-seq data, we chose to use CIBERSORT-inferred values throughout the analysis to allow comparisons with other studies, as CIBERSORT has been widely used (for example, in a pan-TCGA immune landscape analysis⁷⁷).

Gene expression signatures

Published gene expression signatures related to immune-checkpoint response were collected from the literature and validated in our cohort^{18,19,21,34,35,36,37,38,39,40}. Sample-wise scores for these gene signatures were calculated from RNA-seq data using TPM values and following the methodology described in corresponding studies. Genes with unavailable expression data were excluded from calculations of gene signature scores. In two gene signatures^39,40, genes were incorporated independently (that is, weighted) into the published model, but neither the direction nor the coefficient was available, so these signatures were excluded from evaluation. Differences in immune signature scores between responders (CR and PR) and progressors (PD) across all samples, in the ipilimumab-naive subset and in the ipilimumab-treated subset were tested using the Mann–Whitney U-test. The predictive utility of these immune signatures was evaluated with AUC values derived from ROC curves of gene signature scores in the complete cohort, ipilimumab-naive subset and the ipilimumab-treated subset. Results and detailed descriptions of evaluated gene signatures are provided in Supplementary Table 6.

Analysis

Two primary response comparisons were made: (1) responders (defined as having CR or PR as the best RECIST response) versus progressors (defined as having PD as the best RECIST response) and (2) progressors (defined as having non-PD as the best RECIST response) versus non-progressors.

Statistical tests were performed utilizing the Python scipy.stats package. To compare numeric features between response categories, including transcriptome expression, a nonparametric MWW rank-sum test (mannwhitneyu() function) was used to minimize the effects of outliers. For comparison of the proportion between response categories, a chi-squared test (chi2_contingency() function) was utilized. For association of binary variables (for example, association of gene alteration with responders versus progressors), a Fisher’s exact test (fisher_exact() function) was utilized to generate a P value. A conservative adjusted OR was generated by repeating the Fisher’s exact test, adding one to both the number of gene-mutant responders and progressors. All tests were two sided unless otherwise indicated.

Survival analyses were performed utilizing the Python lifelines package⁷⁸. For Kaplan–Meier curve survival analysis, a log-rank test (logrank_test() function) was used to compare survival curves.

Hierarchical clustering was performed using the clustermap() function from the Python seaborn package⁷⁹, with default settings including a Euclidean distance metric and the ‘single’ method of calculating cluster distance (minimization of the nearest point between clusters).

Validation

For validation, we reviewed the literature and found three studies^18,20,44 of advanced melanoma treated with anti-PD1 ICB with response, WES and RNA-seq data. However, one did not have information on previous ipilimumab treatment²⁰, and another¹⁸ had only two patients who were naive to ipilimumab and nine who were treated with ipilimumab with post-ipilimumab tumor biopsies and available WES and NanoString data; thus, we used the remaining cohort⁴⁴ as our primary validation cohort.

To allow appropriate validation, only cutaneous, occult, acral and mucosal samples were included from validation cohorts; specifically, uveal and ocular melanomas were excluded (Riaz cohort, n = 5 excluded). Only patients with evaluable response criteria were included (Riaz cohort, n = 2 excluded). WES, transcriptomic and heterogeneity data were obtained from https://github.com/riazn/bms038_analysis. Fragments per kilobase of transcript per million mapped reads values were converted to TPM to be consistent with our cohort normalization.

Predictive model generation and cross-validation

We used logistic regression for our model to predict PD as the best RECIST response versus non-PD rather than responder versus progressor to better reflect the real-world setting where all outcomes (PD, SD, MR, PR and CR) are possible. We evaluated genomic, transcriptomic and clinical features. Categorical features were converted to binary features for each categorical value. To be conservative, no gene-level mutations or expression values were individually considered. Global genomic tumor characteristics such as TMB, purity, ploidy, heterogeneity and aneuploidy were considered. Features were generated from the transcriptome, including ssGSEA values for gene sets representing Cancer Hallmarks pathways, and MHC-II and MHC-I antigen-presentation genes, as well as gene expression signatures following the methodology in the respective publications, as described above and in Supplementary Table 6. Clinical characteristics including LDH and ECOG at the start of anti-PD1 ICB, the number of metastatic organs, sex, M stage, the number of different metastatic sites, metastatic sites and melanoma subtype were evaluated (Supplementary Table 1). Features were chosen in a forward-selection-based process, where features that were significantly predictive (P < 0.05) when added to the base model were ranked on the basis of the ability of the combined model to discriminate outcomes (using ROC curve AUC as the metric), and the best feature was chosen to be added to the base model. Potential features were evaluated on the basis of a manual review considering biological interpretability and clinical applicability. This process was iterated with the new base model and stopped when no features under consideration were significantly predictive.

The set of tumors with both WES and RNA-seq data was smaller than the set of tumors with only WES data; when the features chosen in model development for ipilimumab-naive tumors resulted in only WES features being chosen, model development was repeated in the superset of tumors requiring only clinical and WES data, and this model in the larger set is reported in the main text.

To estimate the ‘out-of-bag’ AUC, we used k-fold cross-validation (splitting the dataset into k-subsets, training on k − 1 subsets and calculating AUC on the holdout subset) and calculated the mean cross-validation AUC. Given the partially manual review of features, feature selection was not included in cross-validation. For the ipilimumab-treated subset (n = 34), we chose k = 5 folds, and for the larger ipilimumab-naive subset (n = 85), we chose k = 10 folds, to maintain a cross-validation holdout set of >5 tumors. Cross-validation scores were calculated using the cross_val_score function from the Python sklearn package.

To further evaluate the statistical support for our models, we calculated the Akaike information criteria and BIC of each subsequent model after adding an additional feature in forward selection in the ipilimumab-experienced and ipilimumab-naive subgroups (Extended Data Fig. 8c,d), and we also evaluated the addition of TMB as an additional feature to the selected models.

Permutation testing

To test the statistical significance of differences in FDR q values of enriched pathways between ipilimumab-naive versus ipilimumab-experienced subgroups, we performed a permutation analysis. Briefly, we shuffled the ipilimumab-experienced and ipilimumab-naive labels for each tumor, keeping each subgroup size the same and keeping the same number of responders and progressors in each subgroup, and we reran GSEA on each new simulated subgroup. We repeated this 1,000 times to generate a distribution of enriched pathways in each subgroup under the null hypothesis of no relationship between subgroup and pathway enrichment. We then compared our observed outcome within this null distribution to generate an empiric P value. For example, for each pathway enriched in ipilimumab-experienced patients, the proportion of simulations with a difference in log q value between ipilimumab-experienced and ipilimumab-naive subgroups (equivalent to the ratio of q values) greater than or equal to our observed difference, and with an FDR q value in the ipilimumab-experienced subgroup equal to or more extreme than our observed ipilimumab-experienced q value would represent the empiric P value.

We performed a similar permutation test to generate an empiric P value for the predictive models for ipilimumab-experienced and ipilimumab-naive subgroups. Briefly, we permuted the outcome labels (progressors and non-progressors) within each subgroup, and generated an AUC and cross-validation AUC for the predictive model with the specified features (that is, MHC-II, LDH and lymph node metastasis for ipilimumab-experienced tumors; purity, ploidy and heterogeneity for ipilimumab-naive tumors) to generate a null distribution of AUCs and cross-validation AUCs under the null hypothesis that these predictors are not associated with outcomes. By permuting the phenotype rather than the predictors, we preserved the inter-predictor structure. Then, the proportion of simulations with AUC and cross-validation AUC greater than our observed AUC represented an empiric P value.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability

All reasonable requests for raw and analyzed data and materials will be promptly reviewed by the senior authors to determine whether the request is subject to any intellectual property or confidentiality obligations. Patient-related data not included in the paper may be subject to patient confidentiality. Any data and materials that can be shared will be released via a material transfer agreement. All analyzed sequencing data are in supplementary tables or data available online. Raw sequencing data are available in dbGaP (accession number phs000452.v3.p1).

Code availability

Python code is available in packages as described in the manuscript. Code to regenerate figures using supplementary data and tables is available at GitHub at https://github.com/vanallenlab/schadendorf-pd1. Additional reasonable requests for code will be promptly reviewed by the senior authors to verify whether the request is subject to any intellectual property or confidentiality obligations, and shared to the extent permissible by these obligations.

Change history

19 June 2020
An amendment to this paper has been published and can be accessed via a link at the top of the paper.

References

Hodi, F. S. et al. Improved survival with ipilimumab in patients with metastatic melanoma. N. Engl. J. Med. 363, 711–723 (2010).
Article CAS PubMed PubMed Central Google Scholar
Robert, C. et al. Pembrolizumab versus ipilimumab in advanced melanoma. N. Engl. J. Med. 372, 2521–2532 (2015).
Article CAS PubMed Google Scholar
Wolchok, J. D. et al. Overall survival with combined nivolumab and ipilimumab in advanced melanoma. N. Engl. J. Med. 377, 1345–1356 (2017).
Article CAS PubMed PubMed Central Google Scholar
Reck, M. et al. Pembrolizumab versus chemotherapy for PD-L1-positive non-small cell lung cancer. N. Engl. J. Med. 375, 1823–1833 (2016).
Article CAS PubMed Google Scholar
Rosenberg, J. E. et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet 387, 1909–1920 (2016).
Article CAS PubMed PubMed Central Google Scholar
Motzer, R. J. et al. Nivolumab versus everolimus in advanced renal-cell carcinoma. N. Engl. J. Med. 373, 1803–1813 (2015).
Article CAS PubMed PubMed Central Google Scholar
Seiwert, T. Y. et al. Safety and clinical activity of pembrolizumab for treatment of recurrent or metastatic squamous cell carcinoma of the head and neck (KEYNOTE-012): an open-label, multicentre, phase 1b trial. Lancet Oncol. 17, 956–965 (2016).
Article CAS PubMed Google Scholar
Rizvi, N. A. et al. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science 348, 124 (2015).
Article CAS PubMed PubMed Central Google Scholar
Snyder, A. et al. Genetic basis for clinical response to CTLA-4 blockade in melanoma. N. Engl. J. Med. 371, 2189–2199 (2014).
Article PubMed PubMed Central CAS Google Scholar
Van Allen, E. M. et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science 350, 207–211 (2015).
Article PubMed PubMed Central CAS Google Scholar
Cristescu, R. et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. Science 362, eaar3593 (2018).
Article PubMed PubMed Central CAS Google Scholar
Tumeh, P. C. et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature 515, 568–571 (2014).
Article CAS PubMed PubMed Central Google Scholar
Zaretsky, J. M. et al. Mutations associated with acquired resistance to PD-1 blockade in melanoma. N. Engl. J. Med. 375, 819–829 (2016).
Article CAS PubMed PubMed Central Google Scholar
Sade-Feldman, M. et al. Resistance to checkpoint blockade therapy through inactivation of antigen presentation. Nat. Commun. 8, 1136 (2017).
Article PubMed PubMed Central CAS Google Scholar
Gao, J. et al. Loss of IFN-γ pathway genes in tumor cells as a mechanism of resistance to anti-CTLA-4 therapy. Cell 167, 397–404 (2016).
Article CAS PubMed PubMed Central Google Scholar
Peng, W. et al. Loss of PTEN promotes resistance to T cell-mediated immunotherapy. Cancer Discov. 6, 202–216 (2016).
Article CAS PubMed Google Scholar
George, S. et al. Loss of PTEN is associated with resistance to anti-PD-1 checkpoint blockade therapy in metastatic uterine leiomyosarcoma. Immunity 46, 197–204 (2017).
Article CAS PubMed PubMed Central Google Scholar
Roh, W. et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. Transl. Med. 9, eaah3560 (2017).
Article PubMed PubMed Central CAS Google Scholar
Davoli, T. et al. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science 355, eaaf8399 (2017).
Article PubMed PubMed Central CAS Google Scholar
Hugo, W. et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell 165, 35–44 (2016).
Article CAS PubMed PubMed Central Google Scholar
Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G. & Hacohen, N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61 (2015).
Article CAS PubMed PubMed Central Google Scholar
Weber, J. S. et al. Sequential administration of nivolumab and ipilimumab with a planned switch in patients with advanced melanoma (CheckMate 064): an open-label, randomised, phase 2 trial. Lancet Oncol. 17, 943–955 (2016).
Article CAS PubMed PubMed Central Google Scholar
Van Allen, E. M. et al. Whole-exome sequencing and clinical interpretation of formalin-fixed, paraffin-embedded tumor samples to guide precision cancer medicine. Nat. Med. 20, 682–688 (2014).
Article PubMed PubMed Central CAS Google Scholar
Birkeland, E. et al. Patterns of genomic evolution in advanced melanoma. Nat. Commun. 9, 2665 (2018).
Article CAS PubMed PubMed Central Google Scholar
Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).
Article CAS PubMed PubMed Central Google Scholar
Hayward, N. K. et al. Whole-genome landscapes of major melanoma subtypes. Nature 545, 175–180 (2017).
Article CAS PubMed Google Scholar
Miao, D. et al. Genomic correlates of response to immune checkpoint blockade in microsatellite-stable solid tumors. Nat. Genet. 50, 1271–1281 (2018).
Article CAS PubMed PubMed Central Google Scholar
McGranahan, N. et al. Allele-specific HLA loss and immune escape in lung cancer evolution. Cell 171, 1259–1271 (2017).
Article CAS PubMed PubMed Central Google Scholar
Riaz, N. et al. Recurrent SERPINB3 and SERPINB4 mutations in patients that respond to anti-CTLA4 immunotherapy. Nat. Genet. 48, 1327–1329 (2016).
Article CAS PubMed PubMed Central Google Scholar
Shin, D. S. et al. Primary resistance to PD-1 blockade mediated by JAK1/2 mutations. Cancer Discov. 7, 188–201 (2017).
Article CAS PubMed Google Scholar
Sucker, A. et al. Acquired IFN-γ resistance impairs anti-tumor immunity and gives rise to T cell-resistant melanoma lesions. Nat. Commun. 8, 15440 (2017).
Article CAS PubMed PubMed Central Google Scholar
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).
Article CAS PubMed PubMed Central Google Scholar
Liberzon, A. et al. The molecular signatures database hallmark gene-set collection. Cell Syst. 1, 417–425 (2015).
Article CAS PubMed PubMed Central Google Scholar
Auslander, N. et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat. Med. 24, 1545–1549 (2018).
Article CAS PubMed PubMed Central Google Scholar
Fehrenbacher, L. et al. Atezolizumab versus docetaxel for patients with previously treated non-small cell lung cancer (POPLAR): a multicentre, open-label, phase 2 randomised controlled trial. Lancet 387, 1837–1846 (2016).
Article CAS PubMed Google Scholar
Huang, A. C. et al. A single dose of neoadjuvant PD-1 blockade predicts clinical outcomes in resectable melanoma. Nat. Med. 25, 454–461 (2019).
Article CAS PubMed PubMed Central Google Scholar
Jiang, P. et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558 (2018).
Article CAS PubMed PubMed Central Google Scholar
Messina, J. L. et al. 12-Chemokine gene signature identifies lymph node-like structures in melanoma: potential for patient selection for immunotherapy? Sci. Rep. 2, 765 (2012).
Article PubMed PubMed Central CAS Google Scholar
Ayers, M. et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J. Clin. Invest. 127, 2930–2940 (2017).
Article PubMed PubMed Central Google Scholar
Ock, C.-Y. et al. Genomic landscape associated with potential response to anti-CTLA-4 treatment in cancers. Nat. Commun. 8, 1050 (2017).
Article PubMed PubMed Central CAS Google Scholar
Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457 (2015).
Article CAS PubMed PubMed Central Google Scholar
Jerby-Arnon, L. et al. A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade. Cell 175, 984–997 (2018).
Article CAS PubMed PubMed Central Google Scholar
Barbie, D. A. et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112 (2009).
Article CAS PubMed PubMed Central Google Scholar
Riaz, N. et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell 171, 934–949 (2017).
Article CAS PubMed PubMed Central Google Scholar
Weber, J. et al. Baseline tumor T cell receptor (TcR) sequencing analysis and neo antigen load is associated with benefit in melanoma patients receiving sequential nivolumab and ipilimumab. Ann. Oncol. 27, 10470 (2016).
Article Google Scholar
Amaria, R. N. et al. Neoadjuvant immune checkpoint blockade in high-risk resectable melanoma. Nat. Med. 24, 1649–1654 (2018).
Article CAS PubMed PubMed Central Google Scholar
Rodig, S. J. et al. MHC proteins confer differential sensitivity to CTLA-4 and PD-1 blockade in untreated metastatic melanoma. Sci. Transl. Med. 10, eaar3342 (2018).
Article PubMed CAS Google Scholar
Johnson, D. B. et al. Melanoma-specific MHC-II expression represents a tumour-autonomous phenotype and predicts response to anti-PD-1/PD-L1 therapy. Nat. Commun. 7, 10582 (2016).
Article CAS PubMed PubMed Central Google Scholar
McGranahan, N. et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science 351, 1463–1469 (2016).
Article CAS PubMed PubMed Central Google Scholar
Morris, L. G. T. et al. Pan-cancer analysis of intratumor heterogeneity as a prognostic determinant of survival. Oncotarget 7, 10051–10063 (2016).
Article PubMed PubMed Central Google Scholar
Liu, D. et al. Mutational patterns in chemotherapy resistant muscle-invasive bladder cancer. Nat. Commun. 8, 2193 (2017).
Article CAS PubMed PubMed Central Google Scholar
Dewhurst, S. M. et al. Tolerance of whole-genome doubling propagates chromosomal instability and accelerates cancer genome evolution. Cancer Discov. 4, 175–185 (2014).
Article CAS PubMed PubMed Central Google Scholar
Weide, B. et al. Baseline biomarkers for outcome of melanoma patients treated with pembrolizumab. Clin. Cancer Res. 22, 5487–5496 (2016).
Article CAS PubMed PubMed Central Google Scholar
Fransen, M. F. et al. Tumor-draining lymph nodes are pivotal in PD-1/PD-L1 checkpoint therapy. JCI Insight 3, 124507 (2018).
Article PubMed Google Scholar
Blank, C. U. et al. Neoadjuvant versus adjuvant ipilimumab plus nivolumab in macroscopic stage III melanoma. Nat. Med. 24, 1655–1661 (2018).
Article CAS PubMed Google Scholar
Cesano, A. nCounter PanCancer immune profiling panel (NanoString Technologies, Seattle, WA). J. Immunother. Cancer 3, 42 (2015).
Article PubMed PubMed Central Google Scholar
Fisher, S. et al. A scalable, fully automated process for construction of sequence-ready human exome targeted capture libraries. Genome Biol. 12, R1 (2011).
Article PubMed PubMed Central Google Scholar
Gnirke, A. et al. Solution hybrid selection with ultra-long oligonucleotides for massively parallel targeted sequencing. Nat. Biotechnol. 27, 182–189 (2009).
Article CAS PubMed PubMed Central Google Scholar
Cibulskis, K. et al. ContEst: estimating cross-contamination of human samples in next generation sequencing data. Bioinformatics 27, 2601–2602 (2011).
Article CAS PubMed PubMed Central Google Scholar
Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219 (2013).
Article CAS PubMed PubMed Central Google Scholar
Costello, M. et al. Discovery and characterization of artifactual mutations in deep coverage targeted capture sequencing data due to oxidative DNA damage during sample preparation. Nucleic Acids Res. 41, e67 (2013).
Article CAS PubMed PubMed Central Google Scholar
Saunders, C. T. et al. Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs. Bioinformatics 28, 1811–1817 (2012).
Article CAS PubMed Google Scholar
Ramos, A. H. et al. Oncotator: cancer variant annotation tool. Hum. Mutat. 36, E2423–E2429 (2015).
Article PubMed PubMed Central Google Scholar
Shukla, S. A. et al. Comprehensive analysis of cancer-associated somatic mutations in class I HLA genes. Nat. Biotechnol. 33, 1152 (2015).
Article CAS PubMed PubMed Central Google Scholar
Nielsen, M. et al. NetMHCpan, a method for quantitative predictions of peptide binding to any HLA-A and HLA-B locus protein of known sequence. PloS ONE 2, e796 (2007).
Article PubMed PubMed Central CAS Google Scholar
Olshen, A. B., Venkatraman, E. S., Lucito, R. & Wigler, M. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics 5, 557–572 (2004).
Article PubMed Google Scholar
Carter, S. L. et al. Absolute quantification of somatic DNA alterations in human cancer. Nat. Biotechnol. 30, 413–421 (2012).
Article CAS PubMed PubMed Central Google Scholar
Brastianos, P. K. et al. Genomic characterization of brain metastases reveals branched evolution and potential therapeutic targets. Cancer Discov. 5, 1164–1177 (2015).
Article CAS PubMed PubMed Central Google Scholar
Lee, D. D. & Seung, H. S. Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791 (1999).
Article CAS PubMed Google Scholar
Brunet, J. P. et al. Metagenes and molecular pattern discovery using matrix factorization. Proc. Natl Acad. Sci. USA 101, 4164–4169 (2004).
Article CAS PubMed PubMed Central Google Scholar
Gehring, J. S., Fischer, B., Lawrence, M. & Huber, W. SomaticSignatures: inferring mutational signatures from single-nucleotide variants. Bioinformatics 31, 3673–3675 (2015).
CAS PubMed PubMed Central Google Scholar
Gaujoux, R. & Seoighe, C. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics 11, 367 (2010).
Article PubMed PubMed Central CAS Google Scholar
Forbes, S. A. et al. COSMIC: exploring the world’s knowledge of somatic mutations in human cancer. Nucleic Acids Res. 43, D805–D811 (2015).
Article CAS PubMed Google Scholar
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Article CAS PubMed Google Scholar
Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011).
Article CAS PubMed PubMed Central Google Scholar
Durinck, S. et al. BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics 21, 3439–3440 (2005).
Article CAS PubMed Google Scholar
Thorsson, V. et al. The immune landscape of cancer. Immunity 48, 812–830 (2018).
Article CAS PubMed PubMed Central Google Scholar
Davidson-Pilon C. et al. CamDavidsonPilon/lifelines v.0.18.5. Zenodo https://doi.org/10.5281/zenodo.2562267 (2019).
Waksom M. et al. mwaskom/seaborn v.0.8.1 Zenodo https://doi.org/10.5281/zenodo.883859 (2017).

Download references

Acknowledgements

The authors thank the active investigators A. Gesierich (University Hospital Würzburg), J.C. Hassel (NCT Heidelberg), C. Pföhler (Saarland University Medical Center), E. Dabrowski (Ludwigshafen Medical Center), L.A. Schneider (University Medical Center Ulm), C. Weishaupt (University Hospital of Münster), K.G. Griewank (University Hospital Essen), E. Hadaschik (University Hospital Essen), G. Kyriakakis (Laikon General Hospital), F. Meier (NCT Dresden), M.H. Geukes Foppen (The Netherlands Cancer Institute), R. Dummer (University Hospital Zürich), E. Bräunlein (Technical University Munich) and M. Boxberg (Technical University Munich) for providing additional patient samples and clinical data. The authors also thank A. Giobbie-Hurder (Dana-Farber Cancer Institute) and J. Weirather (Dana-Farber Cancer Institute) for useful discussions regarding statistical testing. This work was supported by the Adelson Medical Research Foundation (L.G.), the Conquer Cancer Foundation (David Liu), the Society for Immunotherapy of Cancers (David Liu), the Damon Runyon Cancer Research Foundation (David Liu), the BroadNext10 (E.M.V.), the National Institutes of Health (K08 CA234458 (David Liu), R01 CA227388 (E.M.V.), U01 CA233100 (E.M.V.), T32 GM008313 (M.X.H.)), the Deutsche Forschungsgemeinschaft (German Research Foundation), SCHA 422/17-1, PA 2376/1-1 and HO 6389/2-1 (KFO 337; D.S, A.P.) and National Science Foundation Graduate Research Fellowship Program DGE1144152 (M.X.H.).

Author information

These authors contributed equally: David Liu, Bastian Schilling.

Authors and Affiliations

Dana-Farber Cancer Institute, Boston, MA, USA
David Liu, Derek Liu, Natalie Vokes, Claire A. Margolis, Haitham Elmarakeby, Felix Dietlein, Diana Miao, Ben Izar & Eliezer M. Van Allen
Broad Institute of Harvard and MIT, Cambridge, MA, USA
David Liu, Derek Liu, Livnat Jerby-Arnon, Natalie Vokes, Claire A. Margolis, Jake Conway, Meng Xiao He, Haitham Elmarakeby, Felix Dietlein, Diana Miao, Adam Tracy, Ben Izar, Aviv Regev & Eliezer M. Van Allen
Department of Dermatology, University Hospital Würzburg, Würzburg, Germany
Bastian Schilling
Department of Dermatology, University Hospital, Essen, Germany
Bastian Schilling, Antje Sucker, Elisabeth Livingstone, Lisa Zimmer, Annette Paschen & Dirk Schadendorf
German Cancer Consortium of Translational Cancer Research, German Cancer Research Center, Heidelberg, Germany
Bastian Schilling, Antje Sucker, Elisabeth Livingstone, Lisa Zimmer, Angela Krackhardt, Annette Paschen & Dirk Schadendorf
Harvard Medical School, Boston, MA, USA
Derek Liu, Jake Conway, Meng Xiao He & Diana Miao
Skin Cancer Center Hannover, Department of Dermatology and Allergy, Hannover Medical School, Hannover, Germany
Ralf Gutzmer & Imke Satzger
Department of Dermatology, University Medical Center, Mainz, Germany
Carmen Loquai & Stephan Grabbe
Biophysics Program, Harvard University, Cambridge, MA, USA
Meng Xiao He
First Department of Medicine, National and Kapodistrian University of Athens, Athens, Greece
Helen Gogas
Department of Dermatology, University Hospital Zürich, Zürich, Switzerland
Simone M. Goldinger
Department of Dermatology, University Medical Center Mannheim, Ruprecht-Karl University of Heidelberg, Mannheim, Germany
Jochen Utikal
Skin Cancer Unit, German Cancer Research Center, Heidelberg, Germany
Jochen Utikal
Department of Medical Oncology, The Netherlands Cancer Institute, Amsterdam, the Netherlands
Christian U. Blank
Skin Cancer Center at the University Cancer Centre, Department of Dermatology, Faculty of Medicine, University Hospital Carl Gustav Carus, Technische Universität Dresden, Dresden, Germany
Ricarda Rauschenberg
National Center for Tumor Diseases, Dresden, Germany
Ricarda Rauschenberg
German Cancer Research Centre, Heidelberg, Germany
Ricarda Rauschenberg
Department of Dermatology, Medical Center- University of Freiburg, Faculty of Medicine, University of Freiburg, Freiburg, Germany
Dagmar von Bubnoff
Medizinische Klinik III, Klinikum Rechts der Isar, Technische Universität München, Munich, Germany
Angela Krackhardt
Department of Dermatology, University Medical Center Tübingen, Tübingen, Germany
Benjamin Weide
Department of Dermatology, University Hospital Regensburg, Regensburg, Germany
Sebastian Haferkamp
Department of Dermatology, University Hospital Berlin, Berlin, Germany
Felix Kiecker
Eli Lilly and Co., Indianapolis, IN, USA
Levi Garraway
Massachusetts General Hospital, Boston, MA, USA
Keith Flaherty

Authors

David Liu
View author publications
You can also search for this author in PubMed Google Scholar
Bastian Schilling
View author publications
You can also search for this author in PubMed Google Scholar
Derek Liu
View author publications
You can also search for this author in PubMed Google Scholar
Antje Sucker
View author publications
You can also search for this author in PubMed Google Scholar
Elisabeth Livingstone
View author publications
You can also search for this author in PubMed Google Scholar
Livnat Jerby-Arnon
View author publications
You can also search for this author in PubMed Google Scholar
Lisa Zimmer
View author publications
You can also search for this author in PubMed Google Scholar
Ralf Gutzmer
View author publications
You can also search for this author in PubMed Google Scholar
Imke Satzger
View author publications
You can also search for this author in PubMed Google Scholar
Carmen Loquai
View author publications
You can also search for this author in PubMed Google Scholar
Stephan Grabbe
View author publications
You can also search for this author in PubMed Google Scholar
Natalie Vokes
View author publications
You can also search for this author in PubMed Google Scholar
Claire A. Margolis
View author publications
You can also search for this author in PubMed Google Scholar
Jake Conway
View author publications
You can also search for this author in PubMed Google Scholar
Meng Xiao He
View author publications
You can also search for this author in PubMed Google Scholar
Haitham Elmarakeby
View author publications
You can also search for this author in PubMed Google Scholar
Felix Dietlein
View author publications
You can also search for this author in PubMed Google Scholar
Diana Miao
View author publications
You can also search for this author in PubMed Google Scholar
Adam Tracy
View author publications
You can also search for this author in PubMed Google Scholar
Helen Gogas
View author publications
You can also search for this author in PubMed Google Scholar
Simone M. Goldinger
View author publications
You can also search for this author in PubMed Google Scholar
Jochen Utikal
View author publications
You can also search for this author in PubMed Google Scholar
Christian U. Blank
View author publications
You can also search for this author in PubMed Google Scholar
Ricarda Rauschenberg
View author publications
You can also search for this author in PubMed Google Scholar
Dagmar von Bubnoff
View author publications
You can also search for this author in PubMed Google Scholar
Angela Krackhardt
View author publications
You can also search for this author in PubMed Google Scholar
Benjamin Weide
View author publications
You can also search for this author in PubMed Google Scholar
Sebastian Haferkamp
View author publications
You can also search for this author in PubMed Google Scholar
Felix Kiecker
View author publications
You can also search for this author in PubMed Google Scholar
Ben Izar
View author publications
You can also search for this author in PubMed Google Scholar
Levi Garraway
View author publications
You can also search for this author in PubMed Google Scholar
Aviv Regev
View author publications
You can also search for this author in PubMed Google Scholar
Keith Flaherty
View author publications
You can also search for this author in PubMed Google Scholar
Annette Paschen
View author publications
You can also search for this author in PubMed Google Scholar
Eliezer M. Van Allen
View author publications
You can also search for this author in PubMed Google Scholar
Dirk Schadendorf
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

David Liu, B.S., L.G., A.R., K.F., E.M.V.A. and D.S. conceived and designed the overall study. B.S. collected and reviewed all clinical data. A.S. performed sample processing and shipping. B.S., E.L., L.Z., R.G., I.S., C.L., S.G., H.G., S.M.G., J.U., C.U.B, R.R., D.v.B., A.K., B.W., S.H. and F.K. provided samples and clinical annotation. A.T. oversaw sample processing and sequencing. David Liu, C.A.M., N.V., J.C., D.M., F.D. and M.X.H. designed and performed the mutational and copy number analyses. David Liu, Derek Liu, C.A.M., M.X.H., D.M., L.J.-A. and B.I. designed and performed RNA-seq analysis, including quality control, normalization, signature analysis, CIBERSORT and GSEA. David Liu, H.E. and M.X.H. designed and constructed the predictive modeling. David Liu, B.S., A.P., E.M.V.A. and D.S. interpreted the data. David Liu, B.S., E.M.V.A. and D.S. wrote the manuscript, and all authors reviewed and approved the final manuscript.

Corresponding authors

Correspondence to Eliezer M. Van Allen or Dirk Schadendorf.

Ethics declarations

Competing interests

David Liu reports funding by a postdoctoral fellowship from the Society for Immunotherapy of Cancers, which is funded in part by an educational grant from Bristol-Meyers Squibb (BMS). BMS has had no input into the conception, conduct or reporting of the submitted work. B.S. is on the advisory board or has received honoraria from Incyte, Novartis, Roche, BMS and Merck Sharp & Dohme (MSD), research funding from Pierre-Fabre Pharmaceuticals, BMS and MSD and travel support from Novartis, Roche, BMS, Pierre-Fabre Pharmaceuticals, MSD and Amgen, outside the scope of the submitted work. D.S. reports grants, personal fees and nonfinancial support from BMS, personal fees and nonfinancial support from Roche, grants, personal fees and nonfinancial support from Novartis, nonfinancial support from Regeneron, personal fees from Sanofi, personal fees and nonfinancial support from MSD, personal fees and nonfinancial support from Amgen, personal fees and nonfinancial support from 4SC, personal fees and nonfinancial support from Merck-EMD, personal fees from Array, personal fees and nonfinancial support from Pierre-Fabre, personal fees and nonfinancial support from Philogen, personal fees and nonfinancial support from Incyte and personal fees from Pfizer, outside the scope of the submitted work. E.M.V.A. reports advisory relationships and consulting with Tango Therapeutics, Genome Medical, Invitae, Illumina and Ervaxx; research support from Novartis and BMS; equity in Tango Therapeutics, Genome Medical, Syapse, Ervaxx and Microsoft; and travel reimbursement from Roche and Genentech, outside the submitted work. S.G. reports personal fees from MSD, personal fees from BMS, personal fees from Novartis and personal fees from Roche, outside the scope of the submitted work. L.Z. reports personal fees and others from BMS, personal fees and others from Novartis, personal fees and others from Pierre-Fabre, personal fees and others from MSD, personal fees from Roche, other fees from Amgen and personal fees from Sanofi, outside the scope of the submitted work. F.K. reports personal fees from Amgen, personal fees from BMS, grants and personal fees from Novartis, personal fees from Roche, personal fees from MSD and personal fees from Merck, outside the submitted work. J.U. reports personal fees and nonfinancial support from Amgen, personal fees and nonfinancial support from BMS, personal fees and nonfinancial support from MSD, personal fees and nonfinancial support from Novartis, personal fees and nonfinancial support from Pierre-Fabre, and personal fees and nonfinancial support from Roche, outside the scope of the submitted work. H.G. reports grants and personal fees from BMS, grants and personal fees from Roche, grants and personal fees from MSD, grants and personal fees from Novartis, personal fees from Amgen and personal fees from Pierre-Fabre, outside the scope of the submitted work. R.G. reports personal fees and nonfinancial support from BMS, personal fees and nonfinancial support from Roche Pharma, personal fees and nonfinancial support from Merck Serono, grants, personal fees and nonfinancial support from Amgen, personal fees and nonfinancial support from Pierre-Fabre, personal fees and nonfinancial support from Sanofi Regeneron, personal fees from MSD, grants, personal fees and nonfinancial support from Novartis, personal fees and nonfinancial support from Almirall Hermal, grants and personal fees from Pfizer, personal fees from LEO, personal fees from AstraZeneca, personal fees from Sun Pharma, personal fees from 4SC and grants from Johnson & Johnson, outside the scope of the submitted work. I.S. reports grants and personal fees from Novartis, grants from Pfizer, personal fees from Roche, personal fees BMS and personal fees from MSD, outside the scope of the submitted work. C.L. reports personal fees from Roche, personal fees from Novartis, personal fees from Pierre-Fabre, personal fees from BMS, personal fees from MSD, personal fees from Amgen and personal fees from Idera and Sun Pharma, outside the scope of the submitted work. S.G. reports personal fees from BMS, personal fees from MSD and personal fees from Merck KGaA, outside the scope of the submitted work. E.L. reports personal fees and others from Amgen, personal fees and others from BMS, personal fees and others from MSD, personal fees and others from Novartis, personal fees and others from Roche, personal fees from medac, personal fees from Janssen, other fees from Actelion and other fees from Pierre-Fabre, outside the scope of the submitted work. A.K. reports grants from BMS and Kiadis, personal fees from BMS, Sanofi, Novartis, Roche, Vaccibody and nonfinancial support from Sanofi and BMS, outside the scope of the submitted work. C.B. reports personal fees from BMS, MSD, Roche, Novartis, GSK, Pfizer, Lilly, Pierre-Fabre and GenMab and grants from BMS, Novartis and NanoString, outside the scope of the submitted work. S.H. reports personal fees from Novartis, personal fees from BMS, personal fees from Amgen, personal fees from Pierre-Fabre and personal fees from Roche, outside the scope of the submitted work. B.W. reports grants and personal fees from Philogen, personal fees from Curevac, grants and personal fees from MSD, grants and personal fees from BMS and personal fees from GSK, outside the scope of the submitted work. R.R. reports personal fees from Novartis, personal fees and nonfinancial support from Amgen and nonfinancial support from BMS, outside the scope of the submitted work. D.v.B., A.S. and A.P. have nothing to disclose.

Additional information

Peer review information Saheli Sadanand was the primary editor on this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data

Extended Data Fig. 1

Consort Diagram showing inclusion, exclusion, and quality control criteria for patients/tumors included in analysis.

Extended Data Fig. 2 Survival and genomic copy number characteristics of responders (n = 55 patients) (defined as CR or PR as best response) vs. progressors (n = 65 patients) (defined as PD as best response).

(a) PFS Kaplan Meier survival curves; two-sided log-rank p = 2.1e-28 (b) OS Kaplan Meier survival curves; two-sided log-rank p = 5.1e-18 (c) Proportion of tumor genome with copy number alterations (two-sided MWW p = 0.09). Boxplots: Box limits indicate the interquartile range (25th-75th percentile), with a center line indicating the median. Whiskers show the value ranges up to 1.5xIQR above the 75th or below the 25th percentile, with outliers beyond those ranges shown as individual points. PD = progressive disease. CR = complete response. PR = partial response.

Extended Data Fig. 3 Association of Genomic Alterations in Genes with Response vs. Progression.

For each gene, the association of a genomic alteration with response was tested using a two-sided Fisher’s Exact test, and the odds ratio (OR) and p-value calculated. The adjusted OR is a conservative estimate generated by adding 1 to both mutant gene responders and progressors, thus moving the OR estimate closer to 1. The x- and y-axis are on log scales, but labeled in the original units. (a) Association of gene mutations with response (n = 55 patients) vs. progression (n = 65 patients). (b) Association of gene amplifications with response (n = 55 patients) vs. progression (n = 65 patients).

Extended Data Fig. 4 Amplification of Chromosome 6 Regions in TAP2, HLA-A/B/C amplified tumors.

Patients with amplifications in TAP2 or HLA-A/B/C are shown, with inferred amplifications in the relevant section of Chromosome 6 shown. 6 patients had amplifications in TAP2, 6 patients had amplifications in HLA-A/B/C, and 4 patients had amplifications of a chromosomal region including all four genes.

Extended Data Fig. 5 Hierarchical clustering of the correlation matrix between transcriptional signatures previously associated with immunotherapy response and hallmark gene-sets, genomic, clinical, and transcriptomic features associated with response in our cohort.

The color indicates the Pearson correlation between features, from perfect negative correlation (Pearson r = -1, blue) to perfect positive correlation (Pearson r = 1, red). Most previously hypothesized signatures cluster within an immune-activity related group with immune-related gene-sets. Sample size for each correlation depends on the number of available data points: correlations involving exclusively genomic or clinical data had n = 144 tumor samples, while correlations involving transcriptomic features had n = 121 tumor samples with data available.

Extended Data Fig. 6 Expression of Genomic Features in Responders vs. Progressors in ipilimumab-treated and -naive subsets.

Two-sided Mann-Whitney-Wilcoxon tests were used to compare responders to progressors (ipilimumab-treated: n = 16 responders/20 progressors; ipilimumab-naïve: n = 34 responders/37 progressors). All p-values are unadjusted. Log Nonsyn Mutload: Ipi-treated p = 0.15; Ipi-naïve p = 0.15. Log Clonal Mutload: Ipi-treated p = 0.17; Ipi-naïve p = 0.08. Heterogeneity: Ipi-treated p = 0.51; Ipi-naïve p = 0.057. Ploidy: Ipi-treated p = 0.45; Ipi-naïve p = 0.002. Boxplots: Box limits indicate the interquartile range (25th-75th percentile), with a center line indicating the median. Whiskers show the value ranges up to 1.5xIQR above the 75th or below the 25th percentile, with outliers beyond those ranges shown as individual points.

Extended Data Fig. 7 Survival of patients stratified by progressive disease (PD) or non-PD (SD, MR, PR, or CR) best response (RECIST 1.1) to PD-1 ICB.

(a) PFS (two-sided logrank p = 4.0e-38 (b) OS (two-sided logrank p = 1.2e-19).

Extended Data Fig. 8 Model selection criteria with increasing number of features by forward selection; cross-validation AUC distribution of final models.

(a) Model AUC for the ipi-treated subgroup. Each point is labeled with the additional feature, and the final model includes MHC-II, Lymph node metastasis, and LDH below the median as features. (b) Model AUC for the ipi-naive subgroup. Each point is labeled with the additional feature, and the final model includes ploidy, heterogeneity, and purity as features. (c) Model AIC/BIC for the ipi-experienced subgroup. (d) Model AIC/BIC for the ipi-naïve subgroup. (e) Using 10-fold crossvalidation in the ipilimumab-naïve subset, mean crossvalidation AUC is 0.73, with a SD of 0.25. Each fold split the cohort into a training (n = 74 patients) and test set (n = 10 patients) ((f) In the ipi-experienced subset, mean crossvalidation AUC is 0.85, with a SD of 0.18. Each fold split the cohort into a training (n = 27 patients) and test set (n = 7 patients).

Extended Data Fig. 9 Subgroup Model ROC curves applied to other subgroups.

(a) Ipilimumab-treated subgroup model applied to ipilimumab-naive (n = 74 patients) subgroup; (b) Ipilimumab-naive model applied to ipilimumab-treated (n = 45 patients) subgroup; (c) Ipi-naïve model predicting PD as best response to PD-1 ICB applied to pre-ipilimumab-treated (n = 15 patients) tumors. (top) Comut plot overlaying best response and the ipi-naïve model score. Each column is one of 15 ipilimumab-treated tumors biopsied prior to ipilimumab therapy, ordered by the predictive model score. The 5 highest scoring tumors were all PD. (bottom) Receiver-operator curve for our ipilimumab-naive predictive model applied to pre-ipilimumab tumors (n = 15), with an AUC of 0.71.

Extended Data Fig. 10 Limited validation testing single feature models in validation cohort.

(a) Relationship between response and MHC-II HLA score in the ipilimumab-treated subgroup of a validation cohort. Each column is a patient, and patients are sorted by MHC-II HLA score. (b) Validation of Predictive Model for ipilimumab-treated tumors. In a validation cohort (n = 23), a model using the only common available feature (MHC-II HLA score) had an AUC of 0.65 (log likelihood ratio p = 0.09, empiric p = 0.21, Methods). (c) Relationship between response and MHC-II HLA score in the ipilimumab-treated subgroup of a validation cohort. Each column is a patient, and patients are sorted by MHC-II HLA score. (d) Validation of Predictive Model for ipilimumab-naive Tumors. In a validation cohort (n = 20), a model using the only common available feature (heterogeneity) had an AUC of 0.73 (log likelihood ratio p = 0.17, empiric p = 0.066, Methods).

Supplementary information

Reporting Summary

Supplementary Data 1

Gene-level mutation calls; copy number alterations; immune cell signatures; mutational signatures; and ssGSEA signatures using Hallmarks gene sets.

Supplementary Data 2

RNA-seq TPM matrix.

Supplementary Tables

Supplementary Tables 1–9.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Liu, D., Schilling, B., Liu, D. et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med 25, 1916–1927 (2019). https://doi.org/10.1038/s41591-019-0654-5

Download citation

Received: 25 March 2019
Accepted: 17 October 2019
Published: 02 December 2019
Issue Date: December 2019
DOI: https://doi.org/10.1038/s41591-019-0654-5

This article is cited by

Quantified pathway mutations associate epithelial-mesenchymal transition and immune escape with poor prognosis and immunotherapy resistance of head and neck squamous cell carcinoma
- Yuhong Huang
- Han Liu
- Chao Liu
BMC Medical Genomics (2024)
Mutation status of the KMT2 family associated with immune checkpoint inhibitors (ICIs) therapy and implicating diverse tumor microenvironments
- Dong-Xu Wang
- Jun-Yu Long
- Tao Li
Molecular Cancer (2024)
Conserved methylation signatures associate with the tumor immune microenvironment and immunotherapy response
- Qingqing Qin
- Ying Zhou
- Qiyuan Li
Genome Medicine (2024)
Integrated multi-omic analysis and experiment reveals the role of endoplasmic reticulum stress in lung adenocarcinoma
- Ying Liu
- Wei Lin
- Gaoren Wang
BMC Medical Genomics (2024)
An atlas of cell-type-specific interactome networks across 44 human tumor types
- Zekun Li
- Gerui Liu
- Yang Yang
Genome Medicine (2024)

Subjects

Abstract

Similar content being viewed by others

Main

Results

Genomic and clinical cohort characteristics and melanoma subtypes

Genomic and transcriptomic features associated with response

Correlations between molecular features

Previous exposure to anti-CTLA4 ICB

Integrative predictive modeling of primary resistance

Discussion

Methods

Patient cohort and clinical end points

Samples

Whole-exome and whole-transcriptome sequencing

Quality control and variant calling

Power calculation quality control

Variant calling

TMB and neoantigen load

Copy number variants

Purity and ploidy

Heterogeneity and aneuploidy

Mutational signatures

Transcriptomic analysis

GSEA and ssGSEA

Absolute immune infiltrate and immune subsets

Gene expression signatures

Analysis

Validation

Predictive model generation and cross-validation

Permutation testing

Reporting Summary

Data availability

Code availability

Change history

19 June 2020

References

Acknowledgements

Author information

Authors and Affiliations

Contributions

Corresponding authors

Ethics declarations

Competing interests

Additional information

Extended data

Supplementary information

Rights and permissions

About this article

Cite this article

Share this article

This article is cited by

Search

Quick links