1. Introduction
As one of the most common digestive tract cancers, annual colorectal cancer incidences and fatalities could reach over 1.9 million and 940,000, respectively, accounting for 10% of all cancer cases and deaths.[1] Over the past few decades, physicians seem to be accustomed to describing malignancies in the colon and rectum as colorectal cancer. However, a growing body of research has shown that colon adenocarcinoma and rectum adenocarcinoma (READ) have apparent differences in developmental origin, genetic characteristics, gut microbiota, and risk factors.[2–4] In addition, the incidence of READ in Asia is much higher than that of colon adenocarcinoma. The proportion of early-onset READ cases (diagnosed under 50) increases year by year, resulting in a significant public health concern worldwide.[5,6] Smoking cessation, a healthy diet, and proper exercise can prevent READ. At the same time, colonoscopy screening, stool testing, and other methods can be conducive to the early diagnosis of READ.[1–3,7] Nevertheless, the prognosis of READ patients has not substantially improved due to the high cost of colonoscopies, inadequate medical services, and recurrent diseases.[1,2] Therefore, reliable and reasonable indicators must be identified to assess individual prognostic risk and guide individualized treatment so that the survival of READ patients can be improved.
Over the past decade, continuous efforts have been invested in exploring READ screening strategies and prognostic markers. The proposal and application of the TNM international staging system, C-reactive protein, carcinoembryonic antigen, and other methods have significantly improved READ diagnosis and treatment.[8,9] Notably, these screening strategies and prognostic markers have numerous limitations, and novel effective tools are required to complement and replace them. Fortunately, the continuous accumulation of high-throughput data has provided a deeper and more comprehensive understanding of the genetic characteristics and pathogenesis of READ, making it possible to develop READ prognostic markers based on sequencing data.[10–12]
More and more evidence has shown that fatty acid metabolism profoundly affects cancer occurrence, development, and treatment resistance through adenosine triphosphate generation (beta-oxidation), glycerophospholipid synthesis, protein acylation, and sustained redox homeostasis.[2,3] One clinical study showed significantly varied medium-chain fatty acid concentrations among individuals in the READ progression stage and identified capric acid as a highly effective READ biomarker.[13] Meanwhile, Zhu et al[14] demonstrated that blocking the fatty acid metabolizing enzyme could reduce lipid metabolism in mice and inhibit the occurrence of READ. Therefore, in-depth exploration of fatty acid metabolism in READ is of great significance for identifying potential targets of READ and improving the treatment.
In this study, based on publicly available data in the cancer genome atlas (TCGA) and gene expression omnibus (GEO), we distinguished into fatty acid catabolic and anabolic subtypes through a comprehensive analysis of fatty acid metabolism genes in READ patients, in which patients with fatty acid catabolic subtype had not only a poorer prognosis but also a relatively low proportion of myeloid dendritic cells (mDCs) infiltration within the tumor microenvironment. Subsequently, we applied least absolute shrinkage and selection operator (LASSO) regression and multifactorial Cox regression and integrated the computational results of different cohorts to finally construct a highly specific and sensitive risk prognostic model consisting of 3 fatty acid metabolism genes. The research flow of this paper is shown in Figure 1.
;)
Figure 1.: Flowchart of the research workflow. READ = rectum adenocarcinoma, DEGs = differentially expressed genes, GSEA = gene set enrichment analysis, LASSO = least absolute shrinkage and selection operator.
2. Materials and methods
2.1. Data collection and preprocessing
RNA-seq, microarray data, and corresponding clinical information were obtained from TCGA data portal (https://portal.gdc.cancer.gov/) and the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). Raw counts were normalized to transcripts per million by gene length in the annotation file (GDC.h38 Flattened GENCODE v22 GFF). Then, log2-scale transformed microarrays data were not further normalized. READ patients were screened from the dataset to exclude individuals whose RNA expression did not match clinical information. In addition, we also excluded individuals with genetic mutations from the microarray data. A total of 123 TCGA-READ patients and 109 READ patients from GSE87211[15] were included in the study cohort, and the clinical characteristic is presented (see Table 1).
Table 1 -
Baseline clinical characteristic.
|
TCGA:READ (n = 123) |
GSE87211 (n = 109) |
| Age (yr) |
65 (32–88) |
64 (36–82) |
| Gender |
|
|
| Male |
72 (59%) |
80 (73%) |
| Female |
51 (41%) |
29 (27%) |
| Stage |
|
|
| I/II |
57 (46%) |
34 (31%) |
| III/IV |
66 (54%) |
75 (69%) |
| Event |
|
|
| Alive |
114 (93%) |
82 (75%) |
| Dead |
9 (7%) |
27 (25%) |
| Time (mo) |
4 (0–122) |
63 (0–151) |
READ = rectum adenocarcinoma, TCGA = the cancer genome atlas.
2.2. Consensus clustering
The “HALLMARK_FATTY_ACID_METABOLISM” genes in the molecular signatures database [16] were located and defined as fatty acid metabolism genes. After excluding genes with <10% expression, matrix constructed from fatty acid metabolism genes were used for clustering in Consensus ClusterPlus (v.1.58.0; parameters: reps = 500, pItem = 0.8, pFeature = 1, distance = “spearman”).[17] The best clustering result under different cluster numbers was determined by comparing the consensus cumulative distribution function and Delta Area. The t-distributed stochastic neighbor embedding was executed in Rtsne (v.0.15) to validate the accuracy of the cluster. Finally, clusters were visualized as heatmaps in ComplexHeatmap (v.2.10.0).[18]
2.3. Gene set enrichment analysis (GSEA)
The gene expression level of RNA-seq data (raw counts) and Microarray data were analyzed using the DEseq2 (v.1.34.0) [19] and limma (v.3.50.0),[20] respectively. The genes were ranked sequentially from highest to lowest based on Log2FoldChange for GSEA in clusterProfiler (v.4.2.2),[21] with the signaling pathways with a P value < .05 were considered significant and visualized.
2.4. Metabolic signature identification
Extract the “FATTY ACID METABOLISM” related gene sets from the 225 eigengene dataset built in the IOBR package.[22] Subsequently, based on consensus clustering results, the signature score of different clusters against gene sets was evaluated by a single-sample gene set enrichment analysis,[23] and minimal and maximal values were set to 0 and 1, respectively. In addition, Student t test was applied to analyze whether there were significant differences in the signature score between groups.
2.5. Decoding the tumor microenvironment
Based on the expression matrix, we quantitatively scored and normalized the abundance of immune and stromal cell populations in READ patients via MCP Counter[24] in IOBR. Furthermore, the differences in the relative abundance of the same cell population in different clusters were compared, and the correlation between individual cell abundance and the fatty acid metabolic signature score was analyzed by calculating the Pearson correlation coefficient.
2.6. Kaplan–Meier survival analysis
The matrix containing outcome events and follow up times corresponding to the patient were collated by cluster, and the follow up times were adjusted into monthly units. Kaplan–Meier survival analysis with grouping as categorical information was conducted in survival (v.3.2-13) and presented using survminer (v.0.4.9).
2.7. LASSO regression
The glmnet (v.4.1-3) was applied to the matrix, and the cluster of the patient was used as a regression condition. Variable selection was presented to penalize the data fitting and to calculate both the minimum and standard errors (1se) of the optimal parameters. A model was constructed from the screened variables, and its predictive ability was evaluated by receiver operating characteristic (ROC) curves as well as the area under curve (AUC) in ROCR (v.1.0-11).
2.8. Multivariate cox regression and risk prognostic model construction
To reduce the complexity and increase the interpretability of the model, multivariate Cox regression was applied to the included variables for stepwise regression. The regression coefficients and expression levels corresponding to each variable were then combined to obtain the risk score formula as follows.
Accordingly, patients were divided into 2 groups (high risk and low risk), and the survival statuses of the 2 groups as well as the distribution of each variable’s expression levels were visualized. In addition, we stratified the patients according to age, gender, and tumor stage in turn, thus determining the distribution of individuals in the 2 groups at different levels. Finally, the prognostic risk prediction ability of the model was validated through survival analysis.
2.9. Assessing the independent prognostic value of the model
To more comprehensively assess the independent prognostic value of the model, we calculated risk scores for each patient and integrated them with age, gender and tumor stage for univariate COX regression analysis. In addition, we also analyzed the accuracy of the model based on the ROC curve.
2.10. Statistics and visualization
The statistics information for every experiment was detailed in the figures, figure legends, and figure notes. The violin plot and scatter plot were calculated and visualized by ggstatsplot (v.0.9.1). The other portraits (such as heatmap, the volcano plot, and so on) were calculated by rstatix (v.0.7.0) and represented via ggplot2 (v.3.3.5) or ggtree (v.3.2.1).[25] For the split violin plot, Student t test between 2 categorical variables were performed using the Bonferroni test. The Pearson correlation coefficient was used to check the association between 2 continuous variables in the scatter plot. Log Rank test was used to assess survival differences between groups. For repeated statistical methods or color schemes in the paper, only the first occurrence is annotated, with the following occurrences consistent with the first annotation.
3. Results
3.1. Patient classification based on fatty acid metabolism genes
According to the expression level transcripts per million of fatty acid metabolism genes, we employed consensus clustering under different cluster numbers for the TCGA cohort comprising 123 READ patients (Fig. 2A; Figure S1 A–D, Supplemental Digital Content, https://links.lww.com/MD/I613). Based on cumulative distribution function (Fig. 2B) and Delta Area (Fig. 2C), K = 2 was selected as the best cluster. Two clusters were obtained by t-distributed stochastic neighbor embedding (Fig. 2D), and fatty acid metabolism genes were significantly differently expressed in the 2 clusters (Fig. 2E). The results are consistent with those in the GEO cohort (Figure S1 E–H, Supplemental Digital Content, https://links.lww.com/MD/I613 and Figure S2, Supplemental Digital Content, https://links.lww.com/MD/I614).
Figure 2.: Patient classification in the TCGA cohort. (A) Consensus clustering results (K = 2). (B) The consensus cumulative distribution function (K = 2). (C) The delta area (K = 2). (D) Distribution of cluster A and cluster B in t-SNE. (E) The expression levels of fatty acid metabolism genes between cluster A and cluster B. TPM = transcript per million, t-SNE = t-distributed stochastic neighbor embedding, TCGA = the cancer genome atlas.
3.2. The association between patient prognosis and fatty acid metabolic subtypes
To further explore the potential biological phenomena between cluster A and cluster B, we performed GSEA for different clusters in the TCGA cohort and GEO cohort separately based on the ranked gene sets. Cluster A was significantly enriched in pathways related to fatty acid catabolism, including Fat digestion and absorption, and PPAR signaling pathway (Fig. 3A–B). Cluster B is associated with fatty acid anabolism, such as Terpenoid backbone biosynthesis pathways and Steroid hormone biosynthesis pathways (Fig. 3A–B). Furthermore, in the metabolic signature identification for different cohorts, cluster A scored significantly higher than cluster B in metabolic pathways such as fatty acid beta-oxidation and fatty acid catabolic process, whit the latter preferring fatty acid biosynthesis (Fig. 3C–F). Accordingly, we defined cluster A as the fatty acid catabolic subtype and cluster B as the fatty acid anabolic subtype. More importantly, we found that patients attributed to cluster A had a poorer prognosis (P < .05) (Fig. 3G–H).
Figure 3.: Fatty acid metabolic characteristics and prognosis of patients between different clusters. (A–B) Enriched pathways among different clusters in the TCGA cohort and the GEO cohort (cluster A compared to cluster B). (C–F) Results of metabolic signature identification in the TCGA cohort and the GEO cohort. Component comparison using Student t test and the Bonferroni test (*, P < .05, **, P ≤ .01, ***, P ≤ .001). (G–H) The K-M curve among different clusters in the TCGA cohort and the GEO cohort (TCGA: P = .031; GEO: P = .045; Log Rank test). GEO = gene expression omnibus, K-M = Kaplan–Meier, TCGA = the cancer genome atlas.
3.3. The composition of the tumor microenvironment in different fatty acid metabolic subtypes
Based on the MCPcounter[22] in IOBR, we found that the relative abundance of immune and stromal cells in the tumor microenvironment differed in READ patients with different fatty acid metabolism. Among them, cluster A individuals had lower relative abundance of CD8+ T cells, mDCs, NK cells, endothelial cells, or fibroblasts than cluster B. Notably, the group differences with just mDCs were consistent between the 2 cohorts (Fig. 4A–B; Figure S3A–B, Supplemental Digital Content, https://links.lww.com/MD/I615). In addition, the relative abundance of individual mDCs was significantly negatively correlated with the signature score of fatty acid catabolic pathways including fatty acid beta-oxidation, fatty acid catabolic processes, and significantly increased along with the signature score of fatty acid biosynthesis (Fig. 4C–I; Figure S3C–1I, Supplemental Digital Content, https://links.lww.com/MD/I615).
Figure 4.: Composition of the tumor microenvironment among different clusters in the TCGA cohort. (A–B) Relative abundance of 8 immune and 2 stromal cell populations in cancer tissues. Component comparison using Student t test and the Bonferroni test (NS, P > .05; *, P < .05, **, P ≤ .01). (C–I) Correlation between relative abundance of mDCs and fatty acid metabolic signature score. Pearson correlation coefficient was used to check the association. The gray line is the fitted curve. mDCs = myeloid dendritic cells, TCGA = the cancer genome atlas.
3.4. The correlation of expression of aquaporin 7(AQP7), X inactive specific transcript (XIST), interleukin 4 induced 1 (IL4I1) with fatty acid metabolic subtypes and patient prognosis
In this study, we first identified signature genes (min:28, 1se:16) that could distinguish fatty acid metabolic subtypes of patients by LASSO regression, and constructed a model based on regression coefficients of signature genes. The ROC indicated that the prediction accuracy of the model was acceptable (Fig. 5A–D). Next, we used the selected signature genes to perform multivariate Cox regression with individual prognosis as the effect variable and integrated duplicated results from both cohorts, ultimately identifying 3 key genes, namely, AQP7, XIST, IL4I1, that could predict patient prognosis. The prognostic risk score formula for each patient was as follows: Risk score = (0.458 × the expression value of AQP7) + (−0.117 × the expression value of XIST) + (0.247 × the expression value of IL4I1). High expression of AQP7 (hazard ratio, HR = 2.064 (1.4408–4.5038); P < .01) and IL4I1 (HR = 1.34 (1.13–1.59); P = .034) in READ patients predicted poorer survival while XIST (HR = 0.5216 (0.3758–0.7564); P = .045) predicted better prognosis (Fig. 5E–F).
Figure 5.: Signature genes based on fatty acid metabolic subtypes and prognosis of patients in the TCGA cohort and the GEO cohort. (A–B) Selection of the optimal parameter (λ) (TCGA: left; GEO: right). (C–D) The ROC curve of signature genes. (E–F) The forest plot of AQP7, XIST, and IL4I1 (*, P < .05, **, P ≤ .01). ROC = receiver operating characteristic. AQP7 = aquaporin 7, GEO = gene expression omnibus, IL4I1 = interleukin 4 induced 1, TCGA = the cancer genome atlas, XIST = X inactive specific transcript.
3.5. Risk prognostic model constructed from key genes with accurate predictive utility
Although we have identified key genes that accurately predict the fatty acid metabolic subtype and prognosis of patients, their roles lack effective validation. Therefore, we analyzed the risk scores of patients based on the multivariate Cox regression model to discriminate them into high risk and low risk groups. According to the risk scores’ median, we demonstrated this classification’s effectiveness in different cohorts. The expressions of the 3 key genes in different risk groups exactly correspond to their hazard ratios (Fig. 6A–B). After stratifying patients by age, gender, and tumor stage, we found a higher proportion of patients defined as high risk in individuals with older or advanced tumors (III/IV), while individuals with younger or early-stage tumors (I/II) were more inclined to be low risk. At the same time, gender did not yield consistent results across cohorts (Fig. 6C–D). In addition, survival analyses of different cohorts all indicated that patients in the high risk group had a worse prognosis than those in the low risk group (Fig. 6E–F).
Figure 6.: Risk prognostic model constructed from 3 key genes in the TCGA cohort and the GEO cohort. (A–B) Distribution of patient risk scores, survival times, survival status, and expression levels of key genes. (C–D) Distribution of high- and low-risk individuals among patients stratified by age, gender, and tumor stage. (E–F) The K-M curve of high-risk and low-risk groups (TCGA: P = .0011; GEO: P = .00045; Log Rank test). K-M = Kaplan–Meier, GEO = gene expression omnibus, TCGA = the cancer genome atlas.
Notably, the univariate COX regression results showed that the risk score exhibited a high hazard ratio in both independent cohorts [TCGA: READ: HR = 2.72 (1.43–5.16); GSE87211: HR = 2.72 (1.66–4.46)] and was statistically significant (P < .01) (Fig. 7A–B). More importantly, for the different cohorts, the AUC value of the risk score was 0.91 and 0.708, respectively. In contrast, the tumor stage, which is often used to predict patient prognosis in clinical practice, was only 0.63 and 0.584 (Fig. 7C–D).
Figure 7.: Assessment of the predictive efficacy of a risk prognostic model. (A–B) The forest plot of age, gender, stage and the risk score (**, P ≤ .01). (C–D) The ROC curve of age, gender, stage and the risk score. ROC = receiver operating characteristic.
4. Discussion
Metabolic abnormality is an important feature of cancer cells.[26] It actively promotes the occurrence and development of cancer cells by regulating gene expression, resisting apoptosis, and promoting cell proliferation and migration.[26,27] In-depth studies of the metabolic profile of cancer cells have shown that nutrients such as acetic acid, fatty acid, lactic acid, branched-chain amino acid, serine, and glycine directly participate in cell transformation or support cancer growth through metabolism.[28,29] Not all studies suggest that abnormal intracellular metabolism promotes cancer development. Whether certain substances thought to be key to cancer metabolism must promote the proliferation and migration of cancer cells remains controversial.[30,31] However, it is undeniable that metabolism-based research provides reliable and effective biomarkers for cancer diagnosis and prognosis, and metabolism-based targeted cancer therapy has broad application prospects.
This study found differential expression of fatty acid metabolism genes in READ patients. Based on the expression pattern, patients could be divided into 2 stable subpopulations. Pathway enrichment and metabolic signature identification results for subpopulations showed that the 2 patient subpopulations corresponded to 2 completely different fatty acid metabolic subtypes: fatty acid catabolic subtype and fatty acid anabolic subtype. Subsequent survival analyses indicated that the patient subpopulation corresponding to the fatty acid catabolic subtype had a poorer prognosis. Notably, fatty acid beta-oxidation has been shown to promote cancer cells growth and proliferation in a variety of cancers, including READ.[27,29,32]
In addition to the unique metabolic profile exhibited, the local microenvironment of cancer cells is rich in immune cells and stromal cells, which build a complex ecosystem of cancers and together contribute to cancer development and progression.[27,29,33] In the study, we used the MCPcounter,[22] a validated and robust quantitative method, to assess the relative abundance of immune and stromal cells within the tumor microenvironment of READ patients. Subsequently, we compared the differences in the relative abundance of the same cells between individuals of different fatty acid metabolic subtypes, where the relative abundance of mDCs was significantly lower in patients with fatty acid catabolic subtype than others. Furthermore, the relative abundance of mDCs in both groups was significantly negatively correlated with the signature score of fatty acid catabolic pathway and positively correlated with the signature score of fatty acid anabolic pathway. As specialized antigen-presenting cells of the immune system, dendritic cells (DCs) have the unique ability to attract and activate naïve CD4 and CD8 cells.[34] Among them, DCs-induced cytotoxic T cell responses are an important pathway for individuals to fight against tumors.[34,35] Notably, by comparing the ability of different subtypes of DCs to initiate cytotoxic T cells, Giulia Nizzoli et al[36] found that mDCs produce higher levels of cytotoxic molecules (IL12) compared to plasmacytoid dendritic cells and are the main effector cells for inducing cytotoxic T cell responses. Furthermore, unlike conventional monocyte-derived DCs, mDCs targeting tumors for immunotherapy does not require an extensive culture period, is clinically reproducible, and has demonstrated encouraging immunological and clinical results in completed clinical trials.[37,38] Taken together, we found that READ patients exhibit fatty acid metabolic heterogeneity that significantly affects individual prognosis. We also suggest that targeted metabolic regimens for fatty acid catabolic subtypes with low mDCs infiltration may yield better clinical outcomes, while the potential efficacy of mDCs-mediated immunotherapy is relatively stronger for fatty acid synthetic subtypes.
Further analysis showed that the expression levels of AQP7, XIST, and IL4I1 could predict the prognostic risk of patients. The high expression of AQP7 and IL4I1 indicated a poorer prognosis, and the high expression of XIST indicated a better prognosis. As an aquaglyceroporin, AQP7 promotes the hydrolysis of intracellular triglycerides into glycerol and fatty acids, while enhancing the transport of extracellular glycerol, and intracellular fatty acids, which provide adenosine triphosphate through β-oxidation.[39,40] As an immunosuppressive enzyme, IL4I1 induces cancer escape by inhibiting T-cell proliferation,[41] and our study found that high expression of IL4I1 may alter fatty acid metabolism patterns in READ patients and increase individual prognostic risk. Similarly, IL4I1 was found in the prognostic model of renal clear cell carcinoma constructed by Zhu et al.[42] XIST is a long noncoding RNA that mediates female mammalian development by recruiting protein complexes and is considered a female characteristic gene.[43] Multiple studies[44,45] suggested that women appear to be less likely to develop READ than men, and men are more likely to die from READ than women, which is consistent with the finding that XIST protects READ patients in this study. Subsequently, we established a risk prognostic model based on AQP7, XIST, and IL4I1 and calculated the risk score for each patient. Based on the risk score, we distinguished patients into a high risk group and a low risk group, in which the survival rate of patients in the high risk group decreased significantly with a longer follow up time and was much lower than that of the low risk group. In addition, we found that older and advanced tumor patients were more inclined to obtain a higher risk score, thar is, a poorer prognosis, which is consistent with the effect of age and tumor stage on patient prognosis in clinical practice,[2,3] indicating the validity of the model. More importantly, the univariate COX results suggest that the risk score has a good independent prognostic value, while the AUC value of the risk score is much higher than that of age, gender, and tumor stage, indicating that it is a more accurate predictor of patient prognosis.
This study has certain limitations. First, the fatty acid metabolism genes defined in the study may not be completely accurate due to data incompleteness, data update delay, and other issues in the molecular signatures database itself. Second, the source of data used in the risk prognostic model is limited to the TCGA and GEO databases, and the utility and credibility of the model will further improve if it is applied for prospective clinical trial cohorts.
5. Conclusion
In conclusion, this study revealed the heterogeneity of individual fatty acid metabolism and effectively defined 2 mutually independent fatty acid metabolic subtypes through systematic analysis of fatty acid metabolism genes in READ patients. In addition, the study confirmed that READ patients with fatty acid catabolic subtypes have a poorer prognosis and lower infiltration of mDCs in the tumor microenvironment, which is a good guide for the selection of therapeutic measures. More importantly, we constructed a risk prognostic model based on metabolic subtypes and individual survival status and demonstrated its desirable predictive utility. We expect that this study will help clinicians achieve early assessment and long term monitoring of the prognosis of READ patients and provide a reference for individualized patient treatment.
Acknowledgments
We acknowledge TCGA and GEO database for providing their platforms and contributors for uploading their meaningful datasets.
Author contributions
Conceptualization: Jian Wang, Dong Shang.
Data curation: Jian Wang, Yi Dong.
Formal analysis: Jian Wang, Yi Dong.
Investigation: Jian Wang, Yi Dong.
Methodology: Jian Wang, Yi Dong.
Project administration: Jian Wang, Yi Dong.
Resources: Jian Wang, Yi Dong.
Software: Jian Wang, Yi Dong.
Supervision: Jian Wang.
Validation: Jian Wang, Yi Dong.
Visualization: Jian Wang.
Writing – original draft: Jian Wang.
Writing – review & editing: Jian Wang, Dong Shang.
References
[1]. Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49.
[2]. Dekker E, Tanis PJ, Vleugels JLA, et al. Colorectal cancer. Lancet. 2019;394:1467–80.
[3]. Keum N, Giovannucci E. Global burden of colorectal cancer: emerging trends, risk factors and prevention strategies. Nat Rev Gastroenterol Hepatol. 2019;16:713–32.
[4]. Loree JM, Pereira AAL, Lam M, et al. Classifying colorectal cancer by tumor location rather than sidedness highlights a continuum in mutation profiles and consensus molecular subtypes. Clin Cancer Res. 2018;24:1062–72.
[5]. Deng Y. Rectal cancer in Asian vs. Western countries: why the variation in incidence? Curr Treat Options Oncol. 2017;18:64.
[6]. Devesa SS, Chow WH. Variation in colorectal cancer incidence in the United States by subsite of origin. Cancer. 1993;71:3819–26.
[7]. Song M, Garrett WS, Chan AT. Nutrients, foods, and colorectal cancer prevention. Gastroenterology. 2015;148:1244–60.e16.
[8]. Li J, Ma X, Chakravarti D, et al. Genetic and biological hallmarks of colorectal cancer. Genes Dev. 2021;35:787–820.
[9]. Villéger R, Lopès A, Veziant J, et al. Microbial markers in colorectal cancer detection and/or prognosis. World J Gastroenterol. 2018;24:2327–47.
[10]. Liu Z, Liu Z, Zhou X, et al. A glycolysis-related two-gene risk model that can effectively predict the prognosis of patients with rectal cancer. Hum Genomics. 2022;16:5.
[11]. The Cancer Genome Atlas Network. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012;487:330–7.
[12]. Zuo D, Li C, Liu T, et al. Construction and validation of a metabolic risk model predicting prognosis of colon cancer. Sci Rep. 2021;11:6837.
[13]. Crotti S, Agnoletto E, Cancemi G, et al. Altered plasma levels of decanoic acid in colorectal cancer as a new diagnostic biomarker. Anal Bioanal Chem. 2016;408:6321–8.
[14]. Zhu Y, Gu L, Lin X, et al. Dynamic regulation of ME1 phosphorylation and acetylation affects lipid metabolism and colorectal tumorigenesis. Mol Cell. 2020;77:138–149.e5.
[15]. Hu Y, Gaedcke J, Emons G, et al. Colorectal cancer susceptibility loci as predictive markers of rectal cancer prognosis after surgery. Genes Chromosomes Cancer. 2018;57:140–9.
[16]. Liberzon A, Subramanian A, Pinchback R, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739–40.
[17]. Wilkerson MD, Hayes DN. Consensus cluster plus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3.
[18]. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.
[19]. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
[20]. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47–e47.
[21]. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2:100141.
[22]. Zeng D, Ye Z, Shen R, et al. IOBR: multi-omics immuno-oncology biological research to decode tumor microenvironment and signatures. Front Immunol. 2021;12:687975.
[23]. Hänzelmann S, Castelo R, GSVA GJ. gene set variation analysis for microarray and RNA-Seq data. BMC Bioinf. 2013;14:7.
[24]. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17:218.
[25]. Yu G. Using ggtree to visualize data on tree-like structures. Curr Protoc Bioinformatics. 2020;69:e96.
[26]. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144:646–74.
[27]. DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv. 2016;2:e1600200.
[28]. White E. Exploiting the bad eating habits of Ras-driven cancers. Genes Dev. 2013;27:2065–71.
[29]. Vander Heiden MG, DeBerardinis RJ. Understanding the intersections between metabolism and cancer biology. Cell. 2017;168:657–69.
[30]. Ying H, Kimmelman AC, Lyssiotis CA, et al. Oncogenic Kras maintains pancreatic tumors through regulation of anabolic glucose metabolism. Cell. 2012;149:656–70.
[31]. Sullivan LB, Gui DY, Hosios AM, et al. Supporting aspartate biosynthesis is an essential function of respiration in proliferating cells. Cell. 2015;162:552–63.
[32]. Monaco ME. Fatty acid metabolism in breast cancer subtypes. Oncotarget. 2017;8:29487–500.
[33]. Xiao Z, Dai Z, Locasale JW. Metabolic landscape of the tumor microenvironment at single cell resolution. Nat Commun. 2019;10:3763.
[34]. Steinman RM. The dendritic cell system and its role in immunogenicity. Annu Rev Immunol. 1991;9:271–96.
[35]. Tel J, Schreibelt G, Sittig SP, et al. Human plasmacytoid dendritic cells efficiently cross-present exogenous Ags to CD8+ T cells despite lower Ag uptake than myeloid dendritic cell subsets. Blood. 2013;121:459–67.
[36]. Nizzoli G, Krietsch J, Weick A, et al. Human CD1c+ dendritic cells secrete high levels of IL-12 and potently prime cytotoxic T-cell responses. Blood. 2013;122:932–42.
[37]. Bol KF, Schreibelt G, Rabold K, et al. The clinical application of cancer immunotherapy based on naturally circulating dendritic cells. J ImmunoTher Cancer. 2019;7:109.
[38]. Tacken PJ, de Vries IJM, Torensma R, et al. Dendritic-cell immunotherapy: from ex vivo loading to in vivo targeting. Nat Rev Immunol. 2007;7:790–802.
[39]. Lebeck J. Metabolic impact of the glycerol channels AQP7 and AQP9 in adipose tissue and liver. J Mol Endocrinol. 2014;52:R165–78.
[40]. Lebeck J, Søndergaard E, Nielsen S. Increased AQP7 abundance in skeletal muscle from obese men with type 2 diabetes. Am J Physiol Endocrinol Metab. 2018;315:E367–73.
[41]. Lasoudris F, Cousin C, Prevost-Blondel A, et al. IL4I1: an inhibitor of the CD8+ antitumor T-cell response in vivo. Eur J Immunol. 2011;41:1629–38.
[42]. Liu M, Pan Q, Xiao R, et al. A cluster of metabolism-related genes predict prognosis and progression of clear cell renal cell carcinoma. Sci Rep. 2020;10:12949.
[43]. Brockdorff N, Bowness JS, Wei G. Progress toward understanding chromosome silencing by Xist RNA. Genes Dev. 2020;34:733–44.
[44]. Patel SG, Karlitz JJ, Yen T, et al. The rising tide of early-onset colorectal cancer: a comprehensive review of epidemiology, clinical features, biology, risk factors, prevention, and early detection. Lancet Gastroenterol Hepatol. 2022;7:262–74.
[45]. Motsuku L, Chen WC, Muchengeti MM, et al. Colorectal cancer incidence and mortality trends by sex and population group in South Africa: 2002-2014. BMC Cancer. 2021;21:129.
