A Pan-cancer Landscape of Interactions between Solid Tumors and Infiltrating Immune Cell Populations

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.

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.

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).

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.

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).

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.

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.

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.

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).

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.

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.

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.

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.

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.

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.

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.

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).

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.

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).

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).

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.

No potential conflicts of interest were disclosed.

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

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.