Skip to main content
NCBI home page
NIHPA Author Manuscripts logoLink to NIHPA Author Manuscripts
. Author manuscript; available in PMC: 2019 Apr 29.
Published in final edited form as: Nat Med. 2018 Aug 20;24(10):1550–1558. doi: 10.1038/s41591-018-0136-1

Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response

Peng Jiang 1,2,#, Shengqing Gu 3,#, Deng Pan 4,5,#, Jingxin Fu 6, Avinash Sahu 1,2, Xihao Hu 1,2, Ziyi Li 6, Nicole Traugh 3, Xia Bu 3, Bo Li 1,2,9, Jun Liu 7, Gordon J Freeman 3, Myles A Brown 3,8, Kai W Wucherpfennig 4,5,11,*, X Shirley Liu 1,2,6,8,11,*
PMCID: PMC6487502  NIHMSID: NIHMS1024718  PMID: 30127393

Abstract

Cancer treatment by immune checkpoint blockade (ICB) can bring long-lasting clinical benefits, but only a fraction of patients respond to treatment. To predict ICB response, we developed TIDE, a computational method to model two primary mechanisms of tumor immune evasion: the induction of T cell dysfunction in tumors with high infiltration of cytotoxic T lymphocytes (CTL) and the prevention of T cell infiltration in tumors with low CTL level. We identified signatures of T cell dysfunction from large tumor cohorts by testing how the expression of each gene in tumors interacts with the CTL infiltration level to influence patient survival. We also modeled factors that exclude T cell infiltration into tumors using expression signatures from immunosuppressive cells. Using this framework and pre-treatment RNA-Seq or NanoString tumor expression profiles, TIDE predicted the outcome of melanoma patients treated with first-line anti-PD1 or anti-CTLA4 more accurately than other biomarkers such as PD-L1 level and mutation load. TIDE also revealed new candidate ICB resistance regulators, such as SERPINB9, demonstrating utility for immunotherapy research.

Cancer immunotherapies by immune checkpoint blockade (ICB) aim to help the immune system recognize and attack cancer cells1. The primary targets of ICB treatment are programmed death-ligand 1 (PD-L1), programmed death 1 (PD1) and cytotoxic T-lymphocyte-associated protein 4 (CTLA4). Compared to conventional therapies, ICB can induce durable responses in patients with metastatic cancers. However, a significant limitation of ICB is that only one-third of patients respond to ICB in most cancer types tested2. Combination ICB therapies have shown improved outcomes but also result in more severe side effects than single-agent therapy1. Multiple factors can affect ICB effectiveness2, including the degree of cytotoxic T cell infiltration3, mutation or neo-antigen load4, PD-L1 level5, antigen presentation defects6, interferon signaling7, mismatch repair deficiency8, tumor aneuploidy9 and intestinal microbiota10. However, none of these factors is sufficient to achieve accurate outcome prediction5, and identification of ICB response biomarkers and resistance regulators is a critical challenge in the field.

Gene expression biomarkers, such as Oncotype DX11, MammaPrint12 and Prosigna13, have demonstrated clinical utility in predicting treatment benefits in breast cancer. We hypothesize that transcriptome signatures could also serve as reliable ICB biomarkers. Ideally, a large number of tumor molecular profiles together with the patient clinical outcome could be used to train a reliable multi-gene biomarker. However, current ICB clinical trials have gene expression profiles on only a small number of pre-treatment samples, which are insufficient to train robust prognostic biomarkers3,14,15. Alternatively, there are many public tumor profiling data sets from human and mouse models without immunotherapy, but which are informative regarding tumor immune escape. For example, recent analyses of TCGA and PRECOG data uncovered significant effects of tumor-infiltrating levels of different immune cell types on patient overall survival16–18. Predicting tumor response to ICB requires an understanding of how tumors escape the immune system. Therefore, the public tumor molecular profiles, even without ICB treatment, may still be valuable resources to model immune evasion and derive surrogate biomarkers of ICB response.

Recent work has revealed two distinct mechanisms of tumor immune evasion19–20. Some tumors have a high level of infiltration by cytotoxic T cells, but these T cells tend to be in a dysfunctional state. In other tumors, immunosuppressive factors may exclude T cells from infiltrating tumors21. Therefore, we developed a computational framework, Tumor Immune Dysfunction and Exclusion (TIDE), to identify factors that underlie these two mechanisms of tumor immune escape. TIDE integrated and modeled data from 189 human cancer studies, comprising a total of 33,197 samples. We hypothesized and validated that an accurate gene signature to model the tumor immune escape could serve as a reliable surrogate biomarker to predict ICB response. The web application, source code and analysis results of TIDE are available at http://tide.dfci.harvard.edu.

Results

A statistical interaction test identifies gene signatures of T cell dysfunction

Previous reports showed that cytotoxic T cells could infiltrate a subset of tumors, although they could still fail to control tumor growth if in a dysfunctional state22. We reasoned that by combining transcriptome profiles of treatment-naive tumors with patient survival outcome, we could identify known regulators of T cell dysfunction. For example, in the TCGA melanoma study, we used the average expression level of CD8A, CD8B, GZMA, GZMB and PRF1 to estimate the cytotoxic T lymphocyte (CTL) level in a tumor16. Among metastatic melanoma tumors, a higher CTL level indicates a better patient survival, but only when TGFB1 has a low expression level (Fig. 1a and Supplementary Fig. 1a). This observation corroborates the known role of the cytokine TGFβ (encoded by TGFB1) in promoting tumor immune escape and immunotherapy resistance2,23,24.

Fig. 1 |. The interaction test identifies gene signatures of T cell dysfunction.

Fig. 1 |

Open in a new tab

a, The association between the CTL level and overall patient survival for melanoma tumors with different TGFB1 levels. For each metastatic melanoma tumor in TCGA, the CTL infiltration level was estimated as the average expression level of CD8A, CD8B, GZMA, GZMB and PRF1. The association between the CTL level and overall survival was computed through the two-sided Wald test in the Cox-PH regression. Each Kaplan–Meier plot presents tumors in two groups: ‘High CTL’ (red) have above-average CTL values among all samples, while ‘Low CTL’ (blue) have values below average. Samples were split according to the TGFB1 expression level to show the association between the CTL level and survival outcome. The top panel shows tumors with high TGFB1 expression (one standard deviation above the average), while the bottom panels show the remaining samples. b, The interaction test in a Cox-PH regression to identify genes associated with the T cell dysfunction. The variable CTL represents the level of CTLs in each tumor. The variable V represents the status of a candidate gene. The coefficient d reflects the effect of interaction between the CTL and V on death hazard outcome estimated from the survival data. The graphs represent the association slopes between CTL and death hazard. The black and gold arrows represent the association slopes before and after increasing the level of V. c, Genes with significant T cell dysfunction scores in multiple cancer types. Five data sets, representing five cancer types, had more than 1% of genes passing the FDR threshold 0.1. We display the genes whose T cell dysfunction scores, defined as the z score of d/standard error (s.e.), are significantly positive or negative (two-sided Wald test P values corresponding to an FDR less than 0.1) in at least two cancer types. The orange stars indicate genes of special interest. The number of samples in each cohort is available in Supplementary Table 2b. UCEC, uterine corpus endometrial carcinoma, TNBC, triple-negative breast cancer; AML, acute myeloid leukemia; SKCM, skin cutaneous melanoma; NB, neuroblastoma.

In statistics, two variables interact if the effect of one variable depends on the other variable25. In the previous examples, the effect of CTL on survival outcome depends on the TGFB1 level, which is a typical case of interaction between variables. The interaction of any two variables on survival outcome can be tested by a multiplication term in the Cox proportional hazard (Cox-PH) model26 (Fig. 1b). The coefficient d of the multiplication term indicates the level of the interaction effect, and the Wald test can evaluate its statistical significance26. For example, the TGFB1 expression level has a significantly antagonistic interaction with CTL level, indicating that a higher TGFB1 level in tumors will decrease the beneficial association between CTL and overall survival (Supplementary Table 1a). In contrast to TGFB1, another gene SOX10 expression level has a synergistic interaction with CTL level on patient survival outcome, indicating that a higher SOX10 level in tumors will increase the beneficial association between CTL and survival (Supplementary Table 1b), which is consistent with the known function of SOX10 to promote T cell-mediated tumor killing27,28.

We aim to systematically identify genes such as TGFB1 and SOX10 that influence the function of cytotoxic T cells on patient survival outcome in cancer genomics data cohorts. Using the Cox-PH model, TIDE tests how the interaction between a candidate gene V and the CTL affects death hazard (estimated from survival) (Fig. 1b). The resulting T cell dysfunctional signature is a genomewide vector, where the z score of each gene is the interaction coefficient d divided by its standard error (Supplementary Table 1). Genes with significant z scores are not restricted to genes expressed by T cells but could be expressed in cancer cells (for example, SOX1027,28) or different immune cells associated with T cell dysfunction. In the case of TGFB1, both cancer cells29 and CD4+ FOXP3+ Treg cells30 can express the cytokine TGFβ to inhibit T cell function.

To compute the T cell dysfunction scores in different cancer data sets, we collected hundreds of data sets from the TCGA31, PRECOG17 and METABRIC32 databases, and focused on 73 that had a minimum of 50 samples with both tumor expression profiles on the genome scale and patient survival data (Supplementary Table 2a). Among the data sets, TIDE predicted different numbers of genes to interact with CTL with statistical significance. For example, the P-value distribution for genes in TCGA melanoma was skewed to the left, indicating many significant genes (Supplementary Fig. 1b). However, the peak of significant P values was absent in TCGA glioblastoma. This difference is likely due to differences in T cell infiltration, data quality or sample size. In five data sets, over 1% of genes have significant interaction with CTL to affect survival at a false discovery rate (FDR) cutoff of 0.1: melanoma, neuroblastoma, triple-negative breast cancer, endometrial cancer and acute myeloid leukemia (Supplementary Fig. 1b and Supplementary Table 2b). For visualization, genes with significant dysfunction scores (FDR < 0.1) in at least two cancer types are shown in Fig. 1c (Supplementary Table 3). Although some of the genes are known to regulate T cell-mediated tumor immunity, such as PD-L1, others are likely to be co-expressed with immune-suppressive genes.

The TIDE dysfunction scores are consistent with signatures of tumor immune evasion

We evaluated the quality of TIDE T cell dysfunction scores using published studies of tumor immune evasion in pre-clinical models (Supplementary Table 4). One shRNA screen identified positive or negative hit genes whose knockdown in T cells enhanced or decreased T cell accumulation in mouse tumors, respectively33. Gene expression profiles to study T cell dysfunction are also publicly available, including the transcriptome of exhausted CD8 T cells34, activated regulatory T cells35 and tumors with acquired ICB resistance36. The positive or negative hits are defined as genes upregulated or downregulated in the process of T cell dysfunction or ICB resistance, respectively (Supplementary Tables 4 and 5). We examined whether the TIDE T cell dysfunctional signatures give significantly different scores between positive and negative hit genes in these published studies. We found that TIDE dysfunction signatures averaged from the five clinical cohorts assign positive hits significantly higher dysfunction scores compared to the negative hits (Fig. 2a). Using receiver operating characteristic (ROC) curves, we found that averaging the TIDE dysfunction signatures from the five cohorts gave the best performance (Fig. 2b,c and Supplementary Fig. 2a), suggesting the average profile as a more robust dysfunctional signature.

Fig. 2 |. T cell dysfunction signatures are consistent with published signatures of tumor immune evasion.

Fig. 2 |

Open in a new tab

a, The consistency between T cell dysfunction signatures predicted by the interaction test and published gene signatures of tumor immune evasion. To evaluate the reliability of the T cell dysfunction gene scores, we collected four published gene signatures related to T cell dysfunction and immunotherapy resistance (Supplementary Table 4). We plotted the T cell dysfunction scores averaged across five cancer types (average profile in Fig. 1c) for the positive (red) and negative (blue) hits of each gene signature. The numbers of positive and negative hits for each signature are available in Supplementary Table 4. Within each group, the scattered dots represent all gene values, and the thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. The difference between positive and negative groups was compared through the two-sided Wilcoxon rank-sum test, and P values for signatures of ‘T accum’, ‘T exhaust’, ‘T regulatory’ and ‘ICB resist’ are 6.94×10−3, 1.20×10−8, 1.81 ×10−6 and 1.95 ×10−7, respectively. The range of P values are labeled above each boxplot with asterisks (**P<1 ×10−2; ***P<1 ×10−3). T accum: shRNA screens for regulators of T cell accumulation in tumors; T exhaust: transcriptome of exhausted T cells; T regulatory: transcriptome of CD4 regulatory T cells. ICB resist: transcriptome of murine tumors that resist anti-CTLA4 checkpoint blockade. b, The ROC curves measuring the performance of the average T cell dysfunction scores (average profile in Fig. 1c) in predicting the positive and negative gene hits in each signature in a. c, The area under the ROC curve of the average profile of all five cancer types (black squares) and each of the individual cancer types SKCM, AML, NB, UCEC and TNBC with different dot colors. d, Pearson correlations between the T cell dysfunction scores and the expression profile of exhausted T cells. The correlations were computed across 12,498 genes shared between human and mouse signatures, for all pairwise combinations between five human cancer types and different time points in a mouse model of T cell exhaustion (‘T exh fixed’ in Supplementary Table 4).

Recent studies in mouse tumor models revealed two stages of T cell dysfunction37,38. While anti-PD1 treatment can revive the early-stage dysfunctional T cells, late-stage dysfunctional T cells are resistant to ICB reprogramming. The average profile of TIDE dysfunction signatures derived from the five cancer cohorts shows increasing correlation with the gene expression profiles of dysfunctional T cells in late stages38 (Fig. 2d). This result suggests that the TIDE dysfunction signatures reflect the profiles at the late stage of T cell dysfunction. We also applied gene set enrichment analysis to analyze the functional enrichment of TIDE T cell dysfunction signatures. Immune pathways related to inflammatory and interferon response are highly enriched, while mTORC1 signaling39, protein secretion40 and glycolysis41 that are known to promote CD8 T cell activation are consistently depleted (Supplementary Fig. 2b).

Immunosuppressive cell signatures predict immune escape by T cell exclusion

The previous section described T cell dysfunction signatures in tumors with high cytotoxic T cell infiltration. Next, we explored gene signatures of immune evasion through T cell exclusion in tumors with low T cell infiltration19,20. Several molecular mechanisms might explain the lack of T cell infiltration in the tumor, such as impaired priming of tumor-specific T cells or suppressive cells prohibiting T cell infiltration into the tumor19,20. To model the gene expression signature of T cell exclusion, we examined three cell types reported to restrict T cell infiltration in tumors, namely cancer-associated fibroblasts (CAFs), myeloid-derived suppressor cells (MDSCs) and the M2 subtype of tumor-associated macrophages (TAMs)20. We derived a genome-wide signature of T cell exclusion using expression profiles of these cell types from the Gene Expression Omnibus database42 (Supplementary Table 4). In TCGA melanoma data, tumors whose expression profiles have a higher correlation with the MDSC, TAM or CAF signatures show a significantly lower level of CTLs (Fig. 3a). Moreover, using the average expression profile of MDSCs, TAMs and CAFs to model T cell exclusion, we observed a strong negative correlation between the T cell exclusion scores and the CTL levels across tumors (Fig. 3a). Moreover, the CTL levels and T cell exclusion scores were negatively correlated in all solid tumor data sets (Fig. 3b).

Fig. 3 |. Immunosuppressive cell expression models gene signatures of T cell exclusion.

Fig. 3 |

Open in a new tab

a, Prediction of T cell exclusion scores for tumors using immunosuppressive cell signatures. For each metastatic tumor in the TCGA melanoma data set (blue dots, n = 317), we computed the Pearson correlation between its expression profile and the expression signature of MDSCs, M2 TAMs or CAFs (Supplementary Table 4) or the average of the three expression signatures. In each graph, these values are plotted along the x axis. The y axis shows the CTL level for each tumor (average expression level of CD8A, CD8B, GZMA, GZMB and PRF1). The Pearson correlation (R) between the plotted values is shown in the upper right corner of each plot. The two-sided t-test P values for correlations in MDSCs, TAMs, CAFs and mean are 4.61 ×10−37, 1.58×10−51, 8.84×10−13 and 2.58×10−52, respectively. b, A histogram of the correlations between the CTL levels and the T cell exclusion scores across tumors. The correlations analyzed in the histogram correspond to the R value in the top right corner of a (example of TCGA melanoma) across 43 solid tumor data sets. Gliomas are excluded because of low T cell infiltration levels in most gliomas16. c, Anti-correlation between T cell dysfunction scores and exclusion scores across TCGA melanoma tumors. For each metastatic melanoma tumor (colored dots, n = 317), we computed the Pearson correlation between the sample’s expression profile and the TIDE T cell dysfunction signature (y axis). The same computation was made between the tumor expression profile and the TIDE T cell exclusion signature (x axis). The Pearson correlation between the plotted values is shown in the upper right (two-sided t-test P value = 4.02 ×10−34). The dot color indicates the CTL level in each tumor. d, Anti-correlation between T cell dysfunction scores and exclusion scores across TCGA cancer types. For each TCGA cancer type with normal control samples (n = 17), we calculated the average expression difference between tumor versus normal samples. We then computed the Pearson correlation between that value and the TIDE T cell dysfunction signature (y axis). We also made the same calculation for the TIDE T cell exclusion signature (x axis). The Pearson correlation between the plotted values is shown in the upper right (one-sided t-test P value = 0.042). The CTL level difference between tumor and normal samples is shown by the dot color.

We further examined the associations between the gene signatures of T cell exclusion and T cell dysfunction. For each tumor, the enrichment of a signature is computed as the Pearson correlation between the tumor gene expression profile and the genomewide scores of T cell exclusion or dysfunction signatures. In the five cancer types where we can identify significant T cell dysfunction scores, the level of T cell exclusion in a tumor inversely correlates with the level of T cell dysfunction (Fig. 3c and Supplementary Fig. 3a). We also calculated the signature enrichment based on the differential expression between tumor and normal samples across TCGA cancer types and observed similar negative correlations between T cell exclusion and T cell dysfunction (Fig. 3d and Supplementary Table 6). Kidney renal cell carcinoma (KIRC) has the highest CTL level and the highest enrichment of the T cell dysfunction signature (Fig. 3d and Supplementary Fig. 3b), while lung squamous carcinoma (LUSC) has the highest T cell exclusion signature (Fig. 3d and Supplementary Fig. 3c). Our results are consistent with previous reports of a high CTL level in KIRC and a low CTL level in LUSC16. These results suggest that the KIRC and LUSC tumors utilize distinct immune evasion strategies, with KIRC operating more through T cell dysfunction and LUSC through T cell exclusion. Previous studies reported paradoxical observations that in KIRC the degree of CD8 cytotoxic T cell infiltration is anti-correlated with survival benefits43. Our analysis revealed that KIRC tumors with higher CTL levels tend to have a stronger T cell dysfunction signature, which could impair the ability of cytotoxic T cells to kill cancer cells (Supplementary Fig. 3b).

TIDE signature predicts ICB response

In previous sections, we developed genome-wide expression signatures to measure the level of T cell dysfunction and T cell exclusion in tumors. We next examined whether integrating these two signatures could predict ICB clinical response. Among the five cancer types for which we could compute reliable TIDE signatures (Fig. 1c), only melanoma has publicly available data on tumor expression and clinical outcome of patients treated with anti-PD114 or anti-CTLA43, so it was the focus of our evaluation. We also evaluated TIDE on an anti-PD1 data set that profiled tumor expression profiles across four cancer types using the NanoString assay on a few hundred genes15.

We classified the tumors as CTL-high if the expression levels of all CTL markers (CD8A, CD8B, GZMA, GZMB and PRF1) were higher than their average values in each data set and the remaining tumors as CTL low. In the CTL-high tumors, TIDE correlates the tumor expression data with the T cell dysfunction signature and predicts tumors with high correlation to T cell dysfunction as non-responders (Supplementary Fig. 4a). In CTL-low tumors, it has been reported that ICB can enhance the cytotoxic T cell infiltration44,45, so patients with low tumor CTL might still derive clinical benefits from immunotherapies. Therefore, TIDE correlates the expression data for each tumor with the T cell exclusion signature in CTL-low tumors and predicts those with suppressive cells inhibiting T cell infiltration as non-responders (Supplementary Fig. 4a). Notably, the correlation between tumor expression profiles and TIDE signatures is a single value computed across all human genes (Supplementary Fig. 4b), and therefore not subject to multiple-hypotheses testing and less sensitive to the noise from individual expression or the TIDE signature value. For the pre-treatment transcriptome of each tumor, the Pearson correlation with either T cell dysfunction (in CTL-high tumors) or exclusion (in CTL-low tumors) signatures was defined as the TIDE prediction score (Fig. 4a–c). Correlations with T cell dysfunction or exclusion signatures may have different distributions (Supplementary Fig. 4c). Thus, when merging the prediction scores from two tumor CTL categories, we normalized the correlations by their standard deviations in the TCGA data. Finally, all tumors were ranked by their TIDE scores to predict their ICB response (Fig. 4a–c and Supplementary Fig. 4d).

Fig. 4 |. TIDE signatures predict iCB immunotherapy response.

Fig. 4 |

Open in a new tab

a, A waterfall plot of TIDE prediction scores across 25 melanoma tumors treated with anti-PDI14. The TIDE framework divided tumors into CTL-high or -low categories based on the expression level of CTL marker genes (Supplementary Fig. 4a). Red indicates a tumor that responded to therapy. Blue indicates a non-responder. In each category, we sorted tumors in descending order according to their TIDE prediction scores. b, A waterfall plot of TIDE prediction scores across 42 melanoma tumors treated with anti-CTLA43 in the same way as in a. Besides responders or non-responders, several patients are classified as long-survival in the original study due to the long overall survival time3. c, A waterfall plot of TIDE prediction scores across 33 tumors treated with anti-PD115 in the same way as in a. The gene expression profiles are measured by the NanoString platform. The 33 tumors comprise 9 melanoma, 12 lung adenocarcinoma, 9 lung squamous carcinoma and 3 head and neck tumors. d, ROC curves for the performance of the TIDE prediction score, PD-L1 expression, interferon gamma (IFNG) response and total mutation load in predicting anti-PD1 response among 25 melanoma tumors in a. e, ROC curves for the performance of several signatures in predicting anti-CTLA4 response among 42 melanoma tumors in b. f, ROC curves for the performance of several signatures in predicting anti-PD1 response among 33 tumors in c. g, The area under the ROC curve (AUC) for several signatures in predicting anti-PD1 response among 25 melanoma tumors in a. The signatures are defined in Supplementary Table 7, with TIDE in dark red and other signatures in blue. Besides the genome-wide TIDE signature, ‘TIDE.selected’ is a variation focused on 770 genes with both high expression variation across tumors and significant values in the either T cell dysfunction or exclusion signatures. The performance of a random predictor (AUC = 0.5) is represented by the dashed line. h, AUC for signatures in predicting anti-CTLA4 response among 42 melanoma tumors in b in the same way as in g. i, AUC for signatures in predicting anti-PD1 response among 33 tumors in c in the same way as in g. TIDE AUC metrics are also shown separately for nine melanoma, twelve lung adenocarcinoma (Adeno) and nine lung squamous carcinoma (Squamous) tumors. j, Kaplan–Meier plots of overall survival (OS) for 25 melanoma patients (in a) treated with anti-PD1 with the top (>1) and bottom (<1) TIDE prediction scores. The P value was calculated by testing the association between TIDE prediction scores and overall survival with the two-sided Wald test in a Cox-PH regression. k, Kaplan-Meier plots of overall survival for 42 melanoma patients (in b) treated with anti-CTLA4 in the same way as in j. l, Kaplan-Meier plots of progression-free survival (PFS) for 33 patients (in c) treated with anti-PD1 in the same way as in j.

To evaluate the prediction performance for ICB response, we used ROC to measure the true-positive rates against the falsepositive rates at various thresholds of TIDE prediction scores (Fig. 4d–f). Compared to widely used ICB response biomarkers, tumor mutation load, PD-L1 level and interferon gamma response5,7, the TIDE signature achieved consistently better performance for both anti-PD1 and anti-CTLA4 therapies using both RNA-Seq and NanoString data (Fig. 4d–f and Supplementary Fig. 5a). We also compared TIDE with other ICB response signatures reported in the literature (Supplementary Table 7). Among all candidate biomarkers, we found the TIDE signature to be the best predictor for both anti-PD1 and anti-CTLA4 therapies (Fig. 4g–i and Supplementary Table 8a). The prediction performance of TIDE is also higher than the signatures of T cell dysfunction and ICB resistance discussed in Fig. 2a (Supplementary Fig. 5b). Meanwhile, the TIDE prediction performance is robust against modest variations of CTL cutoff in the definition of CTL-high or -low tumors (Supplementary Fig. 5c). Moreover, a higher tumor TIDE prediction score is associated not only with worse ICB response, but also with worse patient survival under anti-PD1 and anti-CTLA4 therapies (Fig. 4j–l). One explanation for the better performance of TIDE is that TIDE utilized both T cell dysfunction and exclusion signatures to model immune escape in tumors with different CTL levels, while other biomarkers consider only one aspect (Supplementary Fig. 5d–f and Supplementary Table 8b,c). Paradoxically, a previous computational method ImmunoPhenoScore claimed to have 100% accuracy in predicting ICB response in melanoma46, but we observed considerably lower accuracy of ImmunoPhenoScore using the source codes provided by the authors (Fig. 4g–i).

Besides the anti-PD1 RNA-Seq cohort14 analyzed in Fig. 4, a recent study generated RNA-Seq profiles on another melanoma cohort treated with anti-PD145. We focused on 24 patients with genomics profiles (expression and mutation) of pre-treatment tumors and anti-PD1 as the first-line immunotherapy (without previous anti-CTLA4 therapy). While the TIDE prediction score has a similar accuracy to the mutation load (Supplementary Fig. 6a,b), it is significantly predictive of the patient overall survival (‘Ipi naive’ in Supplementary Fig. 6c), demonstrating its prognostic value. We noted that TIDE achieved a lower prediction performance in the Riaz study compared to its performance in the Hugo study (Fig. 4d versus ‘Ipi naive’ in Supplementary Fig. 6). A possible explanation is that the Riaz study45 used the RECIST v1.1 criteria for disease progression, while the Hugo study14 used the immune-related RECIST47, which is more appropriate for predicting immunotherapy response. Further, TIDE is trained using data from ICB-naive tumors and thus not relevant in modeling the tumors that progressed after a first-line ICB2 (‘Ipi progressed’ in Supplementary Fig. 6).

The TIDE dysfunction score predicts regulators of ICB resistance

We hypothesized that some genes with high scores in TIDE signatures might serve not only as biomarkers but also as ICB resistance regulators. We focused on the T cell dysfunction signature for genes regulating T cell dysfunction in tumors. As the T cell dysfunction scores were computed using the data from treatment-naive tumors, we utilized orthogonal data from a mouse model of acquired anti-CTLA4 resistance to identify genes that are associated with ICB resistance36. We ranked genes with significant T cell dysfunction scores in Fig. 1c by the expression fold-change in the ICB-resistant tumors36 and identified Serpinb9 as the most upregulated gene (Fig. 5a,b). In ICB clinical cohorts, the SERPINB9 expression level is consistently lower in responders than non-responders (Supplementary Fig. 7a). Moreover, SERPINB9 expression alone is significantly associated with worse overall survival in two clinical studies of anti-CTLA4 therapy3,48 (Fig. 5c, Supplementary Fig. 7b and Supplementary Table 9).

Fig. 5 |. Validation of SERPINB9 as a regulator of tumor immune escape.

Fig. 5 |

Open in a new tab

a, The log-fold change (log[FC]) of expression between anti-CTLA4-resistant and parental B16 murine tumors for genes with significant T cell dysfunction scores in Fig. 1c. All genes are ranked increasingly with the top one labeled by name. b, The expression value of Serpinb9 between anti-CTLA4-resistant and parental B16 tumors. Within each group, the scattered dots represent Serpinb9 expression values (n = 6 samples in the resistant group, n = 4 samples in the parental group). The thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. The P value, testing the group difference, was calculated with the two-sided Wilcoxon rank-sum test. c, Kaplan-Meier plots of patients with top half and bottom half SERPINB9 expression levels, using the data from an anti-CTLA4 study with 42 patients profiled3. Both progression-free survival and overall survival are shown. The association between SERPINB9 expression and patient survival was tested by the two-sided Wald test in a Cox-PH regression (Supplementary Table 9). d, Western blot of SERPINB9 following genetic knockout and overexpression. For knockout (KO), there are two independent CRISPR guides targeting Serpinb9 and a control non-targeting sequence. Cells were either untreated (left 3 lanes) or treated with 100 ng ml−1 IFNγ to induce Serpinb9 expression (right 3 lanes). VCL is the loading control. For overexpression, the open reading frame of Serpinb9 was cloned into the pEF1a plasmid and overexpressed in B16F10 cells. The protein level was compared between pEF1a backbone- and pEF1a-Serpinb9-transduced cells. The bands of related protein targets are cut and shown. All experiments have been repeated independently two times with similar results. e, The effect of Serpinb9 knockout on T cell-mediated tumor killing. B16F10 cancer cells were co-cultured for three days with cytotoxic T cells at three B16F10 to T cell ratios (3:1, 2:1 or 1:1). Each CRISPR gRNA-transduced GFP positive cell line (Control, KO 1, KO 2) was mixed with the parental GFP-negative cell line at a 1:1 ratio. After co-culture, the ratio of edited GFP+cells to parental cells (GFP−) was determined by flow cytometry. The bar plots present the median value among three cell-culture replicates with standard deviations as the error bars. The results of the two-sided Student t-test, comparing the difference between knockout and control conditions, are available in Supplementary Table 10. f, The effect of Serpinb9 overexpression on T cell-mediated tumor killing. The effect of Serpinb9 overexpression was examined. The bar plots present the median value among two cell-culture replicates with standard deviations as the error bars. The results of the two-sided Student t-test, comparing the difference between overexpression and control conditions, are available in Supplementary Table 10.

SERPINB9 is a member of the serine protease inhibitor (serpin) family. The encoded protein can inactivate granzyme B to protect lymphocytes (for example, T cells and natural killer cells) from granzyme that may leak from the granules49. It is normally expressed in cytotoxic lymphocytes, antigen-presenting cells and immune-privileged sites50–52. Meanwhile, a study using in vitro cell culture models reported that a high SERPINB9 level in cancer cells resulted in resistance to T cell-mediated killing53. To infer which cell type in tumors is the potential source of a high SERPINB9 level, we examined the Protein Atlas database of immunohistochemistry results for 15,000 genes in 20 cancer types54. The SERPINB9 protein is expressed at a higher level in cancer cells of melanoma and several other cancer types as compared to normal tissues (Supplementary Fig. 7c,d). Thus, SERPINB9 may promote the resistance to T cell-mediated killing during ICB therapy through its high expression in cancer cells.

To validate SERPINB9 function in cancer cells, we knocked out Serpinb9 using CRISPR–Cas9 in the murine B16F10 melanoma cell line, which is the parental line of the anti-CTLA4-resistant tumor model previously discussed36. After knocking out Serpinb9 using two different CRISPR guide RNAs (gRNAs), the protein level became undetectable (Fig. 5d and Supplementary Fig. 8). When co-cultured with Pmel-1 cytotoxic T cells, the Serpinb9-knockout B16F10 cells were more sensitive to T cell-mediated killing compared to control cells (Fig. 5e, Supplementary Fig. 9a and Supplementary Table 10). In contrast, B16F10 cells with Serpinb9 overexpression were significantly more resistant to T cell-mediated killing compared to control cells (Fig. 5f, Supplementary Fig. 9b and Supplementary Table 10).

Notably, in B16F10 cells, the SERPINB9 protein level is significantly increased on treatment with IFNγ, a cytokine produced by cytotoxic T cells following antigen-specific activation55 (Fig. 5d). This SERPINB9 induction following IFNγ treatment might be explained by the binding of IRF1, a transcription factor activated by IFNγ signaling56, near the Serpinb9 gene that is observed in public ChIP-Seq data sets (Supplementary Fig. 10a,b). In human melanoma tumors, the expression level of SERPINB9 is highly correlated with IRF1 on both bulk tumor and single-cell levels (Supplementary Fig. 10c,d). These results support that SERPINB9 in cancer cells regulates resistance to T cell-mediated killing, which is essential for ICB response. Interestingly, the SERPINB9 expression level is also consistently upregulated following pathogen infection in curated studies from the NCBI Gene Expression Omnibus42 and Expression Atlas57 databases (Supplementary Fig. 11 and Supplementary Table 11). This result indicates that SERPINB9 is potentially a general regulator of immune evasion utilized by both tumors and pathogens.

Discussion

We developed a computational method called TIDE, which integrates the expression signatures of T cell dysfunction and T cell exclusion to model tumor immune evasion. The TIDE signatures, trained from treatment-naive tumor data, can predict ICB clinical response based on pre-treatment tumor profiles. Furthermore, TIDE predicted regulators of ICB resistance whose inhibition might improve ICB response. We experimentally validated the role of SERPINB9, an inhibitor of the cytotoxic lymphocyte protease GZMB, in tumor immune evasion, which is an essential process of ICB resistance. Although no small-molecule inhibitor of SERPINB9 is available, the Pfizer OASIS database indicates this protein as potentially druggable58.

Of the 73 data sets analyzed in this study, only five gave statistically significant T cell dysfunction signatures from the interaction test (Supplementary Table 2). This result is partly because we considered only 48 out of the 73 data sets where a higher level of tumor-infiltrating cytotoxic T cells is correlated with better survival outcome. In some cancer types, such as renal cell carcinoma, which has a substantial level of CD8 T cell infiltration, higher CTL may not be associated with survival benefits43. Also, depending on the sample size or characteristics of specific data sets, there might not be any statistically significant genes interacting with CTL to influence survival. Since averaging signatures from the five data sets yielded a signature more robust than any individual signature, integrating additional cancer data sets in the future has the potential to further improve the robustness of the T cell dysfunction scores (Fig. 2c). With additional data, cancer-type-specific regulators may be identified on the basis of the biological variations of T cell dysfunction scores across different cancer types.

When using the TIDE model to predict ICB response, we determined a cutoff to classify tumors as CTL-high or CTL-low. We used the average expression of CTL markers (CD8A, CD8B, GZMA, GZMB and PRF1) across all tumors to determine the CTL threshold. However, if matched normal tissues are available, the CTL threshold could also be determined by comparing the CTL marker expression in tumors with marker expression in normal tissues. The TIDE signature consists of genome-wide scores of T cell dysfunction and exclusion. While a genome-wide transcriptome biomarker might be robust for ICB response prediction, RNA-Seq has not been widely adopted in the clinic. A smaller gene panel for qPCR or NanoString assays could be more clinically pragmatic. We demonstrated TIDE performance on an anti-PD1 response data set where baseline tumor expression was measured on the NanoString PanCancer panel (Fig. 4c).

One limitation of our study is that we focused primarily on gene expression biomarkers. However, other biomarker types can also predict T cell infiltration and ICB response. For example, beta-catenin protein level has a negative correlation with CTL in many cancer types (Supplementary Fig. 3d,e and Supplementary Table 12). Moreover, tumors initially responding to ICB may later acquire mutations, such as in B2M, IFNGR1/2 and JAK1/2 genes, to become ICB resistant2. However, our study focuses only on predicting intrinsic ICB resistance. Therefore, more data types and methods are necessary to model the immunotherapy efficacy comprehensively. It is possible that TIDE could be applied jointly with other types of biomarker to achieve a higher prediction performance.

To enable testing of TIDE by clinicians and the public, we created a web application for response prediction using transcriptome profiles at http://tide.dfci.harvard.edu. TIDE has the potential to help oncologists select patients who are more likely to benefit from ICB. It would be of significant interest to test the clinical utility of TIDE in ICB decision-making in a prospective clinical trial. New immune-oncology data are emerging at an increasingly rapid pace. We envision that computational modeling and data integration will be increasingly utilized to refine ICB response biomarkers and identify new immunotherapy targets.

Methods

Methods, including statements of data availability and any associated accession codes and references, are available at https://doi.org/10.1038/s41591-018-0136-1.

Methods

Data collection of clinical genomics studies

We collected cancer data sets with both patient survival durations and tumor gene expression profiles from the TCGA31, PRECOG17 and METABRIC32 databases. If the clinical information is available, we separated the breast cancer data sets into the PAM50 (Prediction Analysis of Microarray 50, Prosigna) subtypes13 of luminal A, luminal B, Her2 positive, basal and triple negative (a variation of basal subtype). This separation is because each PAM50 subtype has a distinct genomics profile59 and degree of cytotoxic T cell infiltration60. Among all TCGA cancers, melanoma has two major tumor types annotated (that is, primary and metastatic); thus, we split melanoma profiles into primary and metastatic subtypes. The PRECOG database provided only survival duration information without other clinical factors; thus, we cannot perform subtype analysis. METABRIC is a breast cancer cohort, and we split all tumors according to the PAM50 subtypes (luminal A, luminal B, HER2, basal and triple negative).

To ensure the robustness of our analysis, we excluded the data sets from microarray platforms with fewer than 15,000 genes or without probes for cytotoxicity T cell markers (CD8A, CD8B, GZMA, GZMB and PRF1). Also, we included only data sets with more than 50 patients and 10% death rate because a low event number may undermine the reliability of Cox-PH survival regression26. Finally, 73 data sets from 3 databases passed our selection criteria (Supplementary Table 2a). The expression values of all genes are normalized by subtracting the mean values across all samples in a data set.

Statistical analysis

The interaction test in multivariate Cox-PH regression was applied to identify gene association with T cell dysfunction phenotype. In statistics, two variables interact if the effect of one variable depends on the status of the other, and a multiplication term in a multivariate linear model can test the interaction effect between two variables25. We applied the Cox-PH survival regression to test how the level of CTL interacts with other genes in the tumor to affect survival outcome. We solve a linear model Hazard = a×CTL + b×V + d×CTL×V + c using the Cox-PH regression26. The CTL level is estimated through the bulk-tumor expression average of CD8A, CD8B, GZMA, GZMB and PRF1. In the Cox-PH model, the death hazard was estimated through the patient survival information. The variable V represents the expression level of a candidate gene in the test. Since we have selected data sets where CTL correlates with favorable survival outcome, the coefficient a is always negative. The association slope between CTL and Hazard is a + d×V (Fig. 1b). If the coefficient d is positive, a higher V level will flatten the slope between CTL and Hazard, indicating a reduced association between the cytotoxic T cell level and better survival outcome. If d is negative, a higher V level will sharpen the slope between CTL and Hazard, indicating an increased association between the cytotoxic T cell level and better survival outcome. The T cell dysfunction score for each gene is defined as the Wald test z score, which is the coefficient d divided by its standard error26 (Fig. 1c and Supplementary Table 1). Of note, the thresholds shown in Fig. 1a and Supplementary Fig. 1a are used only to illustrate the principle of statistical interaction used by the model. When computing the T cell dysfunction scores through regression, we used the continuous variables without any thresholds. Also, we included clinical covariates, such as age, gender and stage (if available), in the regression to control for potential confounding factors.

To identify significant genes in the interaction test, we applied the Benjamini-Hochberg method to convert the two-sided Wald test P values to FDRs61, and selected clinical data sets with more than 1% genes having an FDR smaller than 0.1. This procedure is equal to selecting data sets where the distribution of P values has a significant peak near zero62. For example, the P-value histogram computed using TCGA melanoma data has a spike near zero, indicating that a set of genes significantly interact with CTL to affect survival outcome (Supplementary Fig. 1b). In contrast, the result computed from glioblastoma data does not contain any genes with significant interactions (Supplementary Fig. 1b).

Performance comparison on predicting ICB response

We collected the RNA-Seq data in melanoma for anti-CTLA43 and anti-PD114 therapies with gene expression profiles for 25 and 42 pre-treatment tumors with complete clinical information, respectively. We collected the NanoString data for anti-PD1 therapies with gene expression profiles of 33 baseline tumors in four cancer types15. For each data set, we standardized the transcriptome data across patients by quantile-normalization, and further normalized the expression values of each gene by subtracting the average among all samples. Therefore, a zero value indicates the average expression.

To predict each tumor’s potential to escape T cell-mediated killing, we first classified each tumor into CTL-high or CTL-low categories through the CTL marker expression levels (CD8A, CD8B, GZMA, GZMB and PRF1). Tumors with all positive values (higher than average) are classified as CTL-high, while the rest as CTL-low (Supplementary Fig. 4a). For the CTL-high tumors, we computed the Pearson correlation between tumor gene expression profiles and the T cell dysfunction signature (Supplementary Fig. 4b). For the CTL-low tumors, we computed the Pearson correlation between tumor gene expression profiles and the T cell exclusion signature (Supplementary Fig. 4b). The correlation with T cell dysfunction or exclusion signatures may have different distributions (Supplementary Fig. 4c). Therefore, to make the scale of Pearson correlations comparable between CTL-high and -low tumors, we normalized the correlation values within each sub-category through the standard deviation of correlations computed using the TCGA melanoma data. The scaled correlations were defined as TIDE prediction scores, representing the potential of tumor immune escape (Fig. 4a–c).

We also computed the response prediction from other biomarkers published in the literature. The predicted values of gene expression biomarkers (for example, IFNG, CD8, PDL1 and CRMA) were the average values among all members defined by the original publications (Supplementary Table 7). The predicted values of ImmunoPhenoScore were computed using the source codes provided by the authors46. The predicted value of the tumor SCNA biomarker was downloaded from the original publication for the anti-CTLA4 data set9 and provided by W. Hugo for the anti-PD1 data set14.

The outcome predicted by all biomarkers is a range of values, instead of a binary outcome. For example, total mutation load, CD8 expression level and TIDE prediction score all give one value for each patient tumor instead of a response classification. Therefore, we utilized the ROC curve, which plots the true-positive rates versus false-positive rates at various thresholds of biomarker values (Fig. 4d–f). The area under the ROC curve was used as the quality metric of prediction (Fig. 4g–i).

Gene selection for a focused TIDE signature

We select the most informative genes with both high variance across tumors and significant values in the TIDE signature. We selected 770 genes because that number is compatible with a NanoString platform that could be designed for a clinical assay. In the first step, we computed the standard deviation of expression values across samples for all genes in each TCGA cancer data set and selected 4,150 genes whose standard deviation is higher than the average of all genes in more 80% TCGA data sets. Next, we ranked the 4,150 genes by their maximum absolute values in the TIDE signatures of T cell dysfunction and exclusion. From this ranked list, we selected the top 770 genes, which is the maximum number that can fit on a NanoString assay. The column ‘TIDE.selected’ in Fig. 4g–i shows the TIDE performance on selected genes.

T cell killing assay based on co-culture between B16 and T cells

B16F10 cells were maintained in complete DMEM media (10% FBS and 50 U ml−1 of penicillin/streptomycin). B16F10-Cas963 cells were maintained in complete DMEM media with 2.5–5 μg ml−1 of blasticidin. CD8 T cells isolated from mice were cultured in complete RPMI 1640 media (10% FBS, 20 mM HEPES, 1 mM sodium pyruvate, 0.05 mM 2-mercaptoethanol, 2 mM l-glutamine and 50 U ml−1 streptomycin and penicillin). All cell lines are tested for mycoplasma contamination. Pmel-1 TCR transgenic mice were purchased from Jackson Laboratory (stock no. 005023).

CD8 T cells were isolated from spleen and lymph nodes from Pmel-1 TCR transgenic mice using the EasySep mouse CD8+ T cell isolation kit (STEMCELL no. 19753) according to the manufacturer’s protocol. Freshly isolated CD8 T cells were stimulated with anti-CD3/CD28 beads (ThermoFisher no. 11452D) at a bead to cell ratio of 1:2 to induce differentiation into an effector state. On day 3, recombinant mouse IL-2 (Biolegend, no. 575406) was added to the culture at 20 ng ml−1. T cells were used for co-culture with B16F10 cells after at least six days of in vitro activation. Our animal experiments have complied with all relevant ethical regulations. The study protocol in this study was approved by the Institutional Care and Use Committee at Dana Farber Cancer Institute.

To knockout Serpinb9, CRISPR gRNA sequences targeting Serpinb9 or non-targeting control were cloned into a PLKO3G-GFP vector and confirmed by sequencing. To overexpress Serpinb9, its cDNA was amplified, cloned into a pEF1a-puro vector and confirmed by sequencing. Knockout or overexpression constructs were co-transfected with pCMV-dR8.91 and pCMV-VSV-G (Addgene no. 8454) into HEK293T cells to generate lentivirus. Transfection was performed using TransIT-293 (Mirus, MIR2700) following the manufacturer’s protocol. Lentivirus was collected 48 h later and stored at −80 °C. To generate Serpinb9-knockout cells, B16F10-Cas9 cells were infected with a lentivirus driving expression of a single gRNA overnight to inactivate Serpinb9 genes individually. Infected cells were sorted on the basis of GFP expression by BD FACS Aria II. To generate Serpinb9-overexpressing cells, B16F10-Cas9 cells were infected with either pEF1a-puro backbone or pEF1a-puro-Serpinb9. Infected cells were selected by culturing with 2 μg ml−1 puromycin. Control (non-targeting gRNA or pEF1a-puro backbone), Serpinb9-deficient or Serpinb9-overexpressing B16F10-Cas9 cells were lysed and subjected to western blot analysis with the following antibodies: anti-SERPINB9 (Santa Cruz Biotechnology no. sc-57531) and anti-VCL (Sigma Aldrich no. V9264).

In a competition assay with Serpinb9-knockout cells, Serpinb9-deficient or non-targeting guide control B16F10-Cas9 cells (GFP positive) were mixed with control B16F10-Cas9 cells (GFP negative) at a 1:1 ratio. In a competition assay with Serpinb9-overexpressing cells, pEF1a control or pEF1a-Serpinb9 B16F10-Cas9 cells (GFP negative) were mixed with control GFP-infected B16F10-Cas9 cells (GFP positive) at a 1:1 ratio. Mixed cells were stimulated with 10 ng ml−1 or 100 ng ml−1 of interferon gamma for 24 h to enhance MHC class I expression. These tumor cells were then co-cultured with in vitro-activated Pmel-1 T cells at different effector-to-target ratios in a 6-well plate. Tumor cells were plated at equal density in all wells, and T cells were added at 0, 1/3, 1/2 or 100% of the number of tumor cells. There are two or three cell-culture replicates for each condition. After a three-day co-culture with T cells, the fold-change of Serpinb9-edited B16F10 cells was determined by FACS, comparing the percentage of Serpinb9-edited B16F10 cells to control B16F10 cells (Supplementary Fig. 9). T cells present in these cultures were gated out on the basis of antibodies specific for CD45 (APC–Cy7) (Biolegend, 103115).

Reporting Summary

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

Software availability

The TIDE web application of response prediction is freely accessible with any modern web browser through http://tide.dfci.harvard.edu/. We will keep the website and tool operating and freely accessible for the foreseeable future. The source code for computing the T cell dysfunction score through the interaction test is available under GNU Public License v3 through GitHub at https://github.com/foreverdream2/dysfunction_interaction_test.

Data availability

Users can query our analysis results with gene names: http://tide.dfci.harvard.edu/query/. All of our processed input data, analysis output data and an example script to repeat our major results are available at http://tide.dfci.harvard.edu/download/.

Supplementary Material

1

Acknowledgements

The research was supported by the Cancer Immunologic Data Commons (1U24CA224316-01) grant of the National Cancer Institute (NCI), the Pathway to Independence Award (1K99CA218900-01) grant of NCI (to P.J.), the Specialized Center (1P50CA206963-01) grant of NCI (to G.J.F.), and the Breast Cancer Research Foundation (to X.S.L.). D.P. is a Cancer Research Institute/Robertson Foundation Fellow.

Competing interests

X.S.L. is a cofounder and board member of GV20 Oncotherapy, a scientific advisor of 3DMedCare and a paid consultant for Genentech. K.W.W. is a member of the scientific advisory board for TCR2 and Nextech; he serves as a consultant for Novartis. The laboratory of K.W.W. received sponsored research funding from Astellas Pharma Inc.

Footnotes

Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41591-018-0136-1.

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

References

  • 1.Mahoney KM, Rennert PD & Freeman GJ Combination cancer immunotherapy and new immunomodulatory targets. Nat. Rev. Drug Discov 14, 561–584 (2015). [DOI] [PubMed] [Google Scholar]
  • 2.Sharma P, Hu-Lieskovan S, Wargo JA & Ribas A Primary, adaptive, and acquired resistance to cancer immunotherapy. Cell 168, 707–723 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 3.Van Allen EM et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science 350, 207–211 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 4.Snyder A et al. Genetic basis for clinical response to CTLA-4 blockade in melanoma. New Engl. J. Med 371, 2189–2199 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 5.Nishino M, Ramaiya NH, Hatabu H & Hodi FS Monitoring immune-checkpoint blockade: response evaluation and biomarker development. Nat. Rev. Clin. Oncol 14, 655–668 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 6.Zaretsky JM et al. Mutations associated with acquired resistance to PD-1 blockade in melanoma. N. Engl. J. Med 375, 819–829 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 7.Ayers M et al. IFN-gamma-related mRNA profile predicts clinical response to PD-1 blockade. J. Clin. Investig 127, 2930–2940 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 8.Le DT et al. PD-1 blockade in tumors with mismatch-repair deficiency. N. Engl. J. Med 372, 2509–2520 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 9.Davoli T, Uno H, Wooten EC & Elledge SJ Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science 355, eaaf8399 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 10.Kroemer G & Zitvogel L Cancer immunotherapy in 2017: The breakthrough of the microbiota. Nat. Rev. Immunol 18, 87–88 (2018). [DOI] [PubMed] [Google Scholar]
  • 11.Paik S et al. A multigene assay to predict recurrence of tamoxifen-treated, node-negative breast cancer. N. Engl. J. Med 351, 2817–2826 (2004). [DOI] [PubMed] [Google Scholar]
  • 12.van ‘t Veer LJ et al. Gene expression profiling predicts clinical outcome of breast cancer. Nature 415, 530–536 (2002). [DOI] [PubMed] [Google Scholar]
  • 13.Parker JS et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J. Clin. Oncol 27, 1160–1167 (2009). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 14.Hugo W et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell 165, 35–44 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 15.Prat A et al. Immune-related gene expression profiling after PD-1 blockade in non-small cell lung carcinoma, head and neck squamous cell carcinoma, and melanoma. Cancer Res 77, 3540–3550 (2017). [DOI] [PubMed] [Google Scholar]
  • 16.Rooney MS, Shukla SA, Wu CJ, Getz G & Hacohen N Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 17.Gentles AJ et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat. Med 21, 938–945 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 18.Li B et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol 17, 174 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 19.Gajewski TF, Schreiber H & Fu YX Innate and adaptive immune cells in the tumor microenvironment. Nat. Immunol 14, 1014–1022 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 20.Joyce JA & Fearon DT T cell exclusion, immune privilege, and the tumor microenvironment. Science 348, 74–80 (2015). [DOI] [PubMed] [Google Scholar]
  • 21.Spranger S & Gajewski TF Tumor-intrinsic oncogene pathways mediating immune avoidance. Oncoimmunology 5, e1086862 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 22.Wherry EJ & Kurachi M Molecular and cellular insights into T cell exhaustion. Nat. Rev. Immunol 15, 486–499 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 23.Tauriello DVF et al. TGFβ drives immune evasion in genetically reconstituted colon cancer metastasis. Nature 554, 538–543 (2018). [DOI] [PubMed] [Google Scholar]
  • 24.Mariathasan S et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–548 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 25.Freedman D Statistical Models: Theory and Practice. (Cambridge Univ. Press, Cambridge, 2009). [Google Scholar]
  • 26.Kleinbaum DG Survival analysis, a self-learning text. Biom. J 40, 107–108 (1998). [Google Scholar]
  • 27.Patel SJ et al. Identification of essential genes for cancer immunotherapy. Nature 548, 537–542 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 28.Khong HT & Rosenberg SA The Waardenburg syndrome type 4 gene, SOX10, is a novel tumor-associated antigen identified in a patient with a dramatic response to immunotherapy. Cancer Res 62, 3020–3023 (2002). [PMC free article] [PubMed] [Google Scholar]
  • 29.Thomas DA & Massague J TGF-β directly targets cytotoxic T cell functions during tumor evasion of immune surveillance. Cancer Cell 8, 369–380 (2005). [DOI] [PubMed] [Google Scholar]
  • 30.Woo EY et al. Regulatory CD4+CD25+ T cells in tumors from patients with early-stage non-small cell lung cancer and late-stage ovarian cancer. Cancer Res 61, 4766–4772 (2001). [PubMed] [Google Scholar]
  • 31.The Cancer Genome Atlas Research Network et al. The Cancer Genome Atlas Pan-Cancer analysis project. Nat. Genet 45, 1113–1120 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 32.Curtis C et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature 486, 346–352 (2012). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 33.Zhou P et al. In vivo discovery of immunotherapy targets in the tumour microenvironment. Nature 506, 52–57 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 34.Giordano M et al. Molecular profiling of CD8 T cells in autochthonous melanoma identifies Maf as driver of exhaustion. EMBO J. 34, 2042–2058 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 35.Wakamatsu E, Mathis D & Benoist C Convergent and divergent effects of costimulatory molecules in conventional and regulatory CD4+ T cells. Proc. Natl Acad. Sci. USA 110, 1023–1028 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 36.Twyman-Saint Victor C et al. Radiation and dual checkpoint blockade activate non-redundant immune mechanisms in cancer. Nature 520, 373–377 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 37.Schietinger A et al. Tumor-specific T cell dysfunction is a dynamic antigen-driven differentiation program initiated early during tumorigenesis. Immunity 45, 389–401 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 38.Philip M et al. Chromatin states define tumour-specific T cell dysfunction and reprogramming. Nature 545, 452–456 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 39.Pollizzi KN et al. mTORC1 and mTORC2 selectively regulate CD8+ T cell differentiation. J. Clin. Investig 125, 2090–2108 (2015). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 40.Huse M, Lillemeier BE, Kuhns MS, Chen DS & Davis MM T cells use two directionally distinct pathways for cytokine secretion. Nat. Immunol 7, 247–255 (2006). [DOI] [PubMed] [Google Scholar]
  • 41.Frauwirth KA et al. The CD28 signaling pathway regulates glucose metabolism. Immunity 16, 769–777 (2002). [DOI] [PubMed] [Google Scholar]
  • 42.Barrett T et al. NCBI GEO: archive for functional genomics data sets—update. Nucleic Acids Res 41, D991–D995 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 43.Remark R et al. Characteristics and clinical impacts of the immune environments in colorectal and renal cell carcinoma lung metastases: influence of tumor origin. Clin. Cancer Res 19, 4079–4091 (2013). [DOI] [PubMed] [Google Scholar]
  • 44.Garcia-Diaz A et al. Interferon receptor signaling pathways regulating PD-L1 and PD-L2 expression. Cell Rep 19, 1189–1201 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 45.Riaz N et al. Tumor and microenvironment evolution during immunotherapy with Nivolumab. Cell 171, 934–949 e915 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 46.Charoentong P et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint nlockade. Cell Rep 18, 248–262 (2017). [DOI] [PubMed] [Google Scholar]
  • 47.Wolchok JD et al. Guidelines for the evaluation of immune therapy activity in solid tumors: immune-related response criteria. Clin. Cancer Res 15, 7412–7420 (2009). [DOI] [PubMed] [Google Scholar]
  • 48.Nathanson T et al. Somatic mutations and neoepitope homology in melanomas treated with CTLA-4 blockade. Cancer Immunol. Res 5, 84–91 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 49.Kaiserman D & Bird PI Control of granzymes by serpins. Cell Death Differ 17, 586–595 (2010). [DOI] [PubMed] [Google Scholar]
  • 50.Hirst CE et al. The intracellular granzyme B inhibitor, proteinase inhibitor 9, is up-regulated during accessory cell maturation and effector cell degranulation, and its overexpression enhances CTL potency. J. Immunol 170, 805–815 (2003). [DOI] [PubMed] [Google Scholar]
  • 51.Bladergroen BA et al. The granzyme B inhibitor, protease inhibitor 9, is mainly expressed by dendritic cells and at immune-privileged sites. J. Immunol 166, 3218–3225 (2001). [DOI] [PubMed] [Google Scholar]
  • 52.Hirst CE et al. Perforin-independent expression of granzyme B and proteinase inhibitor 9 in human testis and placenta suggests a role for granzyme B-mediated proteolysis in reproduction. Mol. Hum. Reprod 7, 1133–1142 (2001). [DOI] [PubMed] [Google Scholar]
  • 53.Medema JP et al. Blockade of the granzyme B/perforin pathway through overexpression of the serine protease inhibitor PI-9/SPI-6 constitutes a mechanism for immune escape by tumors. Proc. Natl Acad. Sci. USA 98, 11515–11520 (2001). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 54.Uhlen M et al. A pathology atlas of the human cancer transcriptome. Science 357, eaan2507 (2017). [DOI] [PubMed] [Google Scholar]
  • 55.Schoenborn JR & Wilson CB Regulation of interferon-gamma during innate and adaptive immune responses. Adv. Immunol 96, 41–101 (2007). [DOI] [PubMed] [Google Scholar]
  • 56.Taniguchi T, Ogasawara K, Takaoka A & Tanaka N IRF family of transcription factors as regulators of host defense. Annu. Rev. Immunol 19, 623–655 (2001). [DOI] [PubMed] [Google Scholar]
  • 57.Papatheodorou I et al. Expression Atlas: gene and protein expression across multiple studies and organisms. Nucleic Acids Res 46, D246–D251 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 58.Fernandez-Banet J et al. OASIS: web-based platform for exploring cancer multi-omics data. Nat Methods 13, 9–10 (2016). [DOI] [PubMed] [Google Scholar]
  • 59.Hoadley KA et al. Multiplatform analysis of 12 cancer types reveals molecular classification within and across tissues of origin. Cell 158, 929–944 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 60.Miyan M, Schmidt-Mende J, Kiessling R, Poschke I & de Boniface J Differential tumor infiltration by T-cells characterizes intrinsic molecular subtypes in breast cancer. J. Transl. Med 14, 227 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 61.Benjamini Y & Hochberg Y Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300 (1995). [Google Scholar]
  • 62.Storey JD & Tibshirani R Statistical significance for genomewide studies. Proc. Natl Acad. Sci. USA 100, 9440–9445 (2003). [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 63.Pan D et al. A major chromatin regulator determines resistance of tumor cells to T cell-mediated killing. Science 359, 770–775 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

1
NIHMS1024718-supplement-1.pdf (9.4MB, pdf)

Data Availability Statement

Users can query our analysis results with gene names: http://tide.dfci.harvard.edu/query/. All of our processed input data, analysis output data and an example script to repeat our major results are available at http://tide.dfci.harvard.edu/download/.

RESOURCES