Clinical relevance of tumor immune infiltration

To study the distribution of infiltrating immune cells in the tumor and adjacent/normal tissues, we focused on 18 cancer types for which the mRNA expression profiles of adjacent or normal tissues were available. Consistent with Rooney et al. [6], CD8 T cells are enriched in tumor tissues in kidney cancer (KIRC) and head and neck cancer. In contrast, CD8 T cells appear to be in lower abundance in most other cancers, such as non-small cell lung carcinomas (including adenocarcinoma and squamous cell carcinoma) and colorectal cancer (including colon and rectal adenocarcinoma) (Fig. 2a). Macrophages are significantly enriched in GBM, which is supported by previous observations showing that microglia and macrophages are present in large numbers in the glioma microenvironment [21]. Furthermore, the abundance of tumor-infiltrating B cells is significantly higher than in the adjacent or normal tissues in eight cancer types. Interestingly, B-cell infiltration predicts a significantly better outcome in a subset of these cancers (Fig. 2b). The most dramatic case is GBM, where patients with the top 20 % of B-cell infiltration have a 4.7-month longer median survival time than those with the lowest 20 % (“Methods”). This result suggests that tumor-infiltrating B cells play an important role in the antitumor responses in GBM and lung and bladder cancers.

Open in a separate window

Fig. 2

Distribution of infiltrating immune cells and selective enrichment of B cells in the tumor microenvironment. a The abundance of infiltrating CD8 T cells, macrophages, and B cells in 18 cancer types, with both primary tumor and adjacent (Adj)/normal (Norm) tissue available. Normal tissue was from healthy donors where adjacent tissues were unavailable. Statistical significance was evaluated by Wilcoxon rank sum test. Blue arrowheads point to three cancers with B cells significantly enriched in the primary tumor and associated with clinical outcomes. q values are colored red, blue, or black for significant (false discovery rate≤0.15) enrichment in tumor, adjacent or normal tissue, or insignificance, respectively. b B-cell infiltration level significantly predicted patient survival in selected cancer types. Tumors in the top 20th percentile of B-cell infiltration were compared with those in the bottom 20th percentile. The median survival time for the top 20 % of patients with brain, lung, and bladder cancers was 460, 1778, and 2000 days, respectively, and for the bottom 20 % 345, 976, and 575 days, respectively. Statistical significance and hazard ratios (HR) with 95 % confidence intervals were calculated for all the samples, not just the top and bottom 20 %, using multivariate Cox regression including all six immune cell types, patient age, and clinical stage

We next investigated how infiltrating immune cells influence clinical outcomes. Using multivariate Cox regression adjusted for age, stage, and viral infection status, we identified many significant associations between immune cell abundance and patient outcomes (Fig. 3). While CD8 T-cell infiltration associates with prolonged survival and macrophage infiltration consistently predicts worse outcome, other immune cell types have cancer-specific effects on prognosis (Fig. 3a). These observations corroborate previous reports that cytotoxic T-cell infiltration independently predicts better outcome in liver [22] and rectal [23] cancers. Consistent with established knowledge in melanoma and head and neck cancer [24, 25], we also found infiltrating CD8 T cells to be associated with survival in univariate analysis (Fig. 3b). However, it is not an independent predictor of better outcome after adjusting for other covariates, since CD8 T cells are correlated with neutrophil infiltration in melanoma and human papilloma virus (HPV) infection in head and neck cancer (“Methods”). In fact, infiltrations of B cells and CD4 and CD8 T cells are all significantly higher in HPV-positive than in HPV-negative head and neck tumors (Additional file 1: Figure S4), suggesting that viral antigens result in an elevated lymphocyte response. Besides the effect on patient survival, CD8 T cells may also play an important role in preventing tumor recurrence (Fig. 3c). In melanoma and colorectal and cervical cancers, patients with higher CD8 T-cell infiltration in the primary tumors have a significantly lower risk of developing a second tumor during the follow-up period. Overall, our observations on tumor-infiltrating CD8 T cells are extensively supported by clinical studies [14] (Additional file 4: Table S3), thus providing additional validation to our deconvolution method. Cytolytic activity (CYT) is a previously defined metric of immune-mediated cell destruction [6]. Compared with CYT, our analysis identified many more strong associations between tumor-infiltrating immune cells and patient clinical outcomes (Fig. 3a), presumably because our method takes into account expression data from hundreds of genes instead of only two.

Open in a separate window

Fig. 3

Immune cell infiltration predicts clinical outcome. a Association of tumor infiltrating immune cells with patient survival. For each cancer type, multivariate Cox regression was performed, with covariates including the abundance of six immune cell types, patient age at diagnosis, clinical stage, and viral infection status when available. Each entry on the first six rows of the heatmap represents the hazard ratio (HR) of a corresponding immune cell type, with larger size indicating statistical significance at a false discovery rate (FDR) of 0.15 and color indicating the value of the HR. The last row of the heatmap records the Cox model HRs and statistical significance using cytolytic activity (CYT) scores adjusted for the same covariates. Multiple test correction was performed using q value across cancer types and six immune components. b Kaplan–Meier curves of melanoma (SKCM) and head and neck cancer (HNSC) stratified by infiltration CD8 T-cell abundance. Median survival time for the top 20 % of patients in melanoma and head and neck cancers is 4507 and 1838 days, respectively, and for bottom 20 % 2005 and 862 days, respectively. Statistical significance, hazard ratios, and 95 % confidence intervals were calculated using multivariate Cox regression and all the samples as described above. c CD8 T-cell infiltration in primary tumors (metastatic samples for SKCM) significantly (FDR≤0.15) predicts tumor relapse in selected cancers. Statistical significance was evaluated using logistic regression correcting for patient age and clinical stage

TAM numbers have been reported to be a predictor of worse outcome in many cancers [13]. Consistently, we found that TAM significantly associates with worse outcome in bladder, breast, and ovarian cancers and in lower-grade glioma (Fig. 3a), supporting TAMs as an independent prognostic factor for these cancers. Extending this analysis to the less well-studied chromophobe renal carcinoma (KICH), we detected a significant inverse association between macrophage infiltration and patient survival (Fig. 3a), suggesting that TAMs function in KICH as in other solid tumors.

Potential causes for immune infiltration heterogeneity

We next examined the possible causes of inter-tumor immune infiltration heterogeneity, such as somatic mutations and chemokine expression. Non-synonymous somatic mutations in the tumor genome can generate immunogenic neoantigens that trigger antitumor response through T-cell activation [26–29]. In addition, total mutation load has been suggested as a surrogate for neoantigen count [6]. To understand how the host immune system responds to tumor somatic mutations, we studied the association of infiltrating immune cells and total mutation load. Tumor purity is a confounding factor in this analysis, since purity affects the power to detect somatic mutations [20] and drives the pattern of gene expression [16]. After correcting for purity, we observed positive correlations between the total mutation count and infiltrating immune cells in a subset of cancers (Fig. 4a). In lung squamous cell carcinoma, lower grade glioma, HPV-negative head and neck cancer, and prostate cancer, CD8 T-cell infiltration increases significantly with tumor mutation load, in support of a previous study [6]. In addition, we observed that dendritic cell infiltration is correlated with the total mutation load in breast cancer (Spearman’s ρ=0.11, q=0.037), as is B-cell infiltration (ρ=0.13, q=0.018), suggesting cancer-specific roles for these cell types in antitumor immunity. Similar associations were observed when we used neoantigen load estimated from a previous work [6] instead of total mutation load (Additional file 1: Figure S5), corroborating our conclusion that increased tumor neoepitope load elevates the infiltration of multiple immune cell types.

Open in a separate window

Fig. 4

Potential causes of inter-tumor immune infiltration heterogeneity. a In selected cancer types, counts of total somatic coding mutations positively associated with the level of infiltrating immune cells. The y-axis is the residual of corresponding immune cell abundance after linear regression against tumor purity. Statistical significance was evaluated using partial Spearman’s correlation adjusted for tumor purity. The asterisk indicates only HPV-negative tumors were selected for head and neck cancer. Multiple test correction was performed across cancer types and six immune components. Gray hues indicate previously known results (HNSC, LGG, and LUSC), while other findings are novel in this study. b CD8 T-cell infiltration is associated with microsatellite instability (MSI) status in cancers commonly with MSI. MSI-H high level of microsatellite instability, MSI-L low level of microsatellite instability; MSS microsatellite stable. Statistical significance was evaluated using a Wilcoxon rank sum test. c, d Chemokine/receptor networks for immune infiltration in diverse cancer types. Vertexes are ligands (green) and receptors (purple) and edges indicate known molecular interactions. For each cancer, the partial correlations (corrected for purity) between the chemokine gene expression and infiltration of CD8 T cells (c) or macrophages (d) were calculated. For a pair of interacting chemokine and receptor genes, if both are significantly correlated with immune cell infiltration in one cancer, a colored dot represents the cancer type is placed on the edge connecting the chemokine and receptor. Statistical significance was evaluated using partial Spearman’s correlation at a false discovery rate threshold of 0.01. Heatmap visualizations of the same results are shown in Additional file 1: Figure S6c, d

Besides point mutations, microsatellite instability (MSI) is seen in colorectal, stomach, and endometrial cancers. MSI typically generates small indels across the genome, producing non-self antigens that may be recognized by the host immune system. Consistent with a previous report [7], we found CD8 T cells to be significantly more abundant in MSI-high (MSI-H) tumors compared with MSI-low (MSI-L) tumors in colon cancer (Fig. 4b). Among the remaining three TCGA cancers with available MSI information, we also found higher levels of MSI to be associated with increased CD8 T cells in stomach cancer. A recent study reported that MSI-high colon cancer patients showed significantly better responses to PD-1 blockade therapies [12] and our results suggest that this conclusion may be extended to other gastro-intestinal cancers with MSI.

To further investigate the regulation of immune infiltrates in different cancers, we also systematically studied the expression levels of chemokines and receptors. Most of these molecules were expressed in the microenvironment (Additional file 1: Figure S6a, b). CD8 T-cell level is significantly associated with a subset of chemokine–receptor pairs, including CCL3,4,5–CCR1,5 and XCL1,2–XCR1 (Fig. 4c). On the other hand, different molecules are associated with macrophage abundance in cancer-specific patterns. Macrophage infiltration appears to be related to CXCL12–CXCR4 in thyroid, head and neck, stomach, and colon cancers and to CCL14,CCL23–CCR1 in lung cancers (Fig. 4d). Our results highlight potential bases for inter-tumor heterogeneity in immune cell infiltration and suggest possible means for reducing macrophage recruitment to the tumor microenvironment.

Inclusion criteria for immune cell types

In order to minimize co-linearity in the regression analysis and maximize the robustness of our inference, our study focused on six immune components based on two criteria. First, the reference data contain at least ten independent samples of the immune cell type. Second, if the expression profile of a given cell type is highly correlated (sample-wise Pearson’s r≥0.9) with other cell type(s), we chose the cell type with more samples. The selected cell types represent the finest resolution of immune cell lineages that we can achieve based on the above inclusion standards. Cell types excluded from the inference may affect the highly correlated immune components included. Improved reference immune datasets will be needed to deconvolve individual cell types.

Computational method for immune cell composition deconvolution

We first estimated sample purity for each tumor through DNA SNP array data using our previously developed tool CHAT [15]. Genomic estimations of tumor purity have been validated using diluted series of cancer and blood cell lines with known mixture ratios [16]. To clarify, the quantity inferred from CHAT is the fraction of aneuploid cells. In this work, we used this quantity as a surrogate of tumor purity. Samples with percent on point (PoP) <0.01 were excluded, where PoP is a quality measure reported by CHAT. In each cancer, we selected genes with expression values negatively correlated with tumor purity (Pearson’s r≤−0.2, P≤0.05), denoted as G_p, and intersected with G_i. Our goal in this step is to select informative genes with expression levels strongly affected by tumor purity. It is not important whether the gene is highly or lowly expressed. Pearson’s correlation is suitable because it is a measure of the linear dependence of two random variables, disregarding the magnitude of the observations. The resulting smaller gene set intersection was denoted as G₀, which is cancer-type specific. Meanwhile, we merged the tumor gene expression with the reference immune cell data of all genes using ComBat [17]. According to the principal components analysis plot, ComBat effectively removed the batch effect between different datasets generated using different platforms (Fig. 1b). For each of the six cell types, we then calculated the median expression value in all the samples available for that cell type for each gene in G₀. The resulting dataset (referred to as R) contains six vectors of gene expression values, each for one cell type. We used R to filter G₀ in the following way: for each vector in R, genes with expression values in the top 1 % in G₀ are removed. This is because genes expressed at extremely high levels in the reference dataset will dominate the inference results. Since highly expressed genes have large variance, our inferences become very sensitive to these outliers. We therefore remove them to acquire more robust estimations. The total removed genes were the union of the top 1 % from of the six vectors. The resulting gene set was denoted as G_0f. For each sample, a constrained least square fitting described in a previous study [19] was applied to infer the relative abundance for each of the six immune components (Fig. 1e). The predictions of this method were validated using mixtures of different blood cell types with known ratios. For a given sample, let Y^g denote the gene expression of gene g, where g G_0f. Let X^g_r denote the gene expression of gene g in immune cell type r (r=1,2,…,6) in reference dataset R. Since we model that the given sample is a mixture of the six immune cell types, our task is to find positive coefficients f to minimize the total squared differences:

which is a constrained linear regression problem. The estimation accuracy of f should be affected by the genes used in the fitting, although our model assumes that such an f exists and should be the same for all genes. To note, coefficients f estimated using this method are the relative abundance of immune cells. The scale of the estimation of an individual immune cell type is determined by the variance of the corresponding reference data X_r. Therefore, f are not comparable between cancer types or different immune cells. Source codes for TIMER and downstream statistical analysis as well as related data files are available at http://cistrome.org/TIMER/download.html.

Correction for reference immune cell colinearity

Although the six immune cell types in this study are selected in such ways that the colinearity between cell types is minimized, we found that in THCA and UCS, CD4 and CD8 T-cell signatures are still very similar. Consequently, the inferred CD4 and CD8 T-cell levels are negatively correlated (Pearson’s r≤−0.3), which is an artifact of covariates’ colinearity in the constraint regression. Additional analysis on these two cancers revealed that the negative correlation is driven by a small number of CD4 or CD8 T-cell signature genes that are extremely overexpressed in the tumor samples. We remove the union of the top expressed gene in each tumor sample and re-estimate f. This step is repeated until the correlation between estimated CD4 and CD8 T-cell levels is larger than −0.3. This analysis provides more robust estimations of immune cell abundance in cancer types.

Pathological estimation of neutrophil infiltration in BLCA

For the TCGA data sets, the original samples are unavailable for further studies; however, hematoxylin and eosin (H&E) digital slides have been publicly released. While it is not possible to distinguish T cells and B cells by H&E, neutrophils are morphologically distinctive and their abundance can be estimated. Occasional dendritic cells and macrophages can be identified by H&E but their true abundance is difficult to estimate in the absence of immunohistochemical stains. Slides were reviewed via digital images viewed with the Cancer Digital Slide Archive (http://cancer.digitalslidearchive.net). We chose BLCA because it has a large sample size (n=404), does not have excessive necrosis, and has sufficient neutrophil counts and sample variety to allow for validation by histological evaluation. The pathologist reviewing the slides was blinded to the in silico neutrophil predictions. Samples were stratified into three groups (high, medium, low) with levels of neutrophils relative to the entire collection of samples.

DNA methylation-based total leukocyte estimation

The percentage of total DNA for each tumor was estimated using DNA methylation data as previously described and compared with purity estimates derived from SNP data [20, 37]. In brief, the estimates were derived from loci with tissue-specific methylation that distinguishes the corresponding tissue type from lymphocytes using the TCGA Pancan12 [38] DNA methylation dataset. The leukocyte methylation signature was derived as follows. Each probe was ranked by the difference in mean beta value in buffy coat and corresponding normal samples. We retained the 100 probes with the largest positive difference and the 100 with the largest negative difference between mean DNA methylation in normal tissues and peripheral blood leukocytes, designated NT and BC probes (hypermethylated in normal tissue compared with buffy coat, and vice versa, respectively). Let T_ik denote the beta value for probe k in tumor sample i. Let B_k denote the average beta value of buffy coat samples for each probe. Let T_k denote the minimum observed beta value across all tumor samples for the BC probes and the maximum for the NT probes, which theoretically reflects the ground state of methylation level in the purest tumor. Denote with f_B the fraction of buffy coat (leukocyte) components in the sample, then, assuming a linear relationship, we have the following equation for each probe:

Solving this equation for f_B gives:

The values of f_B for each of the 200 probes in the signature were calculated and a kernel density estimate of f_B was obtained. The leukocyte signature was then calculated as the mode of this density estimate.

Monte Carlo sampling and in silico validation

We validated our predictions on infiltrating immune cell abundance using in silico simulated data. As mentioned, for each cancer we selected a gene set G_0f (length n₀) for least squares fitting. In order to control for the mixing ratios of the six components while maintaining the correlation structure of the real data, we first calculate the gene–gene covariance matrix Σ_0f for all the genes in G_0f using tumor expression data. We then randomly sample six numbers f_1-6, from Uniform(0,1). We calculate μ_0f (length n₀), which is the average of six immune components weighted by f_1-6. Next, we sample a vector of length n₀ from multivariate normal distribution with mean μ_0f and covariance Σ_0f. For each cancer type, we simulated the same number of samples as its sample size in the TCGA data. After applying our method, we compared the estimated immune abundance with true values using Pearson’s correlation. Low quality estimations with Pearson’s r≤0.2 were excluded from the downstream analysis.

Selection of cancer specific genes

For each cancer type, we compare tumor samples with all normal samples collectively. Only genes overexpressed in tumor samples and absent or expressed at lower levels in all normal tissues were selected. Differentially expressed genes were selected based on a FDR ≤0.05 and at least a twofold difference in expression levels. In the case of tumors with established clinical subtypes, such as breast cancer, we selected the top 25 samples for each gene based on their rank of raw read counts, then identified differentially expressed genes within each subtype. The final tumor-specific gene set was the union of all the cancer types (or subtypes).

Statistical analysis

Multivariate Cox regression, log-rank test and Kaplan–Meier estimators were implemented using the R package survival. The association between CD8 T-cell abundance and tumor status was evaluated using logistic regression corrected for age and clinical stage and was implemented using the R package glm. The same analysis was performed for neutrophil abundance and gender associations, corrected for age and smoking history. Partial correlations of immune cell abundance and gene expression of chemokines and receptors, somatic mutation counts, CT gene expression, as well as immunosuppressive molecule expression were calculated using the R package ppcor. Multiple test correction was performed using the R package qvalue [39] and FDR thresholds are applied based on the abundance of signals in the data. In this study, we applied the Pearson correlation to purity and gene expression because it is reasonable to expect that the expression level is linearly associated with tumor purity. For others, we used the Spearman correlation. We applied partial correlation analysis to remove the influence of tumor purity on the involved variables. All other analyses, including linear regression, Fisher’s exact test, Wilcoxon rank sum test, Spearman’s correlation, and hierarchical clustering, were performed using R [40]. Of note, in Figs. 2b and ​and3b,3b, we used the 20 percentile as a cutoff only to help visualize the association of immune infiltration with outcomes and the statistical significance was determined by multivariate Cox regression (Fig. 3a) including all the samples. Our results on survival analysis, neoantigen association, tumor recurrence, and association of checkpoint blockade inhibitory molecules with immune cells are available in Additional file 10: Table S8.

Additional analysis on HNSC and SKCM

One intriguing result we observed is that univariate and multivariate survival analysis results for HNSC and SKCM are not consistent (Fig. 3a, b). For HNSC, we discovered that HPV infection, a recently identified prognostic factor [41], correlates with CD8 T-cell infiltration (Additional file 1: Figure S4). It is likely that the previously observed association of CD8 T cells with survival is contributed to by virus infection. On the other hand, for SKCM, we found that the infiltration level of CD8 T cells is highly correlated with neutrophils (Pearson’s r=0.79) and dendritic cells (r=0.81), indicating that these immune cells work in concert. As highly correlated features confound each other in a multiple regression, we performed principal component analysis on the abundance of the six immune cells. We reanalyzed the Cox model using six principal components (PCs), age and stage as covariates, and found PC1 (hazard ratio (HR)=2.6×10⁻⁴, p=0.0062), PC4 (HR=1.07×10², p=0.033) and PC6 (HR=0.01, p=0.008) to be significantly associated with survival. PC1 was comprised of CD8 T cells, neutrophils, and dendritic cells (by Pearson’s correlation), thus capturing the colinearity in the data. PC4 represented macrophages and predicted worse outcome. There was no clear assignment of PC6 to any immune component(s) and it may represent an unselected immune cell type.

Additional analysis on OV and BRCA

In the survival analysis (Fig. 3a, b), we failed to identify some known prognostic predictors, notably B cells in OV and CD8 T cells in BRCA. A previous study reported that CD20 cells positively associate with survival [42]. We investigated the expression levels of the B-cell markers CD19 and CD20 in OV and discovered that tumor purity is not negatively correlated with gene expression levels for both genes, indicating that aneuploid cells in ovarian cancer may also express B-cell markers. Therefore, cell sorting based on CD19 or CD20, which is not the B-cell component in our analysis, is likely to select cancer cells. CD8 T cells were previously reported to associate with better outcome in BRCA [43], although we did not observe this relationship. This is possibly due to insufficient follow-up time or fewer deaths in the TCGA BRCA data, which thus underpowered our survival analysis.

Patient samples for IHC

De-identified clear cell renal cell carcinoma (ccRCC) formalin-fixed and paraffin-embedded tissue blocks from cases included in the TCGA KIRC cohort were obtained from the department of Pathology at the Brigham and Women’s Hospital. Patients had provided an informed consent for use of specimens and baseline and prospective clinical data for research purposes. The study was approved by the Dana-Farber/Harvard Cancer Center (DF/HCC) institutional review board. In total, five TCGA samples were selected for IHC: TCGA-CZ-5453-01A (negative control), TCGA-CZ-4866-01A, TCGA-CZ-4863-01A, TCGA-CZ-5459-01A, and TCGA-CZ-4862-01A.

We thank Nir Hacohen, Gordon Freeman, Glenn Dranoff, Toni Choueiri, Eliezer Van Allen, and Catherine Wu for their helpful discussions during manuscript preparation and revision.