- Open Access
A lipid metabolism-related genes prognosis biomarker associated with the tumor immune microenvironment in colorectal carcinoma
BMC Cancer volume 21, Article number: 1182 (2021)
Background and aim
Lipid metabolic reprogramming is considered to be a new hallmark of malignant tumors. The purpose of this study was to explore the expression profiles of lipid metabolism-related genes (LMRG) in colorectal cancer (CRC).
The lipid metabolism statuses of 500 CRC patients from the Cancer Genome Atlas (TCGA) and 523 from the Gene Expression Omnibus (GEO GSE39582) database were analyzed. The risk signature was constructed by univariate Cox regression and least absolute shrinkage and selection operator (LASSO) Cox regression.
A novel four-LMRG signature (PROCA1, CCKBR, CPT2, and FDFT1) was constructed to predict clinical outcomes in CRC patients. The risk signature was shown to be an independent prognostic factor for CRC and was associated with tumour malignancy. Principal components analysis demonstrated that the risk signature could distinguish between low- and high-risk patients. There were significantly differences in abundances of tumor-infiltrating immune cells and mutational landscape between the two risk groups. Patients in the low-risk group were more likely to have higher tumor mutational burden, stem cell characteristics, and higher PD-L1 expression levels. Furthermore, a genomic-clinicopathologic nomogram was established and shown to be a more effective risk stratification tool than any clinical parameter alone.
This study demonstrated the prognostic value of LMRG and showed that they may be partially involved in the suppressive immune microenvironment formation.
Colorectal cancer (CRC) has the third highest incidence among malignancies worldwide and is the second most common cause of cancer-related deaths . As the most common gastrointestinal malignancy, CRC is associated with well-known risk factors, including poor dietary patterns, obesity, alcoholic abuse, smoking, and physical inactivity. Besides environmental influences, divergent genomic factors also contribute to the complexity of CRC . Despite clinical application of cancer screening procedures and effective treatments, the morbidity and mortality of CRC remain consistently elevated in most countries, especially in economically underdeveloped regions . Conventional therapeutic approaches, such as surgical resection combined with chemotherapy, have been shown to improve survival and quality of life; however, the overall survival (OS) of CRC patients remains unsatisfactory . Although the 5-year survival rate of patients with early CRC (stage I and II disease) is more than 60%, it is less than 10% for stage IV patients and some stage III patients with metachronous distant metastasis [3, 5]. Hence, a reliable risk prediction model is essential to improve the clinical prognosis of CRC by providing targeted screening and individualized intervention measures . At present, the tumor node metastasis (TNM) staging is still the most widely used method for prognostic evaluation in CRC . Other predictive factors include pathological type, differentiation degree and microvascular/serosa invasion, etc. However, these clinicopathological biomarkers cannot provide precise prognostic guidelines as they do not fully capture disease information. With advances in next-generation sequencing, inter- and intra-tumor molecular heterogeneity has been highlighted, increasing the difficulty of risk stratification . Therefore, molecular characteristics of tumors need to be included in new prognostic risk models.
Accumulating evidence indicates that tumor metabolic heterogeneity is related to clinical outcome, epigenetics status, and treatment resistance [9,10,11]. Oncogene-driven metabolic reprogramming improves cellular fitness to meet bioenergetic, biosynthetic, and redox balance demands, thereby providing a selective advantage during tumorigenesis . Metabolic reprogramming is considered to be a new hallmark of malignant tumors . Although most studies of alterations in cancer metabolism have focused on glucose metabolism (i.e., the Warburg effect), the role of abnormal lipid metabolism in cancer cells has been gradually recognized over the past few years [14, 15]. A rapidly proliferating cancer cell requires more energy than a normal cell and meets its biological needs by activating an endogenous production pathway or increasing its intake . ATP generated by fatty acid oxidation is an important energy source for cancer cells when energy provision is insufficient. Adipocytes and free fatty acids in the hypoxic tumor microenvironment are markedly conducive to cancer proliferation, progression, invasion, and metastasis [17, 18]. It has been reported that cancer cells rely mostly on endogenous adipogenesis rather than uptake of exogenous fatty acids, which is more common in normal cells . Hence, abnormal lipid metabolism, especially fatty acid synthesis and oxidation, is increasingly regarded as an important feature of metabolic reprogramming.
Epidemiological studies have shown that serum triglyceride levels are associated with susceptibility to CRC [20, 21]. Wang et al. confirmed that alterations in the abundances of individual lipids were present in CRC using shotgun lipidomics . The study showed increased expression of lipogenic enzymes involved in de novo adipogenesis (fatty acid synthesis) in CRC, including fatty-acid synthase, acetyl-CoA carboxylase, and carnitine palmitoyltransferase, whereas enzymes involved in fatty acid oxidation showed decreased expression. As a mitochondrial serine/threonine phosphatase, PGAM5 regulates a variety of metabolic pathways in vivo. Research by Zhu et al. showed that blocking PGAM5 would reduce lipid metabolism and inhibit the occurrence of CRC in mice . Gong and his colleagues showed that metabolism reprogramming in tumor-associated fibroblasts significantly enhanced the invasion and metastasis of CRC . Besides, drug resistance to antiangiogenic therapy frequently arises during cancer treatment; the underlying molecular mechanism of this resistance might include lipid metabolism reprogramming . Iwamoto et al. observed that blocking carnitine palmitoyl transferase 1A (CPT1A), a key enzyme in lipid metabolism, could restore sensitivity to antiangiogenic therapy . Thus, they introduced a promising approach to overcoming drug resistance to cancer therapies by combining conventional therapy and targeted lipid metabolism.
Previous studies have confirmed the close relationship between altered lipid metabolism and CRC in tumorigenesis, progression and treatment. A signature of lipid metabolism-related genes (LMRG) was shown to have high prognostic value in papillary thyroid cancer and diffuse glioma [26, 27]. In-depth study of lipid metabolomics characteristics could lead to better understanding of the progression of tumors and provide potential metabolic targets for the development of new treatment methods . However, the prognostic value of LMRG in CRC has not been verified by large-sample studies. The present study aimed to develop a novel risk signature based on LMRG to provide additional information for use in risk assessment and clinical decision-making in CRC.
A flow chart of this research is presented in Fig. 1. Level 3 RNA sequencing data (RNA-seq) data, mutation data, and matched clinical information were obtained from the TCGA CRC cohort (Data Release 25.0 - July 22, 2020, https://portal.gdc.cancer.gov/repository). The data search and selection strategy were as follows. (1) The keywords for cases were “colon and rectum [Primary Site]”, “TCGA [Program]”, “TCGA-COAD, TCGA-READ [Project]”, “Adenomas and Adenocarcinomas [Disease Type]”. (2) The keywords for files were “Transcriptome Profiling [Data Category]”, “Gene Expression Quantification [Data Category]”, “RNA-Seq [Experimental Strategy]”, “HTSeq-FPKM [Workflow Type]”, “.txt Format [Data Format]” . The RNA-seq matrix file was annotated using the human General Transfer Format (hunman.gtf) from the Ensembl database (https://www.ensembl.org/) with the Strawberry Perl software (version 22.214.171.124, https://strawberryperl.com/). The RNA-seq transcriptome data (FPKM) from the TCGA cohort were converted to log2(TPM + 1) to obtain normalized counts.
Raw CEL and clinical data for CRC were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/); the search used the keywords “colon cancer”, “rectum cancer”, “adenocarcinoma”, “Homo sapiens”, and “Expression profiling.” The gene expression microarray dataset GSE39582 was selected and downloaded because it had the largest sample size for CRC . Probe IDs were matched to gene symbols using the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). The mean expression value of the probes was used as the expression value for the gene in question if multiple probes were mapped to a single gene. The robust multichip average (RMA) was used to normalize the raw data from GEO cohort by the R package affy .
Criteria for study exclusion were: (1) patients with unknown survival status, follow-up information, and disease stage; and (2) patients who died within a follow-up period of 30 days. Consequently, 544 cases (500 tumor and 44 normal samples) meeting the criteria were included in the training set, and 542 cases (523 tumor and 19 normal samples) were included in the validation set (Supplementary Table 1). The TCGA cohort was used to establish a risk signature, and the GEO set was used for validation. Both the TCGA and GEO databases are publicly available; therefore, the present study did not require approval from the local ethics committees.
Identification of LMRG
The five LMRG sets were collected from the Molecular Signature Database [27, 31] (Supplementary Table 2). A total of 1044 genes were found to be involved in the lipid metabolism process after removal of overlapping genes.
Construction of risk signature based on LMRG
The shared LMRG in the GEO and TCGA sets were selected for subsequent analysis. Differentially expressed genes (DEGs) between tumor and normal samples were screened in the TCGA cohort using the R package limma and sva with a false discovery rate (FDR) < 0.05 [32, 33]. Gene ontology (GO) and KEGG analyses of DEGs were performed using R package clusterProfiler, org. Hs.eg.db, enrichplot and ggplot2 [34, 35]. The LMRG with prognostic value was identified using univariate Cox analysis. Next, the overlapping genes between DEGs and prognostic-related genes were identified using a Venn diagram for subsequent analysis (http://bioinformatics.psb.ugent.be/webtools/Venn/). Least absolute shrinkage and selection operator (LASSO) Cox regression was used to select the best prediction model based on these mutual genes in TCGA CRC patients . LASSO analysis was performed using R package glmnet and survival, and the optimal value of the penalty parameter was determined by 10-fold cross-validation . Risk signatures were generated from the TCGA and GEO cohorts based on the expression levels of LMRG and the corresponding coefficients. The risk score of each sample was calculated by the following formula: (expGene: the expression level of LMRG in the TCGA or GEO cohort; Coef: the coefficient of LMRG in the LASSO Cox regression model in the training set).
Prognostic value of the risk signature in training and validation group
The patients were stratified into high- and low-risk groups according to the median value of the risk score. Kaplan Meier (K-M) survival curves with the Log-rank test were constructed to demonstrate the prognostic ability of the risk signature. In addition, the area under the curve values (AUCs) of the receiver operating characteristic (ROC) curves for 1-, 3-, and 5-year survival were calculated using R package survivalROC to evaluate the performance of those two signatures .
Gene set enrichment analysis (GSEA)
To explore potential molecular mechanism between the two groups, GSEA was carried out between the high- and low-risk groups. The annotated gene set list, h.all.v7.2.symbols.gmt (Hallmarks), was selected as the reference gene set from the Molecular Signature Database .
Independence of the risk signature from other clinicopathological parameters
To determine the independence of the risk signature from other clinical parameters, univariate and multivariate analyses of the risk score with age, gender, and tumor stage were performed. Forest plots were used to show the independence of the risk score.
Correlation between the risk signature and other clinicopathological parameters
The associations between the risk signature and clinicopathological parameters in TCGA -CRC cohort, including age, gender, tumor stage, pathological T stage, N stage, and M stage, were further explored. Patients were stratified into subgroups of age ≤ 65 years and age ≥ 65 years, female and male, pathological tumor stage I + II and stage III + IV, T1 + T2 and T3 + 4, N0 and N1 + 2, M0 and M1. K-M survival analysis of the aforementioned paired subgroups was performed.
Cancer stem cells are highly dependent on lipid metabolism to maintain their stem cell characteristics. One study showed that cancer stem cells in CRC have higher lipid metabolism levels than tumor cells or normal colonic epithelial cells . Malta et al. had developed a novel analysis tool to assess stemness features based on gene expression . In this study, the mRNA expression-based stemness index (mRNAsi) for each CRC patient was obtained from the results of Malta et al. Besides, CD133 is a marker gene of many tumor stem cells. Hence, the relationships of the risk score with the mRNAsi score and CD133 mRNA were also analyzed.
Nomogram construction and validation
A nomogram integrating the risk signature and other clinicopathological factors was established for prognostic evaluation using the R package rms (https://CRAN.R-project.org/package=rms). The AUCs of the ROC curve were demonstrated to assess the predictive capability of the nomogram. Calibration curves were also established to evaluate the predictive accuracy of the nomogram.
Estimation of relative abundance of immune cell types in different risk groups
The CIBERSORT algorithm (https://cibersort.stanford.edu), which quantifies the relative abundance of immune cells based on specific gene expression profiles, is used to assess the distribution of 22 immune cell types in CRC samples from TCGA cohort . The P-value, correlation coefficient, and root mean squared error were also calculated to evaluate the accuracy of the results for each patient. Samples with P < 0.05 were used to compare immune cell abundance in different risk groups.
The top 20 genes with the highest mutation frequency in TCGA -CRC cohort were analyzed in both high- and low-risk groups using the R package GenVisR . Visualization based on the somatic mutation data in Mutation Annotation was performed using the R package maftoools . Tumor mutational burden (TMB) is known to be associated with the efficacy of adjuvant chemotherapy (fluoropyrimidine plus oxaliplatin regimen) in CRC [44, 45]. TMB is also an independent predictor of immunotherapy response in a variety of tumors [38, 39]. Therefore, the relationship between risk score and TMB was also explored in this study. The TMB score was generated as the total number of somatic mutations divided by the number of exons in each sample . The exon size is often approximately estimated at 38 megabases.
All the statistical analyses and drawings in this study used R (version 3.6.3) or GraphPad Prism (version 8.3.0). Differences between continuous variable were analyzed using t-test. Fisher’s exact test or chi-square test was employed for comparisons of categorical variables. Log-rank test was used to estimate the differences among K-M survival curves. P < 0.05 (two-tailed) was considered significant.
Identification of differentially expressed and prognosis-related genes in LMRG
Expression data of LMRG were extracted from the TCGA and GEO cohorts. A total of 945 shared LMRG were matched. There were 729 DEGs between normal and tumor tissues, including 365 upregulated and 364 downregulated genes, in the training cohort when the cut-off was set to FDR < 0.05 (Fig. 2A). Univariate Cox analysis was employed to filter genes with prognostic value from the 945 intersecting genes. Finally, 57 LMRG were shown to be related to prognosis in the training set. A total of 47 shared LMRG were found to be both DEGs and prognostic genes according to the Venn diagram (Fig. 2B, C).
Functional enrichment analysis identified terms related to lipid metabolism. The top five biological process terms were glycerolipid metabolic process, glycerophospholipid metabolic process, lipid modification, phospholipid metabolic process, and regulation of lipid metabolic process (Fig. 2D). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the 47 shared LMRG were mainly involved in glycerophospholipid metabolism, phosphatidylinositol signaling system, PPAR signaling pathway, fatty acid degradation, and fatty acid metabolism (Fig. 2E).
Construction and validation of gene signature
Next, the expression profiles of the 47 LMRG were used to establish a risk signature using LASSO Cox regression analysis (Fig. 3A, B). Finally, four LMRG, namely PROCA1, CCKBR, CPT2, and FDFT1, were used to establish the optimal lipid metabolism-related risk signature. The risk score for each patient was calculated by the following formula: risk score = (PROCA1*0.03071) + (CCKBR*0.58956) + (CPT2*0.00972) + (FDFT1*0.01381).
The predictive value of this risk signature was evaluated using ROC curves. The AUCs of the signature were 0.6901 at 1 year, 0.6776 at 3 years, and 0.5945 at 5 years in the training set (Fig. 3C). The patients were dichotomized into two risk group according to the median risk score. K-M survival curves showed that high-risk patients had significantly shorter OS than low-risk cases in the training set (Fig. 3D, P < 0.001). As shown in Fig. 3E, more patients died in the high-risk group, whereas the majority survived in the low-risk group. In the training set, principal components analysis (PCA) was used to obtain expression patterns in the low- and high-risk groups. There was no significant difference in risk status between the two groups when PCA was performed with all genes (Fig. 3F). However, the risk groups could be distinguished using the risk signature (Fig. 3G).
In the validation set, the high-risk patients again had significantly shorter OS than those in the low-risk group (Fig. 3H, P < 0.05). The AUCs were 0.6152 at 1 year, 0.6332 at 3 years, and 0.6358 at 5 years (Fig. 3I). The risk score distribution curve and survival status are shown in Fig. 3J.
Independent prognostic value of the risk signature
The results showed that the risk score based on the four-LMRG signature was an independent prognostic factor in the training set, with hazard ratio (HR) = 6.146 (95% confidence interval (CI) = 3.376–11.190 (P < 0.001; Fig. 4A) by univariate Cox analysis, and HR = 4.315 (95% CI = 2.321–8.022; P < 0.001; Fig. 4B) by multivariate Cox analysis. Similar results were obtained in the validation set, showing the independence of the risk signature with HR = 5.822 (95% CI = 2.567–13.205; P < 0.001; Fig. 4C) and HR = 5.395 (95% CI = 2.272–12.809; P < 0.001; Fig. 4D) by univariate and multivariate Cox analysis, respectively.
Relationships between risk score and clinicopathological features
Risk score was significantly associated with age (P < 0.0001; Fig. 5A) but not gender (P = 0.4168; Fig. 5B). In addition, the relationship between obesity and risk score was explored. Obesity was defined as a body mass index (BMI) equal to or greater than 30 . Interestingly, no significant difference in lipid metabolism-related risk score was found between the normal weight and obese groups (P = 0.4168; Fig. 5C).
The results showed that risk score was positively correlated with the degree of tumor progression. For instance, the risk scores in TNM stage III + IV patients were markedly higher than those in tumor stage I + II (P < 0.0001; Fig. 5D). Similar results were obtained in other subgroups, including pathological T stage (I + II vs. III + IV, P < 0.0001; Fig. 5E), N stage (N0 vs. N1 + 2, P < 0.0001; Fig. 5F), M stage (M0 vs. M1, P = 0.0024; Fig. 5G), and tumor location (left vs. right, P < 0.0001; Fig. 5H). There were also more stem cell characteristics in the low-risk group, including significantly higher levels of CD133 mRNA (P = 0.0479; Fig. 5I) and mRNAsi (P < 0.0001; Fig. 5J).
Subgroup analysis was performed to further assess whether the risk signature still had independent prognostic value within specific clinical parameters. The results showed that the risk signature retained its powerful prognostic prediction ability in subgroups based on age (age ≥ 65 or < 65, Fig. 5K and L), gender (male or female, Fig. 5M or 5N), tumor stage (III + IV, Fig. 5O), T stage (III + IV, Fig. 5P), N stage (N1 + 2 or N0, Fig. 5Q and R), M stage (M0, Fig. 5S). The OS of high-risk patients was significantly shorter than that of low-risk patients in the different clinicopathological characteristic subgroups.
Construction and validation of a nomogram combining clinicopathological features and risk signature
A nomogram was constructed based on several factors, namely age, gender, disease stage, and risk score, to provide a visualization tool for clinicians to predict the probability of 1-, 3- and 5-year OS in CRC patients. A total score could be calculated for each patient using this nomogram, where a higher total score indicates a worse outcome (Fig. 6A). This nomogram had high potential for clinical utility, with ROC AUC values of 0.7652 at 1 year, 0.8058 at 3 years, and 0.7972 at 5 years in the training set (Fig. 6B). Moreover, the calibration plots indicated that the actual observation probability was very close to the predictive probability of the nomogram (Fig. 6D-F).
The nomogram for the validation set is shown in Fig. 6D. The AUCs of the ROC curve were 0.7325 at 1 year, 0.7280 at 3 years, and 0.7173 at 5 years (Fig. 6C). The predictive probabilities obtained using the nomogram showed good consistency with the actual observed values (Fig. 6G-I). Hence, the prognostic nomogram model was considered to be robust.
Immune landscapes differ among two risk subgroups in TCGA -CRC cohort
The differences in abundances of immune cells between the high- and low-risk groups were further investigated. The proportions of 22 immune cells in each sample varied markedly (Fig. 7A and B). When the cut-off value was set to P < 0.05, 487 samples of the total set were included in the subsequent analysis. The results showed that monocytes (P = 0.011), M0 (P = 0.001), and M2-like macrophages (P < 0.001) were substantially enriched in the high-risk group, whereas patients in low-risk group had higher levels of plasma cells (P = 0.012), resting memory CD4 T cells (P = 0.008), and activated dendritic cells (P = 0.004) (Fig. 7C). Additionally, we also assessed the correlation between risk score and relative abundances of immune cells. As shown in the correlation heatmap, the results revealed that the risk score was negatively correlated with most of the immune cells (Supplementary Fig. 1A). PD-L1, as an immune checkpoint molecule, is often closely related to immunotherapy response [44, 48]. In this study, the risk score was negatively correlated with PD-L1 mRNA in the training set (R = − 0.13, P = 0.0046) (Fig. 7D).
Several gene signatures related to cancer invasion, including invasiveness signature, epithelial-mesenchymal transition (EMT), pan-fibroblast TGFb response signature (Pan-F-TBRS) and two key immune cells in tumor microenvironment (effector CD8+ T and cancer associated fibroblasts (CAFs)) were chosen to further investigation of the role of LMRG score in tumorigenesis and metastasis [49, 50]. The enrichment score of each gene signature was calculated using sample gene set enrichment analysis (ssGSEA) algorithm. Notably, compared with the low-risk group, the high-risk group had a higher enrichment score of EMT- and CAFs-related signatures, while lower enrichment score of CD8+ T effector, which was a critical component for the antitumor immune response (Fig. 7E; Supplementary Fig. 1B-H).
Mutation landscape of the risk signature in TCGA -CRC cohort
The top 20 genes with the highest mutation frequency in CRC are shown in the waterfall plots in Fig. 8A and B. The high- and low-risk groups had different genetic mutation landscapes. Genes with crucial biological functions in tumorigenesis, including TP53, PIK3CA, and MUC16, showed significant differences in mutation frequency between the two groups (Fig. 8C). Besides, the top 20 most common mutant genes in CRC, including FAT4, FUT9, LRP1B and ZFHX4, showed significant differences in risk score between mutated- and wild-type group (Supplementary Fig. 2A-D).
A total of 429 samples of TMB data were obtained for this analysis. The results of the analysis showed that the TMB score was higher in the low-risk group than in the high-risk group (P = 0.0002, Fig. 8D). We also evaluated the correlation between the risk score and TMB, and found a remarkable correlation between the two variables (r = − 0.19, P = 5.7e-05, Supplementary Fig. 2E). This could be one of the reasons for the better prognosis of patients in the low-risk group.
Distinct biological function pathways characterize high- and low-risk CRC patients in TCGA -CRC cohort
GSEA was also performed to explore whether relevant signaling pathways differed between the two risk groups. The results showed that the high-risk group was associated with Hedgehog signaling, KRAS signaling, Wnt/β catenin signaling, apical junction, epithelial-mesenchymal transition, and angiogenesis. Functional enrichment in the low-risk group was focused on energy-metabolism-related functions, including fatty acid metabolism and oxidative phosphorylation (Fig. 9).
In this study, the lipid metabolism statuses of 1023 CRC samples were analyzed using gene expression data from the public databases, and a risk signature (comprising PROCA1, CCKBR, CPT2, and FDFT1) was constructed. This is the first reported CRC prognostic risk signature based on LMRG. Patients with high risk scores had significantly shorter OS than those with low risk scores, in both training and validation sets. The risk signature could well distinguish high-risk CRC patients. The risk signature, as an independent prognostic indicator, might be a beneficial supplement to TNM staging to provide more accurate prognostic information. A nomogram combining the risk signature and clinical characteristics showed better capacity for risk stratification than any clinical parameter alone. Consistent results were obtained in an external cohort from the GEO database (dataset GSE39582). In general, the prediction model based on the four-LMRG risk signature has potential to be an effective and robust prognostic indicator in CRC.
CRC remains the leading cause of cancer-related deaths worldwide, yet it is one of the most preventable cancers. Appropriate screening strategies can reduce the incidence and mortality of CRC by early detection and elimination of precancerous disease. Current CRC screening guidelines are based mostly on two main risk factors: age and family history . However, the fact is that more than 80% of CRC cases occur in individuals with no positive family history in primary relatives . Obviously, the use of these criteria does not fully take into account the heterogeneity of CRC. Lipid metabolic reprogramming frequently emerges in intestinal tumor cells and has been reported to have crucial biological roles in cell proliferation, energy homeostasis, and signal-transduction [52, 53]. Metabolic status is associated with various clinical outcomes in tumor patients, and a metabolism-related risk signature could serve as a prognostic predictor in cancers . Thus, a predictive model that integrates metabolism-related genomic data and clinicopathological characteristics provides hope for primary and secondary prevention of CRC. The value of this signature is to identify high-risk individuals with CRC so that targeted treatment and more stringent postoperative follow-up can be adopted. It could also be used to distinguish low-risk patients to prevent excessive treatment and unnecessary screening.
Another important finding of this work was that lipid metabolism status was significantly correlated with immune infiltration levels in CRC. According to our investigation, the risk signature based on lipid metabolism was related to abundance of monocyte infiltration, especially M2-like macrophages. The high-risk group had higher levels of M2-like macrophages, suggesting that high-risk patients might have a suppressive tumor immune microenvironment (TIME). In this study, patients with low risk score usually have high-level immune cell infiltration, especially CD8+ T effector, exhibiting a “hot” tumor phenotype. In concordance with previous studies, inflammatory phenotype always has better prognosis [55, 56]. Meanwhile, it is found that significant activation of EMT, TGF-β pathway and CAFs in high-risk group, which are hallmarks of stromal activation, showing an immune-excluded phenotype. In immune-excluded tumors, immune cells, especially CD8+ T effector cells, reside in the stroma surrounding tumor cell making direct contact with CAFs rather than tumor cells, resulting in restricting anti-tumor immune response [57, 58]. Furthermore, CAFs-derived exosomes can strikingly shape cancer-promoting microenvironment through involved in multifarious metabolicl processes . It may partly explain significant difference in the clinical outcomes among different metabolic subgroups.
Tumor-associated macrophages (TAMs) are the most important immune cell component of the TIME. Metabolic reprogramming of TAMs shapes their functional subtype . TAMs are often characterized by M2-like macrophages and have a variety of tumor-promoting effects on the tumor microenvironment. Studies have found different metabolic patterns between pro-inflammatory and anti-inflammatory macrophages [61, 62]. Activated pro-inflammatory macrophages often depend on the glycolytic pathway for energy, whereas immunosuppressed macrophages are more inclined to use fatty acid oxidation . Wu et al. found that fatty acids in the TIME, especially unsaturated fatty acids, might promote the polarization of monocytes to M2-like macrophages with a strong immunosuppressive phenotype . Therefore, the lipid-related metabolism risk signature represents alterations in the TIME of CRC. In addition, the risk signature was associated with the immune checkpoint marker PD-L1, suggesting that it has potential as a metabolic marker for immunotherapy in CRC.
Some metabolic regulators have been considered as oncogenes or tumor suppressors. Carnitine palmitoyl transferase II (CPT2) is a rate-limiting enzyme for mitochondrial fatty acid transportation, with a critical role in regulating fatty acid oxidation. Gastric cancer and CRC patients with lower CPT2 expression level have better disease control rates than those with higher CPT2 expression . Blocking fatty acid oxidation by knocking out CPT1A/CPT2 via CRISPR-mediated has been shown to inhibit the invasive phenotype of radiotherapy-resistant breast cancer, suggesting that CPT2 is a potential metabolic target for breast cancer radiotherapy . Besides, CPT2 levels in tumor tissue were found to be associated with oxaliplatin-based chemotherapy sensitivity. Fatty acid catabolism can be inhibited by knocking out CPT2 or using the CPT2 inhibitor perhexiline. This inhibitory effect promotes induction of tumor cell apoptosis by oxaliplatin or other classic chemotherapy drugs . These results shows that CPT2 has the potential to be a metabolic therapeutic target for CRC.
The results of this study revealed that LMRG have prognostic value in CRC and may be partially involved in the suppressive immune microenvironment formation. LMRG, especially PROCA1, CCKBR, CPT2, and FDFT1, are potential prognostic markers and therapeutic targets for CRC. This lipid metabolism-related prognostic signature represents a potential biomarker for predicting the efficacy of chemotherapy and anti-PD-L1 therapy in CRC.
Availability of data and materials
All data can be obtained from TCGA CRC cohort (https://portal.gdc.cancer.gov/repository) and GEO GSE39583 cohort (https://www.ncbi.nlm.nih.gov/geo/).
the Cancer Genome Atlas
the Gene Expression Omnibus
lipid metabolism-related genes
the principal components analysis
least absolute shrinkage and selection operator
fragments per kilobase million
transcripts per million
receiver operating characteristic
area under the curves
tumor mutation burden
the mRNA expression-based stemness index
differentially expressed genes
false discovery rate
differentially expressed genes
tumor immune microenvironment
Keum N, Giovannucci E. Global burden of colorectal cancer: emerging trends, risk factors and prevention strategies. Nat Rev Gastroenterol Hepatol. 2019;16(12):713–32. https://doi.org/10.1038/s41575-019-0189-8.
Yang Y, Han Z, Li X, Huang A, Shi J, Gu J. Epidemiology and risk factors of colorectal cancer in China. Chin J Cancer Res. 2020;32(6):729–41. https://doi.org/10.21147/j.issn.1000-9604.2020.06.06.
Dekker E, Tanis PJ, Vleugels J, Kasi PM, Wallace MB. Colorectal cancer. Lancet. 2019;394(10207):1467–80. https://doi.org/10.1016/S0140-6736(19)32319-0.
Provenzale D, Ness RM, Llor X, Weiss JM, Abbadessa B, Cooper G, et al. NCCN guidelines insights: colorectal Cancer screening, version 2.2020. J Natl Compr Cancer Netw. 2020;18(10):1312–20. https://doi.org/10.6004/jnccn.2020.0048.
McQuade RM, Stojanovska V, Bornstein JC, Nurgali K. Colorectal Cancer chemotherapy: the evolution of treatment and new approaches. Curr Med Chem. 2017;24(15):1537–57. https://doi.org/10.2174/0929867324666170111152436.
Thomas M, Sakoda LC, Hoffmeister M, Rosenthal EA, Lee JK, van Duijnhoven F, et al. Genome-wide modeling of polygenic risk score in colorectal Cancer risk. Am J Hum Genet. 2020;107(3):432–44. https://doi.org/10.1016/j.ajhg.2020.07.006.
Weiser MR. AJCC 8th edition: colorectal Cancer. Ann Surg Oncol. 2018;25(6):1454–5. https://doi.org/10.1245/s10434-018-6462-1.
Zhang Y, Song J, Zhao Z, Yang M, Chen M, Liu C, et al. Single-cell transcriptome analysis reveals tumor immune microenvironment heterogenicity and granulocytes enrichment in colorectal cancer liver metastases. Cancer Lett. 2020;470:84–94. https://doi.org/10.1016/j.canlet.2019.10.016.
Karasinska JM, Topham JT, Kalloger SE, Jang GH, Denroche RE, Culibrk L, et al. Altered gene expression along the glycolysis-cholesterol synthesis Axis is associated with outcome in pancreatic Cancer. Clin Cancer Res. 2020;26(1):135–46. https://doi.org/10.1158/1078-0432.CCR-19-1543.
Cristea S, Coles GL, Hornburg D, Gershkovitz M, Arand J, Cao S, et al. The MEK5-ERK5 kinase Axis controls lipid metabolism in small-cell lung Cancer. Cancer Res. 2020;80(6):1293–303. https://doi.org/10.1158/0008-5472.CAN-19-1027.
Shen C, Xuan B, Yan T, Ma Y, Xu P, Tian X, et al. m(6)A-dependent glycolysis enhances colorectal cancer progression. Mol Cancer. 2020;19(1):72.
Vander HM, DeBerardinis RJ. Understanding the intersections between metabolism and Cancer biology. Cell. 2017;168(4):657–69. https://doi.org/10.1016/j.cell.2016.12.039.
Thakur C, Chen F. Connections between metabolism and epigenetics in cancers. Semin Cancer Biol. 2019;57:52–8. https://doi.org/10.1016/j.semcancer.2019.06.006.
Luo X, Cheng C, Tan Z, Li N, Tang M, Yang L, et al. Emerging roles of lipid metabolism in cancer metastasis. Mol Cancer. 2017;16(1):76. https://doi.org/10.1186/s12943-017-0646-3.
Huang C, Freter C. Lipid metabolism, apoptosis and cancer therapy. Int J Mol Sci. 2015;16(1):924–49. https://doi.org/10.3390/ijms16010924.
Beloribi-Djefaflia S, Vasseur S, Guillaumond F. Lipid metabolic reprogramming in cancer cells. Oncogenesis. 2016;5(1):e189. https://doi.org/10.1038/oncsis.2015.49.
Lazar I, Clement E, Dauvillier S, Milhas D, Ducoux-Petit M, LeGonidec S, et al. Adipocyte exosomes promote melanoma aggressiveness through fatty acid oxidation: a novel mechanism linking obesity and Cancer. Cancer Res. 2016;76(14):4051–7. https://doi.org/10.1158/0008-5472.CAN-16-0651.
Iwamoto H, Abe M, Yang Y, Cui D, Seki T, Nakamura M, et al. Cancer lipid metabolism confers antiangiogenic drug resistance. Cell Metab. 2018;28(1):104–17. https://doi.org/10.1016/j.cmet.2018.05.005.
Merino SM, Gomez DCM, Moreno RJ, Falagan MS, Sanchez MR, Casado E, et al. Lipid metabolism and lung cancer. Crit Rev Oncol Hematol. 2017;112:31–40. https://doi.org/10.1016/j.critrevonc.2017.02.001.
Rodriguez-Broadbent H, Law PJ, Sud A, Palin K, Tuupanen S, Gylfe A, et al. Mendelian randomisation implicates hyperlipidaemia as a risk factor for colorectal cancer. Int J Cancer. 2017;140(12):2701–8. https://doi.org/10.1002/ijc.30709.
Liu T, Peng F, Yu J, Tan Z, Rao T, Chen Y, et al. LC-MS-based lipid profile in colorectal cancer patients: TAGs are the main disturbed lipid markers of colorectal cancer progression. Anal Bioanal Chem. 2019;411(20):5079–88. https://doi.org/10.1007/s00216-019-01872-5.
Wang Y, Hinz S, Uckermann O, Honscheid P, von Schonfels W, Burmeister G, et al. Shotgun lipidomics-based characterization of the landscape of lipid metabolism in colorectal cancer. Biochim Biophys Acta Mol Cell Biol Lipids. 1865;2020(3):158579. https://doi.org/10.1016/j.bbalip.2019.158579.
Zhu Y, Gu L, Lin X, Liu C, Lu B, Cui K, et al. Dynamic regulation of ME1 phosphorylation and acetylation affects lipid metabolism and colorectal tumorigenesis. Mol Cell. 2020;77(1):138–49. https://doi.org/10.1016/j.molcel.2019.10.015.
Gong J, Lin Y, Zhang H, Liu C, Cheng Z, Yang X, et al. Reprogramming of lipid metabolism in cancer-associated fibroblasts potentiates migration of colorectal cancer cells. Cell Death Dis. 2020;11(4):267.
Cao Y. Adipocyte and lipid metabolism in cancer drug resistance. J Clin Invest. 2019;129(8):3006–17. https://doi.org/10.1172/JCI127201.
Ma B, Jiang H, Wen D, Hu J, Han L, Liu W, et al. Transcriptome analyses identify a metabolic gene signature indicative of dedifferentiation of papillary thyroid Cancer. J Clin Endocrinol Metab. 2019;104(9):3713–25. https://doi.org/10.1210/jc.2018-02686.
Wu F, Zhao Z, Chai RC, Liu YQ, Li GZ, Jiang HY, et al. Prognostic power of a lipid metabolism gene panel for diffuse gliomas. J Cell Mol Med. 2019;23(11):7741–8. https://doi.org/10.1111/jcmm.14647.
Munir R, Lisec J, Swinnen JV, Zaidi N. Lipid metabolism in cancer cells under metabolic stress. Br J Cancer. 2019;120(12):1090–8. https://doi.org/10.1038/s41416-019-0451-4.
Marisa L, de Reynies A, Duval A, Selves J, Gaub MP, Vescovo L, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med. 2013;10(5):e1001453. https://doi.org/10.1371/journal.pmed.1001453.
Gautier L, Cope L, Bolstad BM, Irizarry RA. Affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20(3):307–15. https://doi.org/10.1093/bioinformatics/btg405.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50. https://doi.org/10.1073/pnas.0506580102.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3. https://doi.org/10.1093/bioinformatics/bts034.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51. https://doi.org/10.1093/nar/gkaa970.
Gui J, Li H. Penalized cox regression analysis in the high-dimensional and low-sample size settings, with applications to microarray gene expression data. Bioinformatics. 2005;21(13):3001–8. https://doi.org/10.1093/bioinformatics/bti422.
Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Epigenetics. 2019;11(1):123. https://doi.org/10.1186/s13148-019-0730-1.
Martinez-Camblor P, Pardo-Fernandez JC. Smooth time-dependent receiver operating characteristic curve estimators. Stat Methods Med Res. 2018;27(3):651–74. https://doi.org/10.1177/0962280217740786.
Tirinato L, Liberale C, Di Franco S, Candeloro P, Benfante A, La Rocca R, et al. Lipid droplets: a new player in colorectal cancer stem cells unveiled by spectroscopic imaging. Stem Cells. 2015;33(1):35–44. https://doi.org/10.1002/stem.1837.
Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies Stemness features associated with oncogenic dedifferentiation. Cell. 2018;173(2):338–54. https://doi.org/10.1016/j.cell.2018.03.034.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7. https://doi.org/10.1038/nmeth.3337.
Skidmore ZL, Wagner AH, Lesurf R, Campbell KM, Kunisaki J, Griffith OL, et al. GenVisR: genomic visualizations in R. Bioinformatics. 2016;32(19):3012–4. https://doi.org/10.1093/bioinformatics/btw325.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56. https://doi.org/10.1101/gr.239244.118.
Llosa NJ, Luber B, Siegel N, Awan AH, Oke T, Zhu Q, et al. Immunopathologic stratification of colorectal Cancer for checkpoint blockade immunotherapy. Cancer Immunol Res. 2019;7(10):1574–9. https://doi.org/10.1158/2326-6066.CIR-18-0927.
Lee DW, Han SW, Bae JM, Jang H, Han H, Kim H, et al. Tumor mutation burden and prognosis in patients with colorectal Cancer treated with adjuvant Fluoropyrimidine and Oxaliplatin. Clin Cancer Res. 2019;25(20):6141–7. https://doi.org/10.1158/1078-0432.CCR-19-1105.
Merino DM, LM MS, Fabrizio D, Funari V, Chen SJ, White JR, et al. Establishing guidelines to harmonize tumor mutational burden (TMB): in silico assessment of variation in TMB quantification across diagnostic platforms: phase I of the Friends of Cancer Research TMB Harmonization Project. J Immunother Cancer. 2020;8(1).
Flegal KM, Ogden CL, Fryar C, Afful J, Klein R, Huang DT. Comparisons of self-reported and measured height and weight, BMI, and obesity Prevalence from National Surveys: 1999-2016. Obesity (Silver Spring). 2019;27(10):1711–9. https://doi.org/10.1002/oby.22591.
Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–75. https://doi.org/10.1038/s41575-019-0126-x.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6) A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53.
Herrera M, Berral-Gonzalez A, Lopez-Cade I, Galindo-Pumarino C, Bueno-Fortes S, Martin-Merino M, et al. Cancer-associated fibroblast-derived gene signatures determine prognosis in colon cancer patients. Mol Cancer. 2021;20(1):73. https://doi.org/10.1186/s12943-021-01367-x.
Ladabaum U, Dominitz JA, Kahi C, Schoen RE. Strategies for colorectal Cancer screening. Gastroenterology. 2020;158(2):418–32. https://doi.org/10.1053/j.gastro.2019.06.043.
Aguilar E, Esteves P, Sancerni T, Lenoir V, Aparicio T, Bouillaud F, et al. UCP2 deficiency increases Colon tumorigenesis by promoting lipid synthesis and depleting NADPH for antioxidant defenses. Cell Rep. 2019;28(9):2306–16. https://doi.org/10.1016/j.celrep.2019.07.097.
Ko CW, Qu J, Black DD, Tso P. Regulation of intestinal lipid metabolism: current concepts and relevance to disease. Nat Rev Gastroenterol Hepatol. 2020;17(3):169–83. https://doi.org/10.1038/s41575-019-0250-7.
Jiang P, Sun W, Shen N, Huang X, Fu S. Identification of a metabolism-related gene expression prognostic model in endometrial carcinoma patients. BMC Cancer. 2020;20(1):864. https://doi.org/10.1186/s12885-020-07345-8.
Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug Discov. 2019;18(3):197–218. https://doi.org/10.1038/s41573-018-0007-y.
Lizardo DY, Kuang C, Hao S, Yu J, Huang Y, Zhang L. Immunotherapy efficacy on mismatch repair-deficient colorectal cancer: from bench to bedside. Biochim Biophys Acta Rev Cancer. 1874;2020(2):188447. https://doi.org/10.1016/j.bbcan.2020.188447.
Lan Y, Zhang D, Xu C, Hance KW, Marelli B, Qi J, et al. Enhanced preclinical antitumor activity of M7824, a bifunctional fusion protein simultaneously targeting PD-L1 and TGF-beta. Sci Transl Med. 2018;10(424):eaan5488.
Li S, Liu M, Do MH, Chou C, Stamatiades EG, Nixon BG, et al. Cancer immunotherapy via targeted TGF-beta signalling blockade in TH cells. Nature. 2020;587(7832):121–5. https://doi.org/10.1038/s41586-020-2850-3.
Zhao H, Yang L, Baddour J, Achreja A, Bernard V, Moss T, et al. Tumor microenvironment derived exosomes pleiotropically modulate cancer cell metabolism. Elife. 2016;5:e10250. https://doi.org/10.7554/eLife.10250.
Wang S, Liu R, Yu Q, Dong L, Bi Y, Liu G. Metabolic reprogramming of macrophages during infections and cancer. Cancer Lett. 2019;452:14–22. https://doi.org/10.1016/j.canlet.2019.03.015.
Phan AT, Goldrath AW, Glass CK. Metabolic and epigenetic coordination of T cell and macrophage Immunity. Immunity. 2017;46(5):714–29. https://doi.org/10.1016/j.immuni.2017.04.016.
Bailey JD, Diotallevi M, Nicol T, McNeill E, Shaw A, Chuaiphichai S, et al. Nitric oxide modulates metabolic remodeling in inflammatory macrophages through TCA cycle regulation and Itaconate accumulation. Cell Rep. 2019;28(1):218–30. https://doi.org/10.1016/j.celrep.2019.06.018.
Huang SC, Everts B, Ivanova Y, O'Sullivan D, Nascimento M, Smith AM, et al. Cell-intrinsic lysosomal lipolysis is essential for alternative activation of macrophages. Nat Immunol. 2014;15(9):846–55. https://doi.org/10.1038/ni.2956.
Wu H, Han Y, Rodriguez SY, Deng H, Siddiqui S, Treese C, et al. Lipid droplet-dependent fatty acid metabolism controls the immune suppressive phenotype of tumor-associated macrophages. Embo Mol Med. 2019;11(11):e10698. https://doi.org/10.15252/emmm.201910698.
Wang Y, Lu JH, Wang F, Wang YN, He MM, Wu QN, et al. Inhibition of fatty acid catabolism augments the efficacy of oxaliplatin-based chemotherapy in gastrointestinal cancers. Cancer Lett. 2020;473:74–89. https://doi.org/10.1016/j.canlet.2019.12.036.
Han S, Wei R, Zhang X, Jiang N, Fan M, Huang JH, et al. CPT1A/2-Mediated FAO Enhancement-A Metabolic Target in Radioresistant Breast Cancer. Front Oncol. 2019;9:1201.
We acknowledge TCGA and GEO database and contributors for uploading their meaningful datasets.
The study was supported by the grants of the Wu Jieping Medical Foundation (320.2710.1855).
Ethics approval and consent to participate
Both TCGA and GEO database are publicly available, thus ethical approval was not required for present study.
Consent for publication
All authors declare that there were no potential conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Supplementary Table 1.
Clinical characteristics of the training and validation set.
Additional file 2: Supplementary Table 2.
The lipid metabolism-related gene sets from the Molecular Signature.
Additional file 3: Supplementary Fig. 1.
Immune landscape of the risk signature.
Additional file 4: Supplementary Fig. 2.
Mutation landscape of the risk signature.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Yang, C., Huang, S., Cao, F. et al. A lipid metabolism-related genes prognosis biomarker associated with the tumor immune microenvironment in colorectal carcinoma. BMC Cancer 21, 1182 (2021). https://doi.org/10.1186/s12885-021-08902-5
- Colorectal cancer
- Lipid metabolism-related genes
- Prognostic value
- Tumor immune microenvironment