Exploring the Interplay between Metabolism and Tumor Microenvironment Based on Four Major Metabolism Pathways in Colon Adenocarcinoma

Tumor metabolism plays a critical role in tumor progression. However, the interaction between metabolism and tumor microenvironment (TME) has not been comprehensively revealed in colon adenocarcinoma (COAD). We used unsupervised consensus clustering to establish three molecular subtypes (clusters) based on the enrichment score of four major metabolism pathways in TCGA-COAD dataset. GSE17536 was used as a validation dataset. Single-cell RNA sequencing data (GSE161277) was employed to further verify the reliability of subtyping and characterize the correlation between metabolism and TME. Three clusters were identified and they performed distinct prognosis and molecular features. Clust3 had the worst overall survival and the highest enrichment score of glycolysis. 86 differentially expressed genes (DEGs) were identified, in which 11 DEGs were associated with favorable prognosis and 75 DEGs were associated with poor prognosis. Striking correlations were observed between hypoxia and glycolysis, clust3 and hypoxia, and clust3 and angiogenesis (P < 0.001).We constructed a molecular subtyping system which was effective and reliable for predicting COAD prognosis. The 86 identified key DEGs may be greatly involved in COAD progression, and they provide new perspectives and directions for further understanding the mechanism of metabolism in promoting aggressive phenotype by interacting with TME.


Introduction
Colon adenocarcinoma (COAD) is one of the most diagnosed cancer types worldwide, which represented 6.0% of new cancer cases and was the cause of 5.8% of new cancer deaths in 2020 according to global cancer statistics [1]. Risk factors of COAD are various, which can be generally classified into unhealthy lifestyles (such as smoking, heavy alcohol, and unhealthy diet), disease-caused or drug-induced (such as long-term androgen deprivation therapy, diabetes mellitus and insulin resistance, acromegaly, and renal transplantation under long-term immunosuppression), and family history [2]. Most of COAD individuals are sporadic, about 70%, and ages above 50 years are common [3]. Although screening programs such as colonoscopy, fecal immunochemical test (FIT), and computed tomography (CT) colonography are recommended and popularized for high-risk populations in the recent decades, efficient and personalized therapeutics are still current challenges especially for metastatic patients.
Immunotherapy or immune checkpoint blockade raises the hope for metastatic patients. Targeted therapies such as vascular endothelial growth factor (VEGF) inhibitors and epidermal growth factor receptor (EGFR) inhibitors are examined to effectively hinder angiogenesis and thus impede tumor growth [4]. Immunotherapy and targeted therapy regulate tumor microenvironment (TME) and inhibit tumor progression through targeting key proteins associated with poor prognosis. Tumor metabolism is a pivotal factor in TME modulation, where four critical metabolism pathways, glycolysis, pentose phosphate pathway (PPP), fatty acid oxidation (FAO), and glutaminolysis, are reported to be strikingly involved in tumor metabolic reprogramming contributing to TME modulation and tumor progression [5][6][7][8][9]. Particularly, a number of glycolysis-related prognostic biomarkers have been identified in various cancer types [10][11][12][13]. Nevertheless, the mechanism of the interaction between TME and metabolism pathways has not been well understood.
An increasing number of bioinformatics tools accelerate the development of understanding potential mechanisms of cancer development and the novel prognostic biomarkers for predicting cancer prognosis [14][15][16][17]. Moreover, the development of single-cell RNA sequencing (scRNA-seq) technology allows more accurately parsing the interplay or crosstalk of various pathways in TME. In this study, we sought to reveal the interaction between metabolism pathways and TME based on expression data and single-cell data gleaned from public database. ree molecular subtypes (clusters) were established based on the expression of genes within metabolism pathways. ree clusters showed distinct overall survival and expression patterns of metabolism. We identified key differentially expressed genes (DEGs) between two groups with distinct prognosis. Dramatic correlation between TME hallmarks (angiogenesis and hypoxia) and metabolism pathways in malignant cells were illustrated by single-cell analysis. Our study lays a foundation for further research on the mechanism of tumor metabolism in colon cancer.

Data Preprocessing.
For TCGA-COAD dataset, samples without survival information were excluded, and survival time within 30 days to 10 years was restricted for all samples. 410 samples with expression data remained in TCGA-COAD dataset. For GSE17536 dataset, probes in chip files were converted to gene symbols. Samples were screened according to the same conditions as TCGA-COAD dataset. 170 samples with expression data remained in GSE17536 dataset.

Unsupervised Consensus
Clustering. Firstly, single sample gene set enrichment analysis (ssGSEA) in GSVA R package was applied to calculate ssGSEA score of four metabolism pathways for each sample in TCGA-COAD dataset [21]. en unsupervised consensus clustering was used to cluster samples based on the ssGSEA score of four metabolism pathways through ConsensusClusterPlus R package (v1.54.0) [22]. Kmdist algorithm and "Pearson" distance were used to conduct 500 bootstraps with each bootstrap containing at least 80% samples of TCGA-COAD dataset. Cluster number k was set from 2 to 10, and optimal k was determined according to cumulative distribution function (CDF) and area under CDF curve. e clustering was validated in GSE17536 dataset.

Mutation Analysis.
Copy number variations (CNVs) and single nucleotide variations (SNVs) in TCGA-COAD dataset were analyzed separately. For CNV data, Genomic Identification of Significant Targets in Cancer (GISTIC) 2.0 software was employed to detect arm-level and focal CNVs based on hg38 reference genome [23]. e number of copies greater than 1 was considered the threshold of copy amplification and that less than −1 was considered the threshold of copy loss. R package of maftools was used to assess SNVs [24].

Tumor Microenvironment Analysis.
To evaluate TME features in different clusters, several methodologies were applied. MCPcounter R package and the Proportion of Immune and Cancer cells (EPIC) methodology were used to calculate the enrichment score of immune cells [25,26]. Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) tool was introduced to evaluate immune infiltration and stromal infiltration by calculating immune score, stromal score, and ESTIMATE score based on signatures [27]. Gene sets of macrophage-related pathways including toll-like receptor signaling, natural killer cell mediated cytotoxicity, and antigen processing and presentation were obtained from MSigDB. Signatures of IFN-c and cytotoxicity were obtained from previous research [28,29].

Differential Analysis between Clust12 and Clust3.
DEGs between clust12 and clust3 were detected by limma R package under conditions of |log2 (fold change)| > log2 (1.5) and false discovery rate (FDR) < 0.05 [30]. String (https://cn. string-db.org/) online tool was applied to identify interactions among DEGs, and the interaction network was visualized by CytoScape (3.8.0) [31]. Degree of each DEG was analyzed by Analyze Network tool in CytoScape. Higher degree represents greater importance of the gene in the network. Annotation of gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was outputted by WebGestaltR package [32]. FDR < 0.05 was determined to screen significant GO terms and P < 0.05 was determined to screen significantly enriched KEGG pathways.

2
Journal of Oncology

Identifying Immune Cells in scRNA-Seq Data.
GSE161277 dataset containing scRNA-seq data was included for validating the reliability of molecular subtyping based on four metabolism pathways. Within the dataset, 13 samples, 4 patients, and 5 tissues (adenoma, blood, carcinoma, normal, and paracancer tissues) were included. Seurat R package was implemented to process scRNA-seq data [33]. Firstly, single cells were screened under conditions of each gene expressing at least in three cells and each cell expressing at least 250 genes. en PercentageFeatureSet function was used to calculate the percentages of mitochondria and rRNA. Single cells with mitochondria <35% and unique molecular identifiers (UMI) > 100 within one cell were further screened. Log-normalization was conducted for normalizing single-cell data. FindVariableFeatures function was performed to identify highly variable genes based on "vst." ScaleData function was conducted to scale the data and principle component analysis (PCA) was used to reduce the dimensionality.
en FindNeighbors and FindClusters functions were implemented to cluster cells when dim � 40 and resolution � 0.9. CD45 (PTPRC) was set as the marker to identify immune cells. Finally, 26188 immune cells remained.

Annotating Cell Types in Processed scRNA-seq Data.
e expression data of 26188 immune cells was rescaled and highly variable genes were recalculated. Cells were clustered under dim � 35 and resolution � 0.3 by FindNeighbors and FindClusters functions in Seurat R package [33]. en RunTSNE function was performed to reduce dimensionality for 26188 ells and identify cell subgroups with different samples, tissues, and patients. Classical immune markers were used to annotate cell subgroups. Differential genes were screened by using FindAllMarkers function with parameters of logfc � 0.5 and Minpct � 0.35 (P < 0.05). Only the top 5 differential genes were visualized. clusterProfiler R package was applied to annotate enriched KEGG pathways for differential genes [34].

Constructing Molecular Subtypes Based on Metabolism-Related Pathways.
e workflow of this study is shown in Supplementary Figure S1. For each sample in TCGA-COAD dataset, we calculated the enrichment score of four metabolism pathways, glycolysis, PPP, FAO, and glutaminolysis, by using ssGSEA (Supplementary Table S1). en unsupervised consensus clustering was applied to cluster 410 samples based on the enrichment score of four metabolism pathways. Cluster number k was set from 2 to 10 in CDF analysis. Cluster number k � 3 was confirmed as optimal according to CDF and area under CDF curve (Figures 1(a) and 1(b)). Consensus matrix showed that samples were clearly divided into three clusters or molecular subtypes (Figure 1(c)). ree clusters showed different patterns of the enrichment of four metabolism pathways (Figure 1(d)). Glycolysis was the most enriched in clust3, while clust1 had the highest enrichment of glutaminolysis. e enrichment scores of PPP, FAO, and glutaminolysis were significantly higher in clust1 than in clust3. Over half of genes within four pathways (44 of 74) were differentially expressed among three clusters (Supplementary Figure S2).
Kaplan-Meier survival analysis revealed that three clusters had differential overall survival (P � 0.0077, Figure 1(e)). Using the same methodology in GSE17536 dataset, 170 samples were divided into three clusters with distinct overall survival (P � 0.0049, Figure 1(f )). Clust1 had the worst prognosis, while clust3 had the longest overall survival in both datasets (Figures 1(e) and 1(f )). PCA presented the different distribution of three clusters in two datasets (Supplementary Figure S3). Of the distribution of three clusters in different clinical features, we observed no significant difference in ages, genders, and stages (Supplementary Figure S4).

Differential Genomic Features among ree Subtypes.
To know the genomic features of three subtypes, we applied gistic2 software to analyze the CNV data and visualized CNVs of 22 chromosomes (Figures 2(a)-2(c)). ree clusters showed obvious gain and loss of CNVs where gain of CNVs occurred largely in chromosomes 7, 8, 13, and 20, and loss of CNVs occurred largely in chromosome 18. Although similar CNV patterns were shown in three clusters, still a number of significantly differential CNVs were identified among them (Supplementary Figure S5). Clust2 had the most number of significantly amplified CNVs especially in chromosomes 17 and 20, compared to clust1 and clust3. Relatively lower counts of loss of CNVs were found in clust1 and clust2, but clust3 had more loss of CNVs than gain of CNVs itself (Supplementary Figure S5). In addition, we screened the top 15 mutated genes within four metabolism pathways based on SNV data ( Figure 2(d)). e mutation frequencies ranged from 2% to 5%, and missense mutations contributed the majority of mutations. GLUD2 was the top mutated gene, where missense mutations contributed the most of SNVs.
3.3. Differential TME among ree Clusters. Next we used MCPcounter to evaluate the estimated proportion of 22 immune-related cells in three clusters in TCGA-COAD dataset. 9 of 22 immune-related cells showed a significant difference on the proportion among three clusters (Figure 3(a)). In particular, resting memory CD4 T cells and M0 macrophages consisted of a relatively high proportion among 22 immune-related cells. Clust1 had the highest proportion of resting memory CD4 T cells, while clust3 had the highest proportion of M0 macrophages (P < 0.01, Figure 3 Figure S6). Interestingly, clust1 with the longest overall survival had the least immune infiltration, and clust3 had the highest score of stromal and immune score (P < 0.0001, Figure 3(b)). We then assessed the enrichment of 10 oncogenic pathways [35] and observed that 6 of 10 oncogenic pathways were differentially enriched among three clusters (P < 0.01, Figure 3(c)). Hippo, Notch, RAS, and Wnt signaling pathways were more activated in clust3, while NRF1 and TP53 were more activated in clust1. As M0 macrophages were identified to be significantly differentially enriched among three clusters, we further assessed the immune regulation related to macrophages. We selected three pathways related to macrophages from MSigDB including antigen processing and presentation, toll-like receptor signaling pathway, and natural killer (NK) cell mediated cytotoxicity. GSEA revealed that the three pathways were all most activated in clust3 (P < 0.01, Figures 3(d)-3(f)), which was accordant to the highest proportion of macrophages in clust3. In addition, clust3 also displayed high enrichment of interferon-c (IFN-c) and T cell cytotoxicity (CYT) (P < 0.05, Figures 3(g) and 3(h)). In GSE17536 dataset, similar results were presented (Supplementary Figure S7). Overall, clust3 showed more significantly active immune response than clust1 and clust2. Differential TME among three clusters indicated that metabolism pathways were possibly involved in complicated TME modulation.

Identifying Differentially Expressed Genes among Clusters.
Given that clust1 and clust2 had superior prognosis than clust3, we combined the data of clust1 and clust2 named clust12. By comparing expression profiles between clust12 and clust3 in TCGA-COAD and GSE17536 datasets, we identified a number of DEGs between two groups under conditions of |log2 fold change (FC)| > log2 (1.5) and FDR <0.05 (Figures 4(a) and 4(b)). 11 upregulated and 75 downregulated genes (clust12 versus clust3) were screened in both datasets (Figure 4(c)). Furthermore, based on 86 genes, we used String (https://cn.string-db.org/) to screen their interactions and found 69 interacted genes visualized by CytoScape (Figure 4(d)). Most of the genes were downregulated and only two genes (ADH1C and UGT2A3) within 69 interacted genes were upregulated. Genes with degree >10 were listed (Supplementary Table S2). Subsequently, we conducted GO and KEGG analysis on 69 genes. Within GO analysis, 92 terms of biological process, 20 terms of molecular function, and 13 terms of cellular component were annotated (Supplementary Table S3, FDR < 0.05). e top 10 enriched terms were visualized ( Supplementary  Figures S8(a)-S8(c)). In addition, 9 KEGG pathways were significantly enriched such as ECM-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway (FDR < 0.05, Supplementary Figure S8(d)).

Identifying Eight Cell Types Based on scRNA-Seq Data.
To further verify the reliability of our molecular subtyping based on four metabolism pathways, we introduced scRNA-seq data of COAD and sought to distinguish malignant and nonmalignant cells. Single-cell data was preliminarily processed and screened to meet the standards of each gene expressed at least in three cells, each cell expressing at least 250 genes, each cell containing less than 35% mitochondria, and UMI > 100 for each cell (Supplementary Figures S9(a)-S9(c)). 45892 single cells remained. en PCA was applied to diminish the dimensionality of the data. Cells were clustered into 34 clusters (resolution � 0.9) and CD45 marker was used to identify immune cells with a total number of 26188 ( Supplementary Figures S9(d) and S9(e)). Subsequently, 26188 immune cells were further clustered under resolution � 0.3. 12 subgroups were identified, and t-SNE plots grouped by different samples, patients, and tissues were shown (Figures 5(a)

Malignant Cells Were More Enriched in Clust3.
As the complicity of tumor tissues that may contain normal cells, we employed CopyKat to estimate genomic copy number profiles for more accurately distinguishing malignant and nonmalignant cells (Supplementary Figure S10) [36]. As a result, 8055 malignant cells and 18133 normal cells were distinguished. In the previous section, we identified 11 upregulated genes and 75 downregulated genes by comparing clust12 with clust3 (Figure 4(c)). We used ssGSEA to calculate their enrichment score of upregulated and downregulated genes in malignant and nonmalignant cells, respectively. Significantly differential enrichment of these genes was exhibited in two groups, with the 11 DEGs higher expression in nonmalignant cells and the 75 DEGs higher expression in malignant cells (P < 0.0001, Figure 6). e result demonstrated that the molecular subtyping based on metabolism pathways was reliable, and these dysregulated genes may be highly associated with the regulation of four metabolism pathways.

Metabolism of Malignant Cells Was Associated with TME.
To understand if TME made a difference on the metabolism of malignant cells, we introduced hypoxia and angiogenesis as indicators. Here hypoxia score was calculated by ssGSEA based on genes within hypoxia pathway, and angiogenesis score was calculated based on the signature from Masiero et al. [37]. Pearson correlation analysis revealed positive correlations between either two of the four metabolism pathways in malignant cells, especially between glycolysis and PPP (Figure 7(a)). In addition, we applied Mantel test to evaluate the relation of angiogenesis and hypoxia with four metabolism pathways. We observed that angiogenesis score was negatively correlated with four pathways, while positive correlation was shown between hypoxia and four pathways 6 Journal of Oncology       Journal of Oncology 9 (Figure 7(a)). Notably, significantly high correlation coefficient was observed between hypoxia and glycolysis score (R � 0.62), indicating the strong interaction between two pathways. Table S4). Notably, all detected cell types were more enriched in clust12 compared with clust3 (P < 0.0001, Figure 7(b)). Furthermore, we found that clust12 score was negatively correlated with hypoxia (R � −0.46), glycolysis (R � −0.58), PPP (R � −0.51), and FAO (R � −0.35) particularly (Figure 7(c)). Conversely, clust3 score was positively correlated with hypoxia and angiogenesis especially (R � 0.50 and 0.40, respectively, Figure 7(c)). In addition, we found that the enrichment score of hypoxia was significantly differential among three subtypes, where clust3 had the highest enrichment score (Supplementary Figure S11), which was consistent with the above result. Hypoxia and angiogenesis are important TME hallmarks in cancer. e results suggested that metabolism regulation of malignant cells was possibly regulated by hypoxia and angiogenesis, as well as TME. Additionally, the significantly differential distribution of immune cells between clust12 and clust3 proved that 86 markers (DEGs) were effective and reliable for distinguishing two groups in single-cell data.

Discussion
Tumor metabolism has a major impact in TME modulation [38,39], but the link between tumor metabolism and TME has not been clearly understood in COAD so far. In this study, we dug out key genes interacting with both four major metabolism pathways and TME by analyzing expression data and single-cell data. Initially, three clusters were established through unsupervised consensus clustering based on the enrichment score of four metabolism pathways. Given that clust1 and clust2 had similar overall survival, we combined the expression data of two clusters. en 86 important DEGs were identified between clust12 and clust3, and 69 of 86 DEGs presented close interactions through PPI analysis.
We estimated that these DEGs may play a critical role in tumor metabolism and have a pronounced impact on COAD prognosis. As demonstrated in the single-cell analysis, malignant cells and nonmalignant cells had distinct expression levels of 86 DEGs (11 as clust12 markers and 75 as clust3 markers). Malignant cells had extremely high expression level of 75 DEGs associated with poor prognosis in clust3, which illustrated the important role of 75 DEGs in tumor progression, as well as the effectiveness and reasonability of subtyping based on four metabolism pathways.
Besides glycolysis, PPP was also strongly correlated with hypoxia (R � 0.42), and clust12 score was shown to be negatively associated with PPP (R � −0.51). PPP is one of major pathway for glucose catabolism, and it is reprogrammed accompanied with the activation of glycolysis in tumor cells [49]. In the human colon cancer cell line (HT29), enzymes within PPP pathway are enhanced during cell cycle progression [50], suggesting that PPP is a potential target for inhibiting tumor cell proliferation [51]. High expression of 11 marker genes in clust12 may play an important role for alleviating PPP activation.
Angiogenesis is more activated in clust3 revealed from a correlation coefficient of 0.40 between clust3 score   and angiogenesis. Several targeted drugs display antiangiogenesis effect in colon cancer, such as regorafenib [52], anti-interleukin-6 receptor antibody [53], and SARI [54]. Angiogenesis is regulated by innate immune cells such as TAMs, MDSCs, tumor-associated neutrophils (TANs), mast cells, and NK cells, which produce proangiogenic factors [55], and promote immunosuppressive TME [56]. Other two metabolism pathways, FAO and glutaminolysis, exhibited a relatively weak correlation with either hypoxia or angiogenesis,    indicating that they may be less involved in dysregulated immunity in COAD.

Conclusions
In conclusion, this study depicted a crosstalk between tumor metabolism and TME features for COAD through integrated analysis of expression data and scRNA-seq data. ree molecular subtypes were constructed based on metabolismrelated genes, and they performed significant differences in metabolism, genomic variations, immune infiltration, and enriched pathways, suggesting that the four metabolism pathways may regulate TME and thus affect COAD progression. Importantly, we identified 86 differentially expressed genes from three subtypes that were illustrated in scRNA-seq data to be highly associated with malignancy, hypoxia, and glycolysis. We envisioned that these key genes may provide new directions for further revealing the mechanism of metabolism in TME modulation.

Conflicts of Interest
e authors declare no potential conflicts of interest.

Supplementary Materials
Supplementary Figure S1: the flowchart of this study. Supplementary Figure S2: expression of genes within four metabolism pathways in three clusters. Kruskal test was conducted. ns, no significance. * * P < 0.01, * * * P < 0.001, and * * * * P < 0.0001. Supplementary Figure S3: the distribution of three clusters in PCA plots. Supplementary Figure S4