Abstract
Purpose: Throughout their development, tumors are challenged by the immune system, and they acquire features to evade its surveillance. A systematic view of these traits, which shed light on how tumors respond to immunotherapies, is still lacking.
Experimental Design: Here, we computed the relative abundance of an array of immune cell populations to measure the immune infiltration pattern of 9,174 tumors of 29 solid cancers. We then clustered tumors with similar infiltration pattern to define immunophenotypes. Finally, we identified genomic and transcriptomic traits associated to these immunophenotypes across cancer types.
Results: In highly cytotoxic immunophenotypes, we found tumors with low clonal heterogeneity enriched for alterations of genes involved in epigenetic regulation, ubiquitin-mediated proteolysis, antigen presentation, and cell–cell communication, which may drive resistance in combination with the ectopic expression of negative immune checkpoints. Tumors with immunophenotypes of intermediate cytotoxicity are characterized by an upregulation of processes involved in neighboring tissue invasion and remodeling that may foster the recruitment of immunosuppressive cells. Tumors with poorly cytotoxic immunophenotype tend to be of more advanced stages and bear a greater burden of copy number alterations and frequent alterations of cell cycle, hedgehog, β-catenin, and TGFβ pathways, which may cause immune depletion.
Conclusions: We provide a comprehensive landscape of the characteristics of solid tumors that may influence (or be influenced by) the characteristics of their immune infiltrate. These results may help interpret the response of solid tumors to immunotherapies and guide the development of novel drug combination strategies. Clin Cancer Res; 24(15); 3717–28. ©2018 AACR.
In this study, we define the immunophenotype of 9,174 tumors of 29 solid cancers on the basis of a set of 16 immune cell populations infiltrating them. We identify genomic and transcriptomic features of individual tumors significantly associated with these immunophenotypes. These features will help to shed light on how tumors respond to immunotherapies, as well as to provide rationale for the development of novel therapeutic strategies. Oncogenic traits associated to immunophenotypes with high cytotoxicity may reveal mechanisms of immune resistance that render them unresponsive to checkpoint inhibitors. On the other hand, the features associated to immunophenotypes of lower cytotoxicity could be targeted by existing drugs, potentially improving the outcome of immunotherapies. We provide examples of these translational opportunities in the article.
Introduction
Solid tumors and the immune cells infiltrating them interact in a dynamic equilibrium that shapes the progression of the disease. Tumor evasion of immunosurveillance is a hallmark shared by all types of cancer (1). In this process, the features that protect cancer cells from the action of a cytotoxic immune infiltrate or promote the suppression of this infiltrate are positively selected. Although several mechanisms of immune evasion have been experimentally identified for some tumor types (2–4), to date, there is no comprehensive view of the routes of tumor escape. Some of these mechanisms, such as the activation of immune checkpoints, have been therapeutically exploited, with a remarkable clinical success across a variety of cancer types (4). Therefore, the discovery of features acquired by tumors in response to the immune cells that surround them may open up new strategies to treat the disease.
The genomic and transcriptomic sequences of large cohorts of tumors generated by international projects, such as The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/), have provided the opportunity to identify the distinguishing molecular traits of cancer in unprecedented detail. Furthermore, the abundance of immune cells in the infiltrate of tumors can be estimated using computational methods that exploit the fact that tumor bulk samples are admixtures of cancer cells and cells of their microenvironment (5). Although several recent studies using these approaches have explored some features of the tumors associated with specific characteristics of the immune infiltrate (6–10), a comprehensive pan-cancer landscape of the interactions between the tumor and the immune cells is still lacking.
Here, we first aimed to determine to which extent the cancer type and the tissue of origin shape the immune infiltrate of tumors. To this end, we estimated the immune infiltration pattern of 9,174 tumors and computed its heterogeneity across the 29 solid cancer types they represent. We also compared this immune infiltration pattern to that of 6,544 healthy donors. Our results showed that the cancer type and the tissue of origin do not explain all the variability of the immune infiltrate of tumors. Therefore, we set out to identify intrinsic tumor features that significantly associate with their immune infiltrate. We first grouped the tumors of each cancer type with similar immune infiltration pattern (immunophenotypes). These immunophenotypes represent differentiated scenarios of immune infiltration, and therefore, we reasoned that they could result in the selection of different tumor mechanisms for evading the action of the immune system. To test this hypothesis, we systematically correlated the immunophenotypes with: (i) clinical and pathologic features of the tumors; (ii) genomic characteristics, including putative driver mutations and copy number alterations (CNA; identified by the Cancer Genome Interpreter; ref. 11); and (iii) the activation status of cancer cell processes (after subtracting the contribution of the immune cells to the expression of genes in the bulk sample; ref. 12). As a result, we revealed that solid cancers develop a very distinct set of molecular traits when progressing in different scenarios of immune infiltration. Our findings provide a landscape of the interactions between tumor and immune cells across cancer types and have clear implications for anticancer immunotherapies currently employed in the clinic and their potential combination with targeted drugs.
Materials and Methods
Sample data collection and processing
TCGA data, mutations, CNAs, gene expression (RSEM gene-normalized), clinical annotations (with tumor stage manually annotated using the available pathologic and clinical data), for 9,174 tumors across 29 cancer types were downloaded from the most recent freeze contained in the Firebrowse server (2016_01_28). Triple-negative breast tumors were considered as a separated cancer cohort. The expression data (RPKM) of 8,304 donors across 22 healthy tissues were downloaded from the GTEx portal (version 6). The quantification of virus expression across tumors was retrieved from previous publications (7, 13). We considered a tumor infected if the expression of a virus was higher than the observed across healthy tissues, following the approach described by Rooney and colleagues (7).
A tumor was considered hypermutated if its mutation burden was greater than the median mutation burden of its cohort plus 4.5 times the interquartile range and greater than 500 mutations. The clonal heterogeneity of each tumor was computed as the ratio between the dispersion of the variant allele fraction of the mutations in the sample and its median variant allele fraction (9). We analyzed the RNA sequencing (RNA-seq) data of samples from melanoma patients treated with anti-CTLA4 (14) (RPKM matrix provided by the authors, 42 patients) and anti–PD-1 (FPKM matrix downloaded from GEO: GSE78220, 28 patients; ref. 15) therapies.
Estimation of the abundance of immune cell populations
The relative abundance of 16 immune populations in tumors and healthy tissues were computed from the RNA-seq of each bulk sample. In detail, we used the gene set variation analysis (GSVA; ref. 16), a sample-level enrichment method, which proved to be more appropriate for our aim than available deconvolution-based methods (Supplementary Material). We selected immune cell populations, the relative abundance of which could be computed with confidence from the expression data of the admixture samples (Supplementary Material). Immune populations were represented by nonoverlapping gene sets specifically overexpressed in each cell type, selected from several previous publications (5, 17, 18) (Supplementary Material; Supplementary Table S1A). We also calculated the relative abundance of cytotoxic cells as a proxy of the antitumor activity of the overall immune infiltrate. The gene set to compute the abundance of cytotoxic cells includes genes overexpressed in CD8-activated T cells, γδ T cells, and natural killer cells (Supplementary Material). The GSVA produces normalized enrichment scores ranging from −1 to 1, which represent the abundance of the immune cell population in the sample relative to other tumors of the analyzed cohort. Overall, 31 independent GSVA analyses were carried out for the pan-cancer cohort (all tumor samples pooled together, n = 1), the healthy donors cohort (all healthy samples pooled together, n = 1), and each cancer cohort separately (n = 29).
Comparison of the immune infiltration pattern of tumors and their tissue of origin
We first matched each solid cancer type in TCGA to its closest corresponding healthy tissue included in the GTEx project (Supplementary Material). We then compared the GSVA scores computed for the pan-cancer cohort with those in matching healthy tissue samples (Supplementary Material). We considered that the abundance of a given immune population was different between tumors and healthy samples if the median of GSVA values across tumors of the cancer type differed significantly (Mann–Whitney Q value < 0.1) and more than 0.2 from that measured across matching healthy samples.
Tumor sample clustering
Tumors with qualitatively different immune infiltration pattern were grouped using a hierarchical agglomerative clustering (based on Euclidean distance and Ward linkage) method. To cluster the tumors across the whole pan-cancer cohort, we used the GSVA scores of 16 immune cell populations of the 9,174 tumors analyzed together. To build the immunophenotypes, we clustered the tumors of each cohort separately. We included the GSVA scores of cytotoxic cells and overweighed this signature by a factor of four in the process of clustering (Supplementary Material). The optimal number of clusters of tumors was determined on the basis of the percentage of variance of the data explained by them (Supplementary Material).
Quantification of mutational signatures
The deconstructSigs software (19) was used with default parameters to compute the contribution of previously reported mutational signatures (http://cancer.sanger.ac.uk/cosmic/signatures) to the overall landscape of mutations of each tumor. This calculation was limited to tumor samples with at least 50 mutations. We considered all mutations fitting signatures 3, 6, 10, 15, 18, 20, 21, and 26 (among the 30 available in the current version of http://cancer.sanger.ac.uk/cosmic/signatures) likely caused by altered POLE and BRCA-1/2, and base excision repair and mismatch repair deficiency; and those fitting signatures 2 and 13 caused by APOBEC activity.
Identification of genomic drivers associated to immunophenotypes
We first identified genomic alterations (point mutations, small indels, and gene copy number gains or losses) that are potential drivers of each tumor using the Cancer Genome Interpreter (https://www.cancergenomeinterpreter.org; ref. 11). Then, to identify the genes with driver alterations that significantly associate with a particular immunophenotype, we implemented a regularized logistic regression analysis adjusted by the count of protein-affecting mutations and gene CNAs. Our model specifically addresses the spurious fitting derived from the sparsity of the data and provides empirical significance scores (Supplementary Material). The choice of the regularization algorithm can be interpreted as a Bayesian logistic regression with a Gaussian weakly-informative prior distribution for the parameters of the immune clusters and adjustment covariates. We applied the logistic regression to each cancer cohort separately, and we also performed a pan-cancer analysis in which all the tumors with the equivalent immunophenotype (the ordinal number representing the cytotoxicity level) were pooled together irrespective of their cancer type. Only genes bearing mutations in at least 2.5% of the samples of the cohort and/or CNAs in at least 5% were included.
Expression adjustment
To adjust the gene expression measured in a tumor bulk sample by its immune content, we followed the rationale described in ref. 12. Briefly, we adjusted the expression value of each gene in each tumor sample according to the contribution of CD45 to the expression of that gene observed in the matching healthy tissue (Supplementary Material; Supplementary Fig. S1). Cancers with no matching healthy tissue data available (cholangiocarcinoma, mesothelioma, and uveal melanomas) were excluded from this analysis.
Pathway enrichment analysis
We identified upregulated pathways among tumors of a certain immunophenotype running a gene set enrichment analysis (GSEA; ref. 20) on the adjusted expression data (available for 8,971 tumors across 26 cancer types, due to the exclusion of the three cohorts in which the adjustment of the expression for the leukocyte content was not available; see above). Significance was considered for values of corrected P (which takes into account the size of each gene set plus the multiple test correction, as provided by the method) lower than 0.25, as recommended by the authors (http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html?Interpreting_GSEA). Gene sets were downloaded from the MSigDB database (20). We included broad hallmarks and specific pathways of interest from the curated gene sets/canonical pathways collection. In addition, we manually gathered pathways not found in the aforementioned collections. Overall, we considered a total of 51 pathways (Supplementary Table S1B) with minimal gene set overlap between them (Supplementary Fig. S2A) as well as with the gene sets used to compute the immune cells abundance (Supplementary Fig. S2B).
General model of association between tumor features and immunophenotypes
We grouped the genomic and transcriptomic features of the tumors into several broad terms that summarize our results. For example, the transcriptomic term “invasion and remodeling” comprises the findings for the epithelial-to-mesenchymal transition (MsigDB Hallmarks), extracellular matrix organization (Reactome database), and focal adhesion (KEGG database) pathways. The details of the mapping of pathways and genes to broad terms employed in the model are included in Supplementary Table S2. To build the corresponding summary figure, we considered that each of these broad terms was associated to the highly cytotoxic immune scenario if the median of the immunophenotypes (numbered from 1 to 6 from lowest to highest relative abundance of cytotoxic cells) associated to all the features integrating the term was greater than 4.5, to the poorly cytotoxic immune scenario if the median was lower than 2.5, and to the mid cytotoxic immune scenario otherwise.
Statistical analysis
Unless explicitly stated, the association between pairs of categorical variables was evaluated using the Fisher exact test. The distributions of two sets of any continuous variable were compared using the Mann–Whitney U test. The homogeneity of the distribution of samples across different groups was computed via an entropy score (described in Supplementary Methods). The influence of the immunophenotype on survival was evaluated through the Cox proportional hazard model. A linear regression was used to assess whether the value of a variable (e.g., the mutation burden) increases (or decreases) across immunophenotypes. For the Cox and lineal regression analyses, a model selection procedure was conducted by considering several cluster aggregation levels for the immune phenotypes, resulting in a model with a three-level factor. The selection criteria were the size and significance of the parameter estimates and the likelihood of the model (Akaike information criterion). P values were corrected using the two-stage Benjamini–Hochberg FDR method as appropriate and a corrected P value of 0.1 was considered statistically significant unless otherwise stated.
Results
The immune infiltration pattern of tumors is heterogeneous across and within cancer types
We defined the immune infiltration pattern of each tumor as the relative abundance of an array of 16 cell populations of the adaptive and innate immune system. These populations and the gene signatures representing them were selected upon comparison of several sources (Supplementary Methods). In each case, we chose the most specific cell population that produced biologically sound results (Supplementary Methods; Supplementary Table S1). Also, after rigorous evaluation of the performance of several tools, we employed a sample-level gene set enrichment method (GSVA; ref. 16) to compute the relative abundance of selected cell populations from the tumor RNA-seq bulk data (Supplementary Methods).
We first computed the immune infiltration pattern, that is, the GSVA enrichment scores of the 16 selected immune cell populations, across the 9,174 tumors (pan-cancer cohort; Table 1; Supplementary Table S3A). This analysis revealed that the relative abundance of many immune cell populations is correlated (Supplementary Fig. S3), suggesting at least some degree of coinfiltration of the tumor (21). We then computed the healthy immune infiltration pattern of samples of 22 matching healthy tissues from 6,544 donors using the same method (Supplementary Table S3B) and observed that the immune infiltration pattern of tumors clearly deviated from that of their matched healthy tissue (Supplementary Fig. S4). Next, we grouped the GSVA scores of the 9,174 tumors using a hierarchical clustering. We obtained 17 groups (Supplementary Fig. S5, see Materials and Methods), which reflect distinct profiles of immune infiltration of solid tumors (Supplementary Fig. S6). Tumors of some cancer types (e.g., prostate adenocarcinomas and malignancies affecting the immune privileged brain and eye; refs. 22, 23) appeared mostly in a single cluster. However, tumors of most malignancies were distributed across several clusters with varying degrees of dispersion. Furthermore, tumors of different cancer types grouped together, suggesting that they possess similar immune infiltration patterns. In summary, the tissue of origin and the cancer type of a tumor do not suffice to determine its immune infiltration pattern. The variability of this immune infiltrate could thus be largely explained by specific features of individual tumors.
Summary of the cohorts of tumors employed in the study
Cancer type acronym . | Cancer type full name . | Number of patients (RNA-seq) . |
---|---|---|
ACC | Adrenocortical carcinoma | 78 |
BLCA | Bladder urothelial carcinoma | 404 |
BRCA | Breast invasive carcinoma | 924 |
BRCA-TN | Breast triple-negative invasive carcinoma | 158 |
CESC | Cervical squamous cell carcinoma and endocervical adenocarcinoma | 301 |
CHOL | Cholangiocarcinoma | 36 |
COREAD | Colorectal adenocarcinoma | 370 |
ESCA | Esophageal carcinoma | 182 |
GBM | Glioblastoma multiforme | 163 |
HNSC | Head and neck squamous cell carcinoma | 515 |
KICH | Kidney chromophobe | 65 |
KIRC | Kidney renal clear cell carcinoma | 515 |
KIRP | Kidney renal papillary cell carcinoma | 285 |
LGG | Brain lower grade glioma | 514 |
LIHC | Liver hepatocellular carcinoma | 368 |
LUAD | Lung adenocarcinoma | 511 |
LUSC | Lung squamous cell carcinoma | 485 |
MESO | Mesothelioma | 87 |
OV | Ovarian serous cystadenocarcinoma | 300 |
PAAD | Pancreatic adenocarcinoma | 156 |
PCPG | Pheochromocytoma and paraganglioma | 178 |
PRAD | Prostate adenocarcinoma | 493 |
SARC | Sarcoma | 254 |
SKCM | Skin cutaneous melanoma | 434 |
STAD | Stomach adenocarcinoma | 405 |
THCA | Thyroid carcinoma | 499 |
UCEC | Uterine corpus endometrial carcinoma | 357 |
UCS | Uterine carcinosarcoma | 57 |
UVM | Uveal melanoma | 80 |
Total | 9,174 |
Cancer type acronym . | Cancer type full name . | Number of patients (RNA-seq) . |
---|---|---|
ACC | Adrenocortical carcinoma | 78 |
BLCA | Bladder urothelial carcinoma | 404 |
BRCA | Breast invasive carcinoma | 924 |
BRCA-TN | Breast triple-negative invasive carcinoma | 158 |
CESC | Cervical squamous cell carcinoma and endocervical adenocarcinoma | 301 |
CHOL | Cholangiocarcinoma | 36 |
COREAD | Colorectal adenocarcinoma | 370 |
ESCA | Esophageal carcinoma | 182 |
GBM | Glioblastoma multiforme | 163 |
HNSC | Head and neck squamous cell carcinoma | 515 |
KICH | Kidney chromophobe | 65 |
KIRC | Kidney renal clear cell carcinoma | 515 |
KIRP | Kidney renal papillary cell carcinoma | 285 |
LGG | Brain lower grade glioma | 514 |
LIHC | Liver hepatocellular carcinoma | 368 |
LUAD | Lung adenocarcinoma | 511 |
LUSC | Lung squamous cell carcinoma | 485 |
MESO | Mesothelioma | 87 |
OV | Ovarian serous cystadenocarcinoma | 300 |
PAAD | Pancreatic adenocarcinoma | 156 |
PCPG | Pheochromocytoma and paraganglioma | 178 |
PRAD | Prostate adenocarcinoma | 493 |
SARC | Sarcoma | 254 |
SKCM | Skin cutaneous melanoma | 434 |
STAD | Stomach adenocarcinoma | 405 |
THCA | Thyroid carcinoma | 499 |
UCEC | Uterine corpus endometrial carcinoma | 357 |
UCS | Uterine carcinosarcoma | 57 |
UVM | Uveal melanoma | 80 |
Total | 9,174 |
Identification of six immunophenotypes across cancer types
We reasoned that the cytotoxic component of the immune infiltrate is the main responsible for the selective pressure applied on tumors, which may result in the positive selection of specific features that support the escape of immunosurveillance. Therefore, we incorporated an additional gene signature representing the meta-population of all cytotoxic cells (see Materials and Methods). We then recomputed the immune infiltration pattern of tumors within each cancer cohort separately, to get rid of the confounder of cancer type-specific patterns. Intracohort GSVA scores were then used to group tumors via a hierarchical clustering (Fig. 1A; Supplementary Tables S3C–S3ZE) in which the GSVA score of the cytotoxic cells received a 25% of the overall weight (Supplementary Methods). As a result, we obtained six (see Materials and Methods) cytotoxic-driven clusters (Fig. 1B and C; Supplementary Fig. S7; Supplementary Table S4), which we named immunophenotypes. The weight of the cytotoxic meta-population in the clustering was selected to guarantee that immunophenotypes had different median relative abundance of cytotoxic cells, but that also the abundance of the remaining 16 cell populations contributed to the clustering (see Materials and Methods). As a result, although the immunophenotypes from one to six contained growing abundances of cytotoxic cells, the tumors contained in each of them differ between 25% and 54%, depending on the cancer cohort, with respect to clusters built solely on the basis of the abundance of cytotoxic cells (cyt-clusters; Supplementary Fig. S8). Moreover, the abundance of most cell populations grows monotonically along cyt-clusters, probably due to the aforementioned coinfiltration effect. On the contrary, the immunophenotypes group tumors whose abundance of certain immune populations deviate from the monotonic trend of the increase of cytotoxic cells (Supplementary Fig. S9). We hypothesized that these patterns captured by the immunophenotypes may relate to particular traits of the tumors that would be missed by cyt-clusters.
Tumor immunophenotypes. A, We defined the immune infiltration pattern of a tumor as the relative abundance of a set of immune cell populations estimated via a GSVA. The tumors of each TCGA cohort were then grouped into six immunophenotypes with different cytotoxic content. B, Immunophenotypes of the breast cancer cohort (n = 924) excluding triple-negative breast tumors. Immunophenotypes are numbered in ascending order of their median relative abundance of cytotoxic cells (from one to six and from dark blue to dark red in the bar below the heatmap). C, Proportion of tumors with each immunophenotype across cohorts. Note that immunophenotypes represent tumors with distinct immune infiltrate profiles within each cohort regardless of the overall immune infiltration differences between malignancies. Cancer types are displayed in ascending order of their absolute leukocyte content (measured through CD45 expression of the bulk samples; left boxplot). D, Influence of the immunophenotypes on the overall survival. A log2 HR below zero (vertical dot line) indicates that patients with tumors of higher cytotoxic immunophenotypes exhibit improved survival. P values were calculated using the Cox regression model. E, Association between advanced (III/IV) pathologic stage of the tumors and their immunophenotypes at diagnosis. Events associated with higher cytotoxic immunophenotypes are shown in red, and events associated with lower cytotoxic immunophenotypes are in blue (solid circles for Q values <0.1, empty circles for Q values <0.2, calculated using linear regression; see Materials and Methods).
Tumor immunophenotypes. A, We defined the immune infiltration pattern of a tumor as the relative abundance of a set of immune cell populations estimated via a GSVA. The tumors of each TCGA cohort were then grouped into six immunophenotypes with different cytotoxic content. B, Immunophenotypes of the breast cancer cohort (n = 924) excluding triple-negative breast tumors. Immunophenotypes are numbered in ascending order of their median relative abundance of cytotoxic cells (from one to six and from dark blue to dark red in the bar below the heatmap). C, Proportion of tumors with each immunophenotype across cohorts. Note that immunophenotypes represent tumors with distinct immune infiltrate profiles within each cohort regardless of the overall immune infiltration differences between malignancies. Cancer types are displayed in ascending order of their absolute leukocyte content (measured through CD45 expression of the bulk samples; left boxplot). D, Influence of the immunophenotypes on the overall survival. A log2 HR below zero (vertical dot line) indicates that patients with tumors of higher cytotoxic immunophenotypes exhibit improved survival. P values were calculated using the Cox regression model. E, Association between advanced (III/IV) pathologic stage of the tumors and their immunophenotypes at diagnosis. Events associated with higher cytotoxic immunophenotypes are shown in red, and events associated with lower cytotoxic immunophenotypes are in blue (solid circles for Q values <0.1, empty circles for Q values <0.2, calculated using linear regression; see Materials and Methods).
We next explored whether the different immunophenotypes associated with clinical and pathologic characteristics of the tumors. In coherence with previous studies (8, 24), we observed that the survival of patients tended to improve with the increase of the cytotoxic level of the immunophenotypes in several cancer types (Fig. 1D). An intriguing exception was low-grade gliomas, the clinical outcome of which is worse in the case of immunophenotypes of higher cytotoxicity. This may be explained by the loss of integrity of the blood–brain barrier, which facilitates immune infiltration but also correlates with high tumor size and aggressiveness (Supplementary Fig. S10; refs. 25, 26). Finally, we also observed that tumors of more advanced pathologic stage at diagnosis significantly associated to lower cytotoxicity of the immunophenotypes in several cancer types (Fig. 1E). These observations suggest that tumors preferentially progress in the presence of a poorly cytotoxic immune infiltrate. In contrast, tumors with highly cytotoxic immunophenotype would be partially kept in check by the immune system and progress less frequently to more advanced stages.
Different tumor genomic events associated to immunophenotypes
We reasoned that different selective pressures represented by diverse immunophenotypes could result in the positive selection of distinct tumoral mechanisms of immune evasion. Therefore, we next aimed to identify the genomic features of tumors associated with these immunophenotypes.
The mutational burden of tumors, treated as a continuous variable, positively correlated with the level of cytotoxicity of their immunophenotype only in few cancer types (Fig. 2A). When we specifically evaluated hypermutated tumors and/or tumors with mutational signatures reflective of DNA repair deficiencies (Materials and Methods), we found them to be significantly overrepresented among immunophenotypes of higher cytotoxicity in four cancer types (Supplementary Figs. S11 and S12). The presence of viral infection, another known immunogenic event, inferred from the abundance of likely APOBEC-caused mutations (27) or the expression of known viral genes (7, 13), positively correlated with the level of cytotoxicity of the immunophenotypes in six cancer types (Fig. 2B). Strikingly, hepatitis B–infected hepatocellular carcinomas were significantly overrepresented among poorly cytotoxic immunophenotypes. On the other hand, in agreement with a recent study (28), we observed that the burden of gene CNAs of tumors negatively correlated with the cytotoxicity of their immunophenotypes in a large fraction of cancer types (Fig. 2C). Finally, we observed that tumors with a more heterogeneous clonal composition were associated with lowly cytotoxic immunophenotypes in a significant number of solid cancers (Fig. 2D). This may be the result of tumor immune-edition, in which the action of a highly cytotoxic immune infiltrate results in the restriction of the clonal diversity of tumors.
Genomic features of tumors that significantly associate to immunophenotypes. Association between the mutational (missense and frameshift) burden (A), virus infection (B), the CNAs burden (C), and the clonal heterogeneity (D) of tumors and their immunophenotypes. Events associated to highly cytotoxic immunophenotypes are highlighted in red, and events associated with lowly cytotoxic immunophenotypes are in blue (solid circles for Q values <0.1, empty circles for Q values <0.2, calculated using linear regression; see Materials and Methods). E, Genes with driver alterations significantly associated to immunophenotypes (see Materials and Methods). Associations with Q value <0.1 for individual cancer cohorts or for the pan-cancer cohort were considered significant. Associations with uncorrected P value <0.05 for specific cancer cohorts that appeared significant in the pan-cancer analysis are also shown (HLA-A, -B, and -C mutations are considered together). The color of the circle represents the cancer immunophenotype significantly enriched by driver alterations of the gene (see legend at the bottom), and the size of the circle is proportional to the magnitude of that enrichment. The type of driver alterations (mutations, amplifications, or deletions) is color-coded at the left of the gene name. Cancer types are sorted in decreasing order of number of genes with significant associations (at least one).
Genomic features of tumors that significantly associate to immunophenotypes. Association between the mutational (missense and frameshift) burden (A), virus infection (B), the CNAs burden (C), and the clonal heterogeneity (D) of tumors and their immunophenotypes. Events associated to highly cytotoxic immunophenotypes are highlighted in red, and events associated with lowly cytotoxic immunophenotypes are in blue (solid circles for Q values <0.1, empty circles for Q values <0.2, calculated using linear regression; see Materials and Methods). E, Genes with driver alterations significantly associated to immunophenotypes (see Materials and Methods). Associations with Q value <0.1 for individual cancer cohorts or for the pan-cancer cohort were considered significant. Associations with uncorrected P value <0.05 for specific cancer cohorts that appeared significant in the pan-cancer analysis are also shown (HLA-A, -B, and -C mutations are considered together). The color of the circle represents the cancer immunophenotype significantly enriched by driver alterations of the gene (see legend at the bottom), and the size of the circle is proportional to the magnitude of that enrichment. The type of driver alterations (mutations, amplifications, or deletions) is color-coded at the left of the gene name. Cancer types are sorted in decreasing order of number of genes with significant associations (at least one).
Then, we focused on individual driver genomic alterations identified in each of the tumors across the 29 cohorts (see Materials and Methods) using the Cancer Genome Interpreter (11). To address the statistical issues derived from the sparsity of the data, we implemented a regularized logistic regression analysis adjusting for the number of mutations and CNAs of each tumor (see Materials and Methods). The regression was computed in each cohort separately and also grouping all tumors with equivalent immunophenotypes from different cohorts, thus increasing the statistical power to detect associations that are coherent across malignancies (Fig. 2E; Supplementary Table S5). Some driver alterations that were significantly overrepresented among tumors with highly cytotoxic immunophenotypes (immunophenotypes 5 and 6) have already been described to elicit resistance to immune destruction. These include mutations of members of the HLA gene family and B2M (part of the antigen presentation machinery), mutations of CASP8 (extrinsic apoptosis pathway), and amplifications of the PDL-1/2 genes (negative immune checkpoints; refs. 3, 4, 7, 29). We also identified other genomic alterations significantly overrepresented among these tumors affecting epigenetic regulators (e.g., EP300, ARID1A, ARID2, TRRAP, and CHD8), E3 ubiquitin ligases (e.g., UBR5, CUL3, and TRIP12), and genes encoding cell–cell interaction proteins (e.g. MYH9, MYH10, ANK3, ROCK12, CDH1, and FAT1). These alterations constitute novel potential drivers of immune resistance. On the other hand, genes involved in the regulation of cell cycle, DNA replication, and telomere maintenance (e.g., ATRX, STAG2, AURKA, BAP1, FBXW7, SMARCB1, CDKN2A, CDKN2B, CCND2, CCND3, and CDK6), and others part of the WNT–β-catenin signaling pathway (e.g., CTNNB1, APC, AXIN1, and SMAD4) contained driver alterations that were significantly overrepresented among tumors with poorly cytotoxic immunophenotypes (immunophenotypes 1 and 2). Of note, driver alterations affecting one gene were in some cases associated with different immunophenotypes in diverse cancer types. Mutations of NRAS/HRAS, for example, significantly associated to highly cytotoxic immunophenotypes in head and neck squamous carcinomas and paragangliomas but to poorly cytotoxic immunophenotypes in thyroid carcinomas. NF1 mutations were significantly associated with highly cytotoxic immunophenotypes in ovarian serous cystadenocarcinomas, low-grade gliomas, and glioblastomas, but with poorly cytotoxic immunophenotypes in paragangliomas. TP53 alterations, which significantly associated with various immunophenotypes across 15 malignancies, constituted the most salient example of this phenomenon. Of note, no consistent association between these driver alterations and particular immune cell populations could explain these seemingly contradictory observations across cancers (Supplementary Table S6).
In summary, we identified genomic events strongly associated with specific tumor immunophenotypes, which constitute a catalog of putative genomic drivers of immune evasion. Some of these associations varied depending on the cancer type, which could be explained by context-dependent effects of a gene alteration across malignancies.
Several transcriptional programs appear upregulated across immunophenotypes
To complete the landscape of the mechanisms of immune evasion, we identified cellular programs that appear significantly upregulated in different immunophenotypes. To do this, we first selected a comprehensive collection of cell pathways involved in cancer development with minimal overlap of their gene sets (see Materials and Methods; Supplementary Fig. S2; Supplementary Table S1B). Then, we inferred their differential activation across immunophenotypes from the collective upregulation of the genes integrating them using GSEA (20). Because tumor samples are admixtures of tumor cells and their microenvironment, we applied the GSEA after adjusting the expression level of each gene for the contribution of the immune content (see Materials and Methods). This step prevented the overestimation of upregulation of pathways with genes that are highly expressed in immune cells (e.g., immune-related cytokines or immune checkpoints) or the underestimation of those with genes with lower expression in nontumor cells (e.g., cancer testis antigens; Fig. 3A; Supplementary Fig. S1).
Immunophenotypes with enrichment for upregulated cell pathways across cohorts. A, Bar plot of the proportion of genes of each pathway whose expression changed after adjustment for the immune content of the tumor sample (pale green: decrease, dark green: increase, gray: unchanged; see Materials and Methods). B, Immunophenotypes enriched for selected upregulated cell pathways (Q value provided by GSEA <0.25; see Materials and Methods). Pathways are grouped in the panel according to broad cell processes. The position of each circle represents the median immunophenotype (between 1 and 6) in which the enrichment of the upregulated pathway was detected in each cancer type (black whiskers represent percentiles 25 and 75). The size of the circles follows the number of cancer types in which these significant associations are observed, with the shade of gray proportional to the median GSEA enrichment scores. C, Heatmap detailing the significance of the enrichment of the upregulation of each pathway for immunophenotypes across cohorts (which are summarized in B). Cohorts are sorted in decreasing order of the overall number of pathways with significant enrichment.
Immunophenotypes with enrichment for upregulated cell pathways across cohorts. A, Bar plot of the proportion of genes of each pathway whose expression changed after adjustment for the immune content of the tumor sample (pale green: decrease, dark green: increase, gray: unchanged; see Materials and Methods). B, Immunophenotypes enriched for selected upregulated cell pathways (Q value provided by GSEA <0.25; see Materials and Methods). Pathways are grouped in the panel according to broad cell processes. The position of each circle represents the median immunophenotype (between 1 and 6) in which the enrichment of the upregulated pathway was detected in each cancer type (black whiskers represent percentiles 25 and 75). The size of the circles follows the number of cancer types in which these significant associations are observed, with the shade of gray proportional to the median GSEA enrichment scores. C, Heatmap detailing the significance of the enrichment of the upregulation of each pathway for immunophenotypes across cohorts (which are summarized in B). Cohorts are sorted in decreasing order of the overall number of pathways with significant enrichment.
We observed that tumors with different immunophenotypes presented distinct sets of upregulated cell pathways (Fig. 3B and C; Supplementary Table S7). The upregulation of immune-related genes, such as the HLA I/II families, cancer germline antigens, and proinflammatory cytokines and chemokines, appeared overrepresented among tumors with highly cytotoxic immunophenotypes (5 and 6). These tumors were also enriched for the upregulation of several pathways of the energy metabolism, with the exception of glycolysis. Other pathways that appear upregulated in this scenario, such as the JAK/STAT signaling or a set of suppressive cytokines, may directly curb the pressure of a highly cytotoxic infiltrate (2–4). As expected, the overexpression of the PD-L1/PD-L2–negative immune checkpoints was also enriched among tumors with these immunophenotypes in most cancer types (Supplementary Fig. S13; refs. 4, 30).
The upregulation of pathways and processes related to tumor invasion and the remodeling of neighboring tissues, response to hypoxia, angiogenesis, glycolysis, focal adhesion, inflammation, epithelial-to-mesenchymal transition, and extracellular matrix remodeling, was overrepresented among tumors with immunophenotypes of intermediate cytotoxicity (3 and 4). Of note, paragangliomas and stomach adenocarcinomas did not follow this pattern, with these processes upregulated preferentially within highly cytotoxic immunophenotypes. As previously hypothesized, the discovery of these associations, related to immunophenotypes with high relative abundance of immunosuppressive cells, would have been missed if the tumors had been grouped only on the basis of their cytotoxic cells infiltrate (Supplementary Fig. S14).
Finally, lowly cytotoxic immunophenotypes (1 and 2) were enriched for the upregulation of TGFβ, Notch, and Hedgehog signaling pathways, which have been described to actively prevent immune infiltration (31–33). Furthermore, these immunophenotypes were also strongly associated with the overexpression of genes involved in the progression of the cell cycle, mRNA and protein synthesis, response to DNA damage, and telomere maintenance, all of them related to an accelerated rate of cell growth and proliferation. On the contrary, in head and neck and lung squamous carcinomas and colorectal adenocarcinomas, these pathways were upregulated primarily among tumors with highly or intermediate cytotoxic immunophenotypes.
Three scenarios of tumor immune infiltration
Our results identified immunophenotypes that represent distinct scenarios of immune infiltration in solid tumors. On the basis of these scenarios and the tumor features that we found significantly associated to each of them, we propose a general model that, despite variations between cancer types, describes the progression of tumors in these diverse microenvironments (Fig. 4).
A model of tumor development under different scenarios of immune infiltration. The left panel comprises tumor features associated to three scenarios of immune infiltration integrated, respectively, by the tumors of immunophenotypes 5–6 (high cytotoxicity; top), immunophenotypes 3–4 (intermediate cytotoxicity with large contribution of the immunosuppressive component; center), and immunophenotypes 1–2 (poor cytotoxicity; bottom). Each row of the figure groups several genomic or transcriptomic events (labeled as green squares and purple triangles, respectively) into broad terms (following the mapping in Supplementary Table S2). Circles mark the cancer types in which events included within the term are associated to one of these three scenarios of immune infiltration (see Materials and Methods). The figure only displays the scenario in which associations appear more frequent across cancer types, and thus, results in malignancies with different associations are not shown. The cohorts are sorted in descending order of total significant associations. Cohorts marked with an asterisk lack association with cell pathways due to absence of matching healthy tissue data. SHH, Hedgehog.
A model of tumor development under different scenarios of immune infiltration. The left panel comprises tumor features associated to three scenarios of immune infiltration integrated, respectively, by the tumors of immunophenotypes 5–6 (high cytotoxicity; top), immunophenotypes 3–4 (intermediate cytotoxicity with large contribution of the immunosuppressive component; center), and immunophenotypes 1–2 (poor cytotoxicity; bottom). Each row of the figure groups several genomic or transcriptomic events (labeled as green squares and purple triangles, respectively) into broad terms (following the mapping in Supplementary Table S2). Circles mark the cancer types in which events included within the term are associated to one of these three scenarios of immune infiltration (see Materials and Methods). The figure only displays the scenario in which associations appear more frequent across cancer types, and thus, results in malignancies with different associations are not shown. The cohorts are sorted in descending order of total significant associations. Cohorts marked with an asterisk lack association with cell pathways due to absence of matching healthy tissue data. SHH, Hedgehog.
The first scenario describes the progression of tumors with highly cytotoxic immunophenotypes (Fig. 4, top). This highly cytotoxic infiltrate may be partially explained by the enrichment of tumors for immunogenic events, such as the ectopic expression of cancer testis antigens, virus infection, or a high mutation burden frequently associated to alterations in genes involved in the intrinsic DNA damage response. Other potentially immunogenic events associated to this scenario include the upregulation of several oxygen free radical–producing metabolic pathways (34), and the upregulation of the antigen presentation machinery and IFN signaling. The selective pressure of this highly cytotoxic immunophenotype results in the recurrence, among these tumors, of loss-of-function mutations in genes already identified by previous publications, such as members of the HLA family, B2M and CASP8 (3, 4, 7, 29), and the development of other, also already described, mechanisms such as the upregulation of the PD-L1/2 genes (partially explained by their genomic amplification), inhibitory chemokines, and the JAK/STAT signaling pathway (35). Other genomic alterations overrepresented among these tumors, which had not been previously described, affect a variety of chromatin-regulatory factors, E3 ubiquitin ligases, and cell–cell interaction genes involved in several tumor processes that could potentially drive resistance to the effective immune infiltrate (36–41).
In the second scenario (Fig. 4, center), tumors with immunophenotypes of intermediate cytotoxicity are enriched for the upregulation of processes involved in the invasion and remodeling of neighboring tissues, such as the epithelial-to-mesenchymal transition, focal adhesion, and extracellular matrix remodeling. The upregulation of other hallmarks of tumor progression, such as angiogenesis, inflammation, and hypoxia, are also overrepresented among them. These processes have been associated with the infiltration of suppressive immune cell populations (42–45), such as macrophages that appear with abnormally high relative abundance within the immune infiltration pattern of these tumors (Supplementary Fig. S9). Of note, the high correlation between macrophage abundance and the expression of CD163 expression indicates that this metric predominantly captures the relative abundance of M2-polarized macrophages (Supplementary Material). Tumor cells that progress in this scenario may thus, at least partly, keep the immunosurveillance in check through mechanisms that foster the recruitment of an immunosuppressive infiltrate. The upregulation of glycolysis, also overrepresented among these tumors, may favor this suppressive infiltrate via the depletion of glucose in the microenvironment. Indeed, recent studies have shown that the differentiation of immune cells into cytotoxic populations relies on the availability of glucose rather than other carbon sources (46, 47).
The third scenario describes the progression of tumors with poorly cytotoxic immunophenotypes (Fig. 4, bottom), which present high aneuploidy. These tumors also tend to overexpress, and/or present genomic alterations of, genes involved in cell processes that are key to cell growth and proliferation, such as the regulation of cell cycle, cell division, telomere maintenance, and mRNA and protein synthesis. Some of these molecular traits may contribute to the disruption of the cancer-immunity cycle, and thus drive immune evasion, whereas others may constitute opportunistic events fostered by the lowly cytotoxic immune infiltrate (48, 49). Examples of the former are the mutations of the IDH1 gene, in brain tumors, and of genes of the WNT–β-catenin pathway, as well as the upregulation of the TGFβ and Hedgehog signaling pathways, all of which have been reported to reduce immune cell infiltration in certain conditions (31, 32, 35, 50, 51). These tumors also exhibit worse prognosis across several cancer types, which is partially explained by their further advanced pathologic stage. The lower antitumor pressure exerted by these immunophenotypes also allows the growth of more tumor clones, resulting in more heterogeneous malignancies with greater metastatic potential (52).
Discussion
Here, we describe a comprehensive landscape of tumor immune infiltration pattern across 29 solid cancers and identify tumor features that are significantly associated with distinct immunophenotypes. Although the estimation of the immune infiltration pattern based on the expression of the bulk sample presents several limitations, the size of the cohorts included in the study provides the power to detect these associations. New methods that address some of these limitations will contribute to refine these analyses (21). Our study presents certain methodologic particularities that set it apart from previous works. First, to our knowledge, we have evaluated the most comprehensive, to date, catalog of genomic and transcriptomic characteristics of tumors that may be positively selected under the pressure applied by the immune system. Furthermore, in the case of genomic events, we have focused specifically on those most likely under positive selection, that is, likely drivers according to the Cancer Genome Interpreter (11). Second, we have adjusted the values of gene expression for the contribution of the immune content in the bulk sample, thus enhancing the detection of the upregulation of tumor-specific processes. Third, we have represented the immune infiltration as an array of immune populations, which revealed specific patterns associated to tumor characteristics that would have been missed otherwise. Finally, because we conducted a pan-cancer analysis, our results provide a landscape of the distribution of diverse mechanisms of escape from immunosurveillance across a comprehensive set of solid malignancies.
Certain genomic and transcriptomic features of tumors that we have found associated to particular immunophenotypes have been reported by previous studies with approaches similar to ours (6–8, 10). However, as a result of the methodologic particularities of our study (described above), we have also been able to identify novel features of the tumors that are significantly associated to particular scenarios of immune infiltration. Some of these features correspond to alterations of genes that constitute new candidates to drivers of immune resistance. Both confirmed and novel findings are summarized in the final section of Results. Furthermore, we have extended to the pan-cancer realm some drivers of immune evasion that had previously been identified in particular cancer types (35). Among these are the overexpression of TGFβ, as well as cell proliferation-related processes. To our knowledge, this is the first study that, by virtue of its completeness, is able to delineate three broad pan-cancer scenarios of immune infiltration and the mechanisms that allow tumors to evolve in each of them.
It is not clear from our results to what point the three scenarios of immune infiltration delineated in the final section of Results represent tumors with separate evolutionary pathways or different stages of their progression. In the latter case, tumors with highly cytotoxic immunophenotype would be at their early stages, and equipped with mechanisms to curb a robust immune pressure. They would then progressively evolve to invade neighboring tissues while their immunophenotype shifts toward a more immunosuppressive infiltration pattern. Tumors with a poor cytotoxic infiltrate would represent advanced stages of a malignancy in full progression and virtually unchecked by the host's immune system.
The most immediate result of our analysis is the identification of tumor features that are positively selected due to their interaction with the immune infiltrate. In this regard, we have described a general approach to classify tumors on the basis of their immunophenotype. More long-term outcomes of this work include understanding how these somatic genomic events and transcriptomic programs may actively shape the immune infiltrate. This knowledge shall have important implications for the design of novel therapeutic approaches (36). Currently, the most widely exploited immunotherapies involve immune checkpoint inhibitors. We explored transcriptomics data of small cohorts of patients with metastatic melanomas that received anti-CTLA4 and anti–PD-1 therapies (42 and 28 patients, respectively; refs. 14, 15). Patients with tumors with poorly cytotoxic immunophenotypes tend to respond worse to the anti-CTLA4 treatment, whereas patients with tumors that overexpress genes involved in the WNT–β-catenin, extracellular matrix organization, and angiogenesis pathways respond worse to anti–PD-1 (see Supplementary Note and Supplementary Fig. S15). Further analysis of larger cohorts of patients treated with immunotherapies will continue to clarify which characteristics of the tumors and/or their microenvironment are better predictors of their response. We also envision that the inclusion of drugs able to counteract mechanisms of immune-resistance or immunosuppression (such as the ones identified in our study) in clinical trials will reveal their potential for synergies with existing immunotherapies. Some of our findings, such as those supporting the use of drugs targeting the cell cycle to favor a transformation of the immune infiltrate, and TGFβ inhibitors or epigenetic modulators to improve the response to negative immune checkpoint inhibitors (33, 37, 48), constitute promising examples to explore in this direction. In summary, this study provides a comprehensive landscape of the characteristics of solid tumors that may influence (or be influenced by) the characteristics of their immune infiltrate. These results may help interpret the response of solid tumors to immunotherapies and guide the development of novel drug combination strategies.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: D. Tamborero, R. Dienstmann, A. Gonzalez-Perez
Development of methodology: D. Tamborero, C. Rubio-Perez, R. Dienstmann
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Tamborero, C. Rubio-Perez, F. Muiños, R. Sabarinathan, J.M. Piulats, A. Muntasell, R. Dienstmann
Writing, review, and/or revision of the manuscript: D. Tamborero, C. Rubio-Perez, F. Muiños, J.M. Piulats, A. Muntasell, R. Dienstmann, N. Lopez-Bigas, A. Gonzalez-Perez
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Tamborero
Study supervision: D. Tamborero, R. Dienstmann, N. Lopez-Bigas, A. Gonzalez-Perez
Acknowledgments
D. Tamborero is supported by project SAF2015-74072-JIN, which is funded by the Agencia Estatal de Investigacion (AEI) and Fondo Europeo de Desarrollo Regional (FEDER). C. Rubio-Perez is supported by an FPI fellowship (BES-2013-063354) from the Spanish Ministry of Economy and Competitiveness. N. Lopez-Bigas acknowledges funding from the European Research Council (consolidator grant 682398). A. Gonzalez-Perez is supported by a Ramón y Cajal contract (RYC-2013-14554), which also covers the costs of this publication. The results shown here are in whole or part based upon data generated by the TCGA Research Network. Data generated by the Genotype-Tissue Expression (GTEx) Project, supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS, were also used. We thank Elena Gros (PhD) and Ignacio Melero (MD, PhD) for the scientific review and helpful discussions of the study, Fiorella Ruiz for the annotation of the TCGA clinical data, and Dvir Aran (PhD) for the support in the expression adjustment methodology. We also thank Eliezer Van Allen (MD, PhD) for providing the transcriptomic data of melanoma patients treated with anti-CTLA4 therapy.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
References
Supplementary data
Supplementary Methods, Supplementary Note, and Supplementary Figures
Immune-phenotypes
Results of the GSEA enrichment