Integrative analysis of imaging and transcriptomic data of the immune landscape associated with tumor metabolism in lung adenocarcinoma: Clinical and prognostic implications.

Although metabolic modulation in the tumor microenvironment (TME) is one of the key mechanisms of cancer immune escape, there is a lack of understanding of the comprehensive immune landscape of the TME and its association with tumor metabolism based on clinical evidence. We aimed to investigate the relationship between the immune landscape in the TME and tumor glucose metabolism in lung adenocarcinoma.
Methods: Using RNA sequencing and image data, we developed a transcriptome-based tumor metabolism estimation model. The comprehensive TME cell types enrichment scores and overall immune cell enrichment (ImmuneScore) were estimated. Subjects were clustered by cellular heterogeneity in the TME and the clusters were characterized by tumor glucose metabolism and immune cell composition. Moreover, the prognostic value of ImmuneScore, tumor metabolism and the cell type-based clusters was also evaluated.
Results: Four clusters were identified based on the cellular heterogeneity in the TME. They showed distinct immune cell composition, different tumor metabolism, and close relationship with overall survival. A cluster with high regulatory T cells showed relative hypermetabolism and poor prognosis. Another cluster with high mast cells and CD4+ central memory T cells showed relative hypometabolism and favorable prognosis. A cluster with high ImmuneScore showed favorable prognosis. Multivariate Cox analysis demonstrated that ImmuneScore was a predictive prognostic factor independent of other clinical features. Conclusions: Our results showed the association between predicted tumor metabolism and immune cell composition in the TME. Our studies also suggest that tumor glucose metabolism and immune cell infiltration in the TME can be clinically applicable for developing a prognostic stratification model.

Introduction

Cancer cells generate an immunosuppressive tumor microenvironment (TME), which suppresses T cell infiltration and activation through recruitment of the immunosuppressive cellular populations as well as T cell anergy, exhaustion, and senescence 1, 2. Accumulation of regulatory T cells (Treg) and myeloid-derived suppressor cells in the TME is well-demonstrated in immunosuppressive populations, which eventually impacts patients' prognosis 3, 4. Tumor-associated macrophages also interfere with induction of anti-tumor responses by changing their subpopulation 5. Recent studies suggest that one of the key factors leading to these immunosuppressive TMEs is abnormal metabolic milieu caused by cancer cells. Aerobic glycolysis of cancer cells restricts glucose consumption by T cells, which dampens their functions as glucose is essential for the differentiation and effector function of T cells 6-8. Lactic acid derived by aerobic glycolysis of cancer cells changes the normal immune reactions by polarizing tumor-associated macrophages 9. In spite of the recent cell-level evidence that cancer cell metabolism affects immune cells, how the metabolic milieu systematically impacts heterogeneous immune cells of the TME has not been fully understood. So far, no comprehensive analysis dealing with the association of the immune landscape of the TME, tumor metabolism and the clinical course has been conducted based on clinical data.

To this end, we conducted studies with clinical data of lung adenocarcinoma (LUAD) from The Cancer Genome Atlas (TCGA) projects 10, hypothesizing that the immune landscape of the TME is closely associated with tumor metabolism and also affects patients' prognosis. Using transcriptome data, we estimated the immune cell composition of the TME and built a model for predicting tumor metabolism. The model was based on gene networks related to glucose uptake measured by 18F-Fluorodeoxyglucose (FDG) positron emission tomography (PET) data 11-13. We clustered subjects according to the cellular components of their TME and characterized each cluster in terms of tumor metabolism, prognosis, and the immune cell composition. Moreover, the prognostic predictive value of overall immune cell enrichment score and tumor metabolism was investigated. Here, we present the association of the immune landscape of the TME and tumor metabolism, which can be a clue for cancer immunoediting mechanisms as well as a clinically feasible prognostic biomarker.

Materials and Methods

Study design

This study was an observational study using transcriptome and FDG PET image data. We developed a model for estimating the tumor metabolism index (TMI) of all LUAD samples (n = 490; Table ) of TCGA projects, and the comprehensive immune landscape was evaluated by calculating immune cell enrichment scores. The subjects were divided into four distinct clusters based on cell composition of the TME, and the characteristics of each cluster, including the differences in TMI and each immune cell type enrichment, were investigated. The prognostic value of the immune landscape, TME-based clusters and tumor metabolism was also evaluated.

Coexpression gene networks correlated with tumor metabolism

The preprocessing of both transcriptome and FDG PET data is described in Supplementary Material. To generate the model for estimating tumor metabolism, a microarray dataset with a matched image dataset (GSE28827, n = 20) 13, 14 was combined with a large microarray dataset composed of 226 LUAD samples (GSE31210) 15, 16 as the training set. We used Weighted Gene Coexpression Network Analysis package (WGCNA) 17 to identify gene network modules from the training set. We created a correlation matrix on the basis of the Pearson's correlation coefficient for all pair-wise genes across all samples. The power—the key parameter for the weighted network—was selected to optimize both scale-free topology and sufficient node connectivity and we chose a threshold of 6 in this study. The correlation matrix was transformed into an adjacency matrix (matrix of connection strength) using the power function, and pair-wise topological overlap between genes was calculated. We identified network modules using a hierarchical clustering method with topological overlap dissimilarity as the distance measure 17. The modules were detected by the dynamic tree cut algorithm in WGCNA package, defining a height cutoff value of 0.99, deep split as 2, and minimum module size cutoff value of 50. Genes that were not assigned to any module were classified to color gray.

The expression matrix of the twenty samples with FDG PET data from the training set was extracted to conduct a correlation test. To summarize the module expression profile, we calculated the module eigengene, the first principal component of the expression matrix of the corresponding module. The Pearson's correlation analysis was performed between the module eigengene and the maximum standardized uptake value (SUV) measured by FDG PET. The modules with p-value corrected by false discovery rate < 0.05 were regarded as significant tumor metabolism-associated modules.

The tumor metabolism-associated modules were independently validated for TCGA data by a subset of LUAD samples including corresponding FDG PET data (n = 17 of 490). The preprocessed mRNA transcriptome data of TCGA was transformed into log2 scale after adding value one. The gene network of TCGA data was constructed using the same 4010 genes of the training set by imputing missing genes with k-nearest neighbors algorithm (k = 10) 18. The Pearson's correlation between the module eigengene of tumor metabolism-associated modules and maximum SUV was tested to see whether the modules identified from the training set were also significantly correlated with maximum SUV in TCGA data.

Tumor metabolism estimation model

The model based on neural networks was trained to predict maximum SUV from transcriptome data. The inputs of the model were gene expression values of the tumor metabolism-associated modules, which included 480 and 146 genes. Training data consisted of 20 subjects who have FDG PET image data and 226 subjects without image data. For the semi-supervised training, we used a simplified ladder network model for this task 19. To use unlabeled data, denoising autoencoder was used 20. It could be used for dimension reduction as well as feature extraction from high-dimensional data. The encoder of the denoising autoencoder extracts 64 features from gene expression data, and then another neural network layer that connects the features to one-dimensional output, maximum SUV, was added. This network was used for the supervised learning. The decoding of autoencoder was performed as an auxiliary task. Sixty-four features extracted by the encoder were decoded by a neural network layer connected to the original data dimension and the denoising autoencoder was trained to reconstruct the original expression data. These supervised learning and auxiliary unsupervised training, denoising autoencoder, were combined for predicting maximum SUV (Figure ). The output of the estimation model was defined as TMI. The encoder and decoder simply consisted of one layer, respectively. For optimizing the neural networks, Adam optimizer was used 21. Neural networks were constructed by a deep learning library, Keras (version 1.2.2, https://keras.io).

After the training, the model was tested for 17 TCGA data that have corresponding FDG PET images (Figure . The output of the model was defined as TMI. TMI of all LUAD samples was obtained for further analyses (Figure .

Cell types enrichment analysis

To evaluate the heterogeneous cellular landscape of TME, cell types enrichment scores were evaluated. A gene signature-based method for inferring cell types from tissue transcriptome profiles, the xCell tool (http://xcell.ucsf.edu/), was used 22. It infers 64 immune and stroma cell types of the TME. The composite score of immune cells, ImmuneScore, was also obtained.

Two-dimensional TME landscape mapping and clustering

Sixty-four cell types enrichment scores were assigned to each tumor sample via the xCell analysis. To visualize cellular heterogeneity in TME, a dimension reduction method, t-distributed stochastic neighbor embedding (t-SNE), was employed 23. Cell types enrichment scores were normalized by mean and standard deviation across the samples and they were projected onto 2D axes, defined as the 2D TME landscape map.

Using the values of two axes of the TME landscape map, a dendrogram representing a pairwise similarity was constructed. The similarity was determined by the average distance between elements of each hierarchical cluster. Cutting the dendrogram at a given height provides a number of clusters.

Statistical analysis

The comparison of TMI between groups (LUAD subtypes or cell type-based clusters) was conducted by one-way ANOVA, followed by the post hoc Tukey's test. The survival rate of the groups was depicted with the Kaplan-Meier's method and compared with the log-rank test. The prognostic properties of TMI, cell types-based clusters and ImmuneScore were evaluated by univariate and multivariate Cox regression analysis. To define the risk groups, continuous variables including TMI were dichotomized using the median value of each variable in TCGA data.

Results

Transcriptome-based tumor metabolism estimation

The tumor metabolism-associated modules were selected and then a neural network model for estimating FDG uptake was trained. Two modules among ten gene network modules were significantly correlated with tumor metabolism (magenta module: r = 0.60, p < 0.05; brown module: r = -0.567, p < 0.05; p-values corrected for false-discovery rate). These modules were also significantly correlated with tumor metabolism for TCGA samples including FDG PET (r = 0.54, p = 0.02 and r = -0.58, p = 0.01 for magenta and brown modules, respectively; Figure ). Both modules were associated with extracellular matrix organization processes (Table ). A tumor metabolism estimation model was built by genes of the two modules. The mean absolute error of the predicted maximum SUV was 1.52 for the training set and 2.79 for TCGA data. TMIs were estimated for all samples of TCGA data (n = 490) using the model.

TMIs were compared between molecular subtypes of LUAD (subtyping is described in Supplementary Material). The TMI of the bronchioid subtype was significantly lower and that of the squamoid subtype was significantly higher than the other subtypes (F = 211.7, p < 1×10-15; 3.87±1.09 vs. 4.60±1.11 vs. 6.45±0.96 for TMI of bronchioid, magnoid and squamoid subtypes, respectively; Figure ). The patients were divided into two groups according to TMI, and the group with high TMI showed significantly worse survival (p = 0.0019, Figure ). We compared TMI with two other metabolic signatures representing glucose metabolic activity (detailed methods are described in Supplementary Material). The glycolysis and gluconeogenesis pathway score based on mean expression of selected genes 24 was positively correlated with TMI (r = 0.26, p < 1×10-8; Figure ). We also obtained enrichment scores of the glycolysis pathway using Reactome 25. A significant positive correlation between glycolysis enrichment score and TMI was also found (r = 0.32, p < 1×10-12; Figure ).

Clusters based on tumor microenvironment cell types and tumor metabolism

Sixty-four different cell types enrichment scores and ImmuneScore were computed for TCGA data. Using these scores, a 2D TME landscape map based on t-SNE was drawn with TMI (Figure ). A trend that relatively hypometabolic tumors were predominantly located towards low values on both axes was identified (Figure ). Subjects were divided into four clusters (C1-C4) according to their similarity in TME cell composition (Figures ). There was a significant effect of the cell types-based clusters on TMI (F = 29.4, p < 1×10-15; Figure ). Post-hoc analyses revealed that the TMI of C2 was significantly lower than those of the other clusters. In contrast, the TMI of C3 was significantly higher than those of the other clusters (Table ).

Immune landscape of cell type-based clusters

Immune cell-type enrichment scores were compared between clusters. Tumors with high ImmuneScore tended to locate on the left portion of the 2D TME landscape map (Figure ), and ImmuneScore was significantly different between clusters (F = 43.84, p < 1×10-15; Figure ). Post-hoc analysis demonstrated that the C4 showed significantly higher ImmuneScore than the other clusters (Table ). The 2D TME landscape map with each cell type enrichment score (Figure ) and distribution of scores of clusters (Figure ) are presented. Additionally, a heatmap of immune cell type enrichment scores is illustrated to identify the comprehensive immune landscape (Figure ). These illustrations reveal that immune cell enrichment scores of most immune cells including CD4+, CD8+ T cells, macrophages, and B cells were relatively higher in C4 than the other clusters. A cluster with high TMI, C3, showed a significantly higher score of Tregs. Another cluster with low TMI, C2, showed relatively higher scores of mast cells and CD4+ central memory T cells. Most of the immune cell enrichment scores were low in C1; however, basophils and Th1, Th2 cells were relatively enriched in C1.

Survival analysis

Overall survival was associated with the cell type-based clusters (p = 0.029; Figure ). C2 was associated with good prognosis compared with C1 (HR = 0.58, 95% C.I. 0.38-0.89; p = 0.01). C4 also showed a trend of favorable prognosis (HR = 0.60, 95% C.I. 0.38-1.01; p = 0.05) (Table ). The high ImmuneScore group showed significantly better prognosis (HR = 0.65, 95% C.I. 0.48-0.88; p = 0.005; Figure ). LUAD subtypes were also associated with overall survival (p < 1×10-4; Figure ). Magnoid and squamoid subtypes showed significantly poorer prognosis compared to bronchioid subtypes (HR = 1.99, 95% C.I. 1.35-2.92 and HR = 1.84, 95% C.I. 1.29-2.62; p = 0.0005 and 0.0008 for magnoid and squamoid subtypes, respectively). TMI and tumor stage also showed significant association with overall survival (HR = 2.62, 95% C.I. 1.91-3.60; p < 0.0001 for stage III or IV tumors; HR = 1.61, 95% C.I. 1.19-2.19; p = 0.002 for high TMI group). Multivariate Cox analysis demonstrated that tumor stage and ImmuneScore were significant independent predictors for overall survival (HR = 2.31, 95% C.I. 1.67-3.20, p < 0.0001 for tumor stage; HR = 0.69, 95% C.I. 0.48-0.98, p = 0.04 for ImmuneScore; Table ). In addition, C2 was also associated with favorable prognosis with marginal significance (HR = 0.63, 95% C.I. 0.40-1.01; p = 0.05). All samples were plotted according to their ImmnueScore and TMI, and prognostic and metabolic features are summarized (Figure ).

To verify the study results, we validated the association between TMI and immune landscape as well as the prognostic prediction with another publicly available independent gene expression dataset including 183 lung adenocarcinoma samples (see additional methods in Supplementary Material) 26. There was significant effect of cell type-based clusters on TMI (F = 10.7, p = 1.67 x 10-6; Figure ). Additionally, ImmuneScore was significantly different between cell type-based clusters (F = 72.6, p < 1 x 10-15; Figure ). From the TCGA data, C2 showed the lowest TMI and C3 showed the highest TMI. The ImmuneScore of C1 was lower and that of C4 was higher than those of the other clusters. The results of survival analysis were also consistent with the results from the TCGA data. C2 and C4 were associated with favorable prognosis (HR = 0.39, 95% C.I. 0.18-0.85 for C2 and HR = 0.56, 95% C.I. 0.31-1.00 for C4). High ImmuneScore demonstrated significantly favorable overall survival (HR = 0.60, 95% C.I. 0.38-0.98). High TMI was associated with a trend of poor overall survival though it did not reach statistical significance (HR = 1.42, 95% C.I. 0.88 - 2.28) (Figure ).

Discussion

As the TME and metabolic milieu serve critical roles in cancer immune escape leading to tumor growth and aggressiveness 8, 27, we sought to demonstrate the association between immune cell type composition in the TME and tumor metabolism. The key feature of the altered metabolism within cancer tissue is increased aerobic glycolysis, so-called 'Warburg Effect' 28, and it changes the TME to favor cancer cell growth. An important mechanism of cancer immune escape is glucose competition between the cancer cells and tumor-infiltrating lymphocytes 6, 7. As T cell activation requires upregulated glucose metabolism similar to cancer cells to permit proliferation and differentiation into effector T cells 29, the glucose competition between cancer cells and tumor-infiltrating lymphocytes can influence immune cell differentiation and composition in the TME. Our study demonstrated that a cluster with higher Treg cell scores (C3) showed tumor hypermetabolism, while a cluster with higher mast cell and CD4+ central memory T cell scores (C2) showed relative hypometabolism. Immune tolerance stimuli mediated by Treg cells drives metabolic reprogramming of immune systems to fatty acid oxidation instead of glycolysis 30. Thus, enriched Treg cells in the TME could result from relatively poor nutrient conditions due to elevated glycolysis by cancer cells. The previous study also reported that inhibition of glycolysis promoted Treg cell differentiation 31. This cluster with high Treg cell scores also showed relatively poor prognosis, which corresponded with the previous studies reporting its association with higher risk of recurrence 32, 33. On the other hand, the association of mast cell infiltration with favorable prognosis in various cancer types has been reported 34-36. Our study also showed relatively low glucose metabolism and favorable prognosis in the mast cell-enriched cluster, which corresponds with these previous clinical studies. Another cluster with a favorable prognosis (C4) demonstrated overall high immune cell enrichment scores including CD8+ T cells and macrophages. This result corresponded with the recent meta-analysis result that CD8+ T cell infiltration is associated with favorable prognosis in various cancer types 37. Of note, recent studies based on histopathologic and FDG PET data reported an inverse correlation between tumor-infiltrating lymphocytes and cancer glucose metabolism, which supports a close interconnection between the TME and cancer cells 38, 39. According to the systematic analyses of all immune cell types, the elevated glucose metabolism in tumor samples is not solely correlated with immune cell depletion in the TME. Instead, the association between immune cell composition and tumor metabolism suggested that the metabolic milieu leads to a broad spectrum of immune landscape dynamics. We suggest glucose deprivation in the TME may cause a compositional change of immune cells as well as cellular exhaustion and dysfunction so that it eventually creates an immunosuppressive microenvironment and could be closely related to patients' prognosis.

The tumor metabolism prediction model we developed in this study used the transcriptome data of the genes from two SUV-related modules as inputs (see additional data in Supplementary Material). Functional enrichment analysis demonstrated that both modules were functionally associated with extracellular matrix organization (Table ). It is now well known that the extracellular matrix plays a very important role in tumor initiation and progression; moreover, it also regulates cancer metabolism in many complicated ways 40. Focal adhesion signaling mediates the transmission of extracellular matrix signals into the tumor cells, which promotes the activation of certain pathways that increase glycolysis. Therefore, we think that the genes included in two modules are closely associated with the regulation of tumor metabolism, which could make our model estimate tumor metabolism with decent performance.

Our results could be applied to developing an integrative biomarker based on tumor metabolism and immune microenvironment. ImmuneScore was a predictive prognostic factor independent of other clinical features. TMI was also associated with prognosis though it did not reach statistical significance in the multivariate analysis. Nonetheless, maximum SUV in lung cancer has been regarded as an important prognostic imaging parameter 41, 42. These two axes, ImmuneScore and TMI, could divide LUAD according to immune landscape to some degree as shown in Figure . C1 and C4 could be clustered by low and high immune cell enrichment, respectively. Though C2 and C3 showed similar ImmuneScores, they could be divided by different metabolisms. A new pathologic classification system based on immune cell distribution and density was recently suggested as a potential predictive marker for cancer immunotherapy as well as a prognostic marker 43, 44. Moreover, the role of tumor metabolism in immunotherapy response evaluation and prediction has been increasingly emphasized 45, 46. As a future research, it needs to be investigated whether ImmuneScore, TMI, and the immune cell type clusters have predictive values for evaluating immunotherapy response. Monitoring dynamics of tumor metabolism and immune landscape along with immunotherapy will provide a clue for understanding immunotherapy response and resistance.

There are some limitations to our study. First, the design of this study was retrospective; thus, datasets that contain both matched transcriptome and FDG PET image data were limited. To overcome these data issues, we combined large transcriptome data without image data for developing the tumor metabolism estimation model. To clarify our results, further study with a well-controlled prospective design is warranted. Second, an additional analysis using FDG PET images and histopathologic data will make our results more solid and reliable. The key variables analyzed in this study were immune cell enrichment scores and TMI, and they were estimated based on mRNA transcriptome data. Even though our tumor metabolism prediction model showed robust performance and TMI demonstrated prognostic and clinical characteristics similar to FDG uptake reported in previous literature, comprehensive multimodal analysis including transcriptome, imaging, and histopathologic data can support our results. We developed TMI, which predicts tumor SUV; however, TMI can have different scales than real SUV according to independent external gene expression data. As shown in Figure , the model relatively underestimates SUV, though the output showed significant positive correlation. Thus, our TMI prediction model should be analyzed cautiously in other gene expression data and should be used for comparison between subjects from the same experimental batch.

In this study, we found the association between the immune landscape of the TME and tumor metabolism. The immune cell composition in the TME was heterogeneous across LUAD samples. It was closely associated with tumor metabolism and also patient prognosis. As tumor metabolism was elevated, immune cell subpopulations were prone to enrichment of immunosuppressive cellular subtypes. ImmuneScore and TMI were associated with patient prognosis. Additionally, the comprehensive analysis of the immune landscape of the TME and tumor metabolism may support development of an integrated biomarker for immunotherapy responses as well as prognosis.

Supplementary figures and tables.

Click here for additional data file.

Table 1

Multivariate Cox regression analysis for overall survival

		Univariate analysis			Multivariate analysis
Variables		Hazard Ratio (95% CI)	p-value		Hazard Ratio (95% CI)	p-value
Age	Older than 65	1.13 (0.83 - 1.53)	0.45
Stage	≥III	2.62 (1.91 - 3.60)	<0.0001		2.31 (1.67 - 3.20)	<0.0001
TMI	Hypermetabolic tumor	1.61 (1.19 - 2.19)	0.002		1.23 (0.85 - 1.77)	0.27
Microenvironment cell type clusters	C2	0.58 (0.38 - 0.89)	0.01		0.63 (0.40 - 1.01)	0.05
	C3	0.90 (0.59 - 1.37)	0.61		0.83 (0.53 - 1.29)	0.41
	C4	0.62 (0.38 - 1.01)	0.05		0.82 (0.48 - 1.38)	0.45
ImmuneScore	High ImmuneScore	0.65 (0.48 - 0.88)	0.005		0.69 (0.48 - 0.98)	0.04
Molecular Subtypes	Magnoid	1.99 (1.35 - 2.92)	0.0005		1.42 (0.93 - 2.18)	0.11
	Squamoid	1.84 (1.29 - 2.62)	0.0008		1.46 (0.96 - 2.22)	0.08

CI: confidence interval