Abstract Hepatocellular carcinoma (HCC) is one of the most common malignant tumors. This study was aimed to identify a lipid metabolism-related signature associated with the HCC microenvironment to improve the prognostic prediction of HCC patients. Clinical information and expression profile data were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, including the GEO dataset [28]GSE76427. The gene expression profile of lipid metabolism was downloaded from Molecular Signatures Database (MSigDB) database. The infiltrating immune cells were estimated by the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE), MCP-counter, and TIMER algorithms. Functional analysis, including Gene ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene set enrichment analysis (GSEA) were performed to elucidate the underlying mechanisms. The prognostic risk model was performed by Least absolute shrinkage and selection operator (LASSO) algorithm and Cox regression analysis. Two distinct subgroups of survival were identified. Better prognosis was associated with high immune score, high abundance of immune infiltrating cells, and high immune status. GO and KEGG analysis showed that differentially expressed genes (DEGs) between the two subgroups were mainly enriched in immune related pathways. GSEA analysis suggested that the expression of lipid metabolism related genes (LMRGs) was related to dysregulation of immune in the high-risk group. Risk models and clinical features based on LMRGs predicted HCC prognosis. This study indicated that the lipid metabolism-related signature was important for the prognosis of HCC. The expression of LMRGs was related to the immune microenvironment of HCC patients and could be used to predict the prognosis of HCC. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-024-76578-5. Keywords: Hepatocellular carcinoma, Tumor microenvironment, Lipid metabolism, LMRGs Subject terms: Data mining, Data processing, Databases, Genome informatics, Software, Biomarkers Introduction Hepatocellular carcinoma (HCC) has a high incidence and accounts for the third highest number of cancer-related deaths in the world^[29]1. There are several treatments for hepatocellular carcinoma, including chemotherapy, radiofrequency ablation, liver transplantation and surgery. However, it still has a high mortality rate. Therefore, it is urgent to explore the characteristics of HCC to develop new therapeutic methods. The tumor microenvironment (TME) is a dynamic system composed of immune cells, mesenchymal cells, cancer cells, complex cytokine and chemokine environment^[30]2, which plays a crucial role in the occurrence and development of tumors^[31]2,[32]3. In the tumor microenvironment, tumor-infiltrating immune cells are the major non-tumor components and have been shown to play an important role in predicting the prognosis of patients. Therefore, TME is essential to inhibit the occurrence, invasion and metastasis of tumors and to effectively manage immune responses. Recently, reprogramming of lipid metabolism is considered as a new marker of tumor malignancy^[33]4. More and more studies showed that the disorder of lipid metabolism played an important role in the occurrence, development and treatment of tumors^[34]5,[35]6. In addition, previous studies showed that lipid metabolism related genes (LMRGs) had a potential prognostic potential in a variety of tumor types^[36]7,[37]8. Therefore, targeting lipid metabolism is a novel tumor therapy strategy, and the establishment of a prognostic risk model is a feasible strategy to evaluate the prognosis. Multiple risk models have been developed to investigate the prognostic value of genes associated with tumor microenvironment, immune cell infiltration, and energy metabolism^[38]9. However, the relationship between lipid metabolism and immune microenvironment in HCC remains unclear. In this study, the influence of lipid metabolism on survival time of HCC patients was investigated through comprehensive analysis of LMRGs. In addition, we established a risk score model based on LMRGs to evaluate the prognostic value of LMRGs for HCC. This study provides a new way to explore the molecular mechanism and develop therapeutic strategies of HCC. It is important to acknowledge that our study relies on publicly available datasets, which may introduce biases due to patient heterogeneity. Furthermore, experimental validation of the findings is necessary to confirm the clinical applicability of the proposed risk model. Despite these limitations, our study offers a robust framework for future research in lipid metabolism and tumor microenvironment interactions in HCC. Materials and methods Data collection and collation Clinical information and expression profile data were downloaded from The Cancer Genome Atlas^[39]10 ([40]https://cancergenome.nih.gov) and Gene Expression Omnibus (GEO)^[41]11 databases ([42]https://www.ncbi.nlm.nih.gov/geo/). Inclusion criteria were as follows: (a) samples diagnosed with HCC; (b) samples with clinical information maps and gene expression matrices; (c) at least complete clinical data, including survival time, survival status, age and sex. Exclusion criteria are as follows: (a) samples with incomplete clinical data; (b) normal tissue specimens; (c) samples with bias in performance values. (d) no expression of more than half of the genes. A lot of 115 GEO samples ([43]GSE76427)^[44]12served as the experimental queue and 310 TCGA samples served as the validation queue. Datasets of 742 LMRGs were obtained from the Molecular Signatures Database^[45]13 (MSigDB, [46]http://www.gsea-msigdb.org/gsea/msigdb). Gene expression data were normalized using log2(TPM + 1) transformation. Batch effects between datasets were corrected using the ComBat method in the R package. A flowchart of the investigation is depicted in Fig. [47]1. Fig. 1. [48]Fig. 1 [49]Open in a new tab The flowchart of this study. Acquisition of prognostic genes In this study, survival, a R package, was used to integrate survival time, survival status and gene expression data, and the prognostic significance of each gene was evaluated by univariate Cox regression analysis. Forty-one prognostic LMRGs were obtained. Immune analyses The infiltrating immune cells were estimated by the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE), MCP-counter, and TIMER algorithms^[50]14. TIMER focuses on six key immune cells, including B cells, CD4^+ T cells, CD8^+ T cells, dendritic cells, neutrophils, and macrophages. MCP-counter assesses the abundance of eight immune and stromal populations, such as T cells, CD8^+ T cells, cytotoxic lymphocytes, NK cells, B cells, monocytes, fibroblasts, and endothelial cells. ESTIMATE provides stromal and immune scores, offering a quantitative measure of stromal and immune cell infiltration. Together, these algorithms offer a comprehensive evaluation of immune cell infiltration in the tumor microenvironment. Functional annotation and enrichment analyses Differentially expressed genes (DEGs) were identified using Limma, an R package, with the screening criteria set as |log2(FC)| >1.5 and FDR < 0.05. Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using the ‘cluster Profiler (version 3.14.3)’ R package to enrich related pathways and visualized in Metascape database. Moreover, gene set enrichment analysis (GSEA) was used to analyze the difference in both risk groups. In this study, we selected a series of biologically important KEGG pathways for mapping, including map01100 Metabolic pathways, map00562 Inositol phosphate metabolism, map04070 Phosphatidylinositol signaling system, map03320 PPAR signaling pathway, map01212 Fatty acid metabolism, map00120 Primary bile acid biosynthesis, map00061 Fatty acid biosynthesis, map04664 Fc epsilon RI signaling pathway, map04146 Peroxisome, map04666 Fc gamma R-mediated phagocytosis^[51]15–[52]17. Establishment and verification of risk model Least absolute shrinkage and selection operator (LASSO) analysis (‘glmnet’ R package) was further conducted to downsize the candidate genes and establish the prognostic signature. And 23 genes were obtained. Risk score in the training cohorts was calculated as: RiskScore=-0.800174294887228*ME1-0.0193829285839926*PPT1 + 0.525213565672368*AACS-0.0632880430826397*ALOX5AP + 0.403786265154411*CSNK2A2 + 1.28123770533637*MED17-0.635686604105556*LGMN + 0.280306338382094*LIPH + 2.05978842340822*SMPD4-0.00615923541586774*UGCG + 0.173805803425654*MTMR9-0.962375286733725*GPX4-0.691208094966126*ARNT + 0.602012000516297*ESRRA + 2.02316300371252*ACSL6 + 0.564963574887128*OSBPL9 + 1.92836220858469*PITPNM2-0.0260981734112701*PLA2G1B + 0.701949656008654*PRKD1 + 0.472682631205348*MOGAT1-0.131410651385763*PMVK-0.624200606686374*FAR1 + 0.448175501350242*CSNK1G2. To determine the value of the constructed features in predicting outcomes, we plotted a survival curve using the ‘survival’ R package. In addition, ROC analyses were performed at 365, 1095, and 1825 time points using the pROC package to obtain the effectiveness of the final AUC results in predicting the prognosis of patients. Validation of the expression levels of the hub genes via the human protein atlas The Human Protein Atlas (HPA, [53]https://www.proteinatlas.org/) provides information on the distribution and expression of proteins in a variety of normal human tissues, tumor tissues, cell lines, and blood cells. We searched the expression of key genes in different hepatoma cell lines. To further verify the protein expression levels of ME1 encodes malidase 1 (ME1), Palmitoyl-protein thioesterase 1 (PPT1), and Legumain (LGMN) in liver cancer and normal tissues, immunohistochemistry (IHC) data was downloaded from the HPA. Statistical analysis R (version 4.2.1), GraphPad Prism (version 8.0) and sangerbox^[54]18 tools was used for data analysis and visualization. Student’s t-test was used between two groups, and one-way ANOVA analysis was used when there were three or more groups. All P values < 0.05 were considered statistically significant. Results Clinical characteristics between training and validation cohorts The study compared demographic and clinical characteristics between the training cohort (n = 310) and validation cohort (n = 115). Age distribution showed no significant difference (P = 0.1772), with 46.4% of the training cohort and 39.1% of the validation cohort under 60 years old. There was a significant difference in sex distribution (P = 0.0214), with 69.7% males in the training cohort and 80.9% in the validation cohort. Survival status also differed significantly (P = 0.0295), with 69.4% of the training cohort and 80% of the validation cohort deceased. TNM staging showed 71.6% of the training cohort and 78.3% of the validation cohort were in stage 1/2. In the training cohort, 30.6% had hepatitis B, 16.1% had hepatitis C, 5.5% had NAFLD, and 22.9% had cirrhosis (Table [55]1). Table 1. Clinical data of patients in GEO and TCGA databases. Training cohort (n = 310) Validation cohort (n = 115) P-value n/% n/% Age 0.1772 < 60 144/46.4% 45/39.1% >=60 166/53.6% 70/60.9% Sex 0.0214 Male 216/69.7% 93/80.9% Female 94/30.3% 22/19.1% Survival status 0.0295 Alive 215/69.4% 92/80% Dead 95/30.6% 23/20% TNM Stage 1/2 222/71.6% 90/78.3% Stage 3/4 73/23.5% 25/21.7% Discrepancy 1/0.4% / Not available 14/4.5% / Chronic liver disease Hepatitis B: 95/30.6% / Hepatitis C: 50/16.1% / Non-alcoholic Fatty liver disease 17/5.5% / Cirrhosis Yes 71/22.9% / No 239/77.1% / [56]Open in a new tab Acquisition and analysis of DEGs We identified 1239 different prognostic genes by univariate Cox analysis (Supplementary Table [57]1) and identified 742 LMRGs in the MSigDB database. Intersection targets were collected by Venn diagram, and 41 prognostic genes related to lipid metabolism were obtained (Fig. [58]2A). Subsequently, function and pathway analysis of 41 genes were performed. PPI network was constructed (Fig. [59]2B). PPI network and GO analysis results showed that genes were mainly enriched in lipid biosynthetic process, monocarboxylic acid metabolic process and lipid transport (Fig. [60]2C). PPI network and KEGG results showed that genes were mainly enriched in metabolic pathways, inositol phosphate metabolism, phosphatidylinositol signaling system and PPAR signaling pathway (Fig. [61]2D). Fig. 2. [62]Fig. 2 [63]Open in a new tab Analysis of DEGs. (A) Venn diagrams of DEGs. This panel shows the overlap of DEGs between different comparison groups. Each circle represents a set of DEGs identified from a specific condition or comparison, and the overlapping regions represent DEGs shared between the groups. (B) PPI network construction of DEGs. Nodes represent proteins, and edges represent predicted functional associations between them, providing insights into the interaction patterns among the DEGs. (C) GO function analysis of DEGs. (D) KEGG pathway enrichment analysis of DEGs. This graph is based on information from the KEGG database ([64]www.kegg.jp/kegg/kegg1.html). Acquisition of risk model based on LMRGs Then, a risk signature model was established to evaluate the predictive value of LMRGs for HCC prognosis. LASSO analysis was applied to identified the potential 23 genes with optimal lambda value in the risk model (Fig. [65]3A) (Table [66]2). Three of these genes, including ME1, PPT1 and LGMN were showed by heat map (Fig. [67]3B). The results showed that there were significant differences in the expression of these three key genes in different HCC cell lines (Supplementary Fig. [68]1). In addition, the IHC results for the expression of ME1, PPT1 and LGMN proteins in the HPA database are shown in Supplementary Fig. [69]2. The expression of ME1 is low in normal tissues and elevated in liver cancer tissues. The expression of PPT1 is high in both normal and liver cancer tissues, and the expression of LGMN is higher in liver cancer tissues than in normal tissues. The established risk model successfully divided HCC patients into high and low-risk group. The expression levels of the three genes in the high-risk group were higher than those in the low-risk group, and the overall survival (OS) rate in the high-risk group was lower than that in the low-risk group (Fig. [70]3C). The receiver operating characteristic (ROC) curve evaluated the results of the risk model. Time-dependent ROC analysis showed that the established risk model had accurate prediction ability within 5 years, with 1-year, 3-year and 5-year Area Under ROC Curve (AUC) of 0.85, 0.84 and 0.90, respectively (Fig. [71]3D), which suggested the risk score signature were well predictive. Fig. 3. [72]Fig. 3 [73]Open in a new tab Construction of risk model in the training cohort. (A) LASSO analysis with minimal lambda. This plot represents the coefficient profiles of the LMRG-based prognostic genes across a range of lambda values. The top panel shows how the gene coefficients shrink as the lambda increases, and the bottom panel shows the partial likelihood deviance. The vertical dashed line indicates the lambda value that minimizes the deviance, representing the optimal model. (B) Distribution of risk score and survival status of HCC patients in the low and high risk groups. The top section displays the distribution of risk scores among the HCC patients, separated into low- (blue) and high-risk (red) groups. The middle section depicts patient survival status (alive = green, dead = black) along the time axis. The bottom heatmap shows the expression levels of three key genes (ME1, PPT1, LGMN) in the low- and high-risk groups. (C) Overall survival curve of the HCC patients in the two groups. Survival curves for the HCC patients in the low-risk (blue) and high-risk (red) groups. (D) Time-dependent ROC curve of the risk model. Time-dependent ROC curves at 1, 3, and 5 years are shown to assess the predictive accuracy of the risk model. The AUC values for each time point indicate high sensitivity and specificity in predicting the prognosis of HCC patients. Table 2. 23 LMRGs related prognostic genes. Name Full name Coefficient HR P-value ME1 NADP-dependent malic enzyme 0.000188 0.405202 0.000509 PPT1 Palmitoyl-protein thioesterase 1 0.002283 0.337324 0.002242 AACS Acetoacetyl-CoA synthetase 0.005474 7.233331 0.003901 ALOX5AP Arachidonate 5-lipoxygenase-activating protein 0.005149 0.528584 0.005758 CSNK2A2 Casein kinase II subunit alpha 0.00652 4.584595 0.007672 MED17 Mediator of RNA polymerase II transcription subunit 17 0.013388 42.05104 0.012088 LGMN Legumain 0.011403 0.428704 0.01262 LIPH Lipase member H 0.016302 56.29267 0.014762 SMPD4 Sphingomyelin phosphodiesterase 4 0.017657 114.4033 0.016043 UGCG Ceramide glucosyltransferase 0.016685 0.458272 0.017096 MTMR9 Myotubularin-related protein 9 0.020788 4.743746 0.017428 GPX4 Phospholipid hydroperoxide glutathione peroxidase 0.023942 0.333637 0.018181 ARNT Aryl hydrocarbon receptor nuclear translocator 0.013325 0.348481 0.018578 ESRRA Steroid hormone receptor ERR1 0.024407 2.973173 0.022396 ACSL6 Long-chain-fatty-acid–CoA ligase 6 0.028662 39.39252 0.025242 OSBPL9 Oxysterol-binding protein-related protein 9 0.030145 42.01394 0.028413 PITPNM2 Membrane-associated phosphatidylinositol transfer protein 2 0.033791 103.6533 0.0359 PLA2G1B Phospholipase A2 0.027966 0.605435 0.038197 PRKD1 Serine/threonine-protein kinase D1 0.047215 2.518711 0.039906 MOGAT1 2-acylglycerol O-acyltransferase 1 0.035093 2.310221 0.040747 PMVK Phosphomevalonate kinase 0.037869 0.423983 0.044161 FAR1 Fatty acyl-CoA reductase 1 0.02525 0.277131 0.045524 CSNK1G2 Casein kinase I isoform gamma-2 0.051574 15.79398 0.045638 [74]Open in a new tab Patients in the two subtypes showed different immune status Subsequently, immunoassays were performed to explore the immune differences between the two subtypes. We performed a GSEA analysis to assess the relative expression differences between the two pathways. GSEA analysis showed that T cell receptor signaling pathway, Fc epsilon Ri signaling pathway, Antigen processing and presentation and B cell receptor signaling pathway expressed lowly in the high-risk group (Fig. [75]4A). These results indicated that the expression of LMRGs was related to dysregulation of immune, which might be involved in the poor prognosis of HCC patients. Fig. 4. [76]Fig. 4 [77]Open in a new tab Immune analyses in the two subgroups. (A) GSEA diagram visualizing GSEA analysis results. The GSEA diagram shows significant enrichment of pathways related to immune signaling, including T cell receptor signaling, B cell receptor signaling, and antigen processing and presentation. These pathways were enriched in the high-risk group compared to the low-risk group. (B) Stromal score, immune score and ESTIMATE score calculated by using ESTIMATE algorithm. (C) The abundance of immune infiltrating cells was evaluated by TIMER algorithm. Box plots show the abundance of various immune cell types, including B cells, CD4^+ T cells, CD8^+ T cells, neutrophils, macrophages, and dendritic cells. (D) The abundance of immune infiltrating cells was evaluated by MCP Counter algorithm. Violin plots show the expression levels of various immune cells, including CD8^+ T cells, cytotoxic lymphocytes, B cells, NK cells, monocytic cells, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts. The ESTIMATE algorithm showed that the low-risk patients had significantly higher immune score (P < 0.001), ESTIMATE score (P < 0.01) compared with the control group, with no significant difference found in stromal score (P = 0.21) (Fig. [78]4B). Furthermore, TIMER algorithm showed that the abundance of B cell (P < 0.05), CD4 T cell (P < 0.01), neutrophil (P < 0.05), macrophage (P < 0.05) and dendritic cell (P < 0.01) in the low-risk group was significantly higher than the high-risk group, while no statistical significance was showed with respect to CD8 T cell (P = 0.23) (Fig. [79]4C). In addition, the MCP Counter algorithm showed that the low-risk group had higher T cells (P = 1.4e-3), monocytic cell (P = 2.6e-3) and myeloid_ dendritic_ cells (P = 1.1e-3) (Fig. [80]4D). These results indicate that there were significant differences in immune status between the two subtypes. Correlation between prognostic models and clinicopathological features In addition, the relationship between risk scores and clinical features was explored, and the independence of the constructed risk models was assessed by subgroup analysis and regression analysis. There was no significant difference in the risk score between gender (Fig. [81]5A), age (Fig. [82]5B) and TNM stage (Fig. [83]5C), indicating that the risk score was not significantly associated with clinical features. In addition, the survival curve showed that when patients were grouped by gender (Fig. [84]5D-E), age (Fig. [85]5F-G) and TNM stage (Fig. [86]5H-I), the risk model still had strong predictive power, and patients with lower risk scores had better prognosis. These results indicated that the established risk model could better predict the prognosis of HCC patients. Fig. 5. [87]Fig. 5 [88]Open in a new tab Correlation between risk scores and clinical features in the experimental cohort. (A-C). Risk score distribution in relation to gender, age, and TNM stage: (A) Scatter plot showing that no significant difference (ns) was found between the risk scores of male and female HCC patients. (B) Scatter plot illustrating a significant difference (*) in risk scores between patients younger than 60 years and those aged 60 years or older. (C) Scatter plot demonstrating no significant difference (ns) in risk scores across different TNM stages (I-IV) of HCC. D-I. Survival curves stratified by clinical features: (D, E). Kaplan-Meier survival curves of HCC patients stratified by gender (male and female), comparing low- and high-risk groups. (F, G) Kaplan-Meier survival curves for patients younger than 60 years and those aged 60 years or older, showing the survival probability over time for the low- and high-risk groups. (H, I). Kaplan-Meier survival curves for patients with early-stage (I, II) and late-stage (III, IV) TNM classification, comparing the survival outcomes of low- and high-risk patients. Correlation between prognostic models and clinicopathological features in the verification cohort To verify the correlation between risk score and clinical features, TCGA data was used as validation. Similar result was also observed in the verification cohort. There was no significant correlation between risk scores and different gender (Fig. [89]6A) and age (Fig. [90]6B), while in TNM stage, the risk of stage II and III patients was higher than that of stage I patients (Fig. [91]6C). In addition, the survival curve showed that when patients were grouped by gender (Fig. [92]6D-E), age (Fig. [93]6F-G) and TNM stage (Fig. [94]6H-I), the risk model still had strong predictive power, and patients with lower risk scores had better prognosis. These results indicated that the established risk model could better predict the prognosis of HCC patients. Fig. 6. [95]Fig. 6 [96]Open in a new tab Correlation between risk scores and clinical features in the verification cohort. (A-C) Risk score distribution in relation to gender, age, and TNM stage: (A) Scatter plot showing that no significant difference (ns) was found between the risk scores of male and female HCC patients. (B) Scatter plot illustrating a significant difference (*) in risk scores between patients younger than 60 years and those aged 60 years or older. (C) Scatter plot demonstrating no significant difference (ns) in risk scores across different TNM stages (I-IV) of HCC.No significant difference was identified in patients with different gender (A), age (B), and TNM staging (C). (D-I). Survival curves stratified by clinical features: (D, E) Kaplan-Meier survival curves of HCC patients stratified by gender (male and female), comparing low- and high-risk groups. (F, G). Kaplan-Meier survival curves for patients younger than 60 years and those aged 60 years or older, showing the survival probability over time for the low- and high-risk groups. (H, I). Kaplan-Meier survival curves for patients with early-stage (I, II) and late-stage (III, IV) TNM classification, comparing the survival outcomes of low- and high-risk patients. Discussion A lot of studies evidence suggested that disorder of tumor lipid metabolism led to tumor progression and local immunosuppression^[97]19. Due to the limited clinical efficacy of liver transplantation, chemotherapy, surgical excision and other methods, it was urgent to find new biomarkers, which could predict the survival of HCC patients, but also be used to guide anti-cancer therapy^[98]8. Recent studies have deepened our understanding of the impact of abnormal lipid metabolism on the immune microenvironment in liver cancer patients, shedding light on its relevance for prognosis and therapeutic development^[99]8. Our study highlights the significant impact of lipid metabolism disorders on HCC progression and the immune microenvironment. The identification of distinct HCC subtypes based on lipid metabolism characteristics and the development of our LMRG-based prognostic model contribute valuable insights into the role of lipid metabolism in cancer prognosis. Previous studies have utilized similar bioinformatics approaches to develop potential cancer biomarkers, which support and contextualize our findings. For example, one study developed a novel CTA-related feature for gastric cancer, identifying ELOVL4 as a potential therapeutic target^[100]20. Another study developed a novel cancer-associated fibroblast-related signature to predict clinical outcomes and immune responses in cervical cancer^[101]21. In this study, we identified two distinct subtypes of HCC patients that exhibited significant differences in lipid metabolism characteristics. Immunoassay analyses revealed that patients with poor prognosis were characterized by lower immune scores, lower ESTIMATE scores, and overall immunosuppression. Further functional analysis showed that the expression of LMRGs was related to dysregulation of immune in the high-risk group. In addition, we established a prognostic risk model for HCC based on LMRGs, which could accurately predict the prognosis of HCC patients. This model provides clinicians with a valuable tool to make more informed treatment decisions, particularly in terms of targeted therapy for HCC. The prognostic risk score generated by our model enhances the accuracy of survival prediction and could be instrumental in guiding clinical decision-making. A substantial number of genes have previously been identified as potential cancer biomarkers, but our study highlights the importance of integrating mathematical models into patient risk assessment^[102]22. By combining survival time, survival status, and gene expression data, we first identified 41 survival-related genes using the Cox method. Through further intersection analysis of LMRGs and HCC survival-related genes, we identified lipid metabolism-related DEGs, including ME1, PPT1, and LGMN. ME1 encodes malic enzyme 1, which catalyzes the conversion of malic acid to pyruvate and contributes to NADPH formation. Previous studies reported that ME1 was mainly involved in lipid biosynthesis to mediate lipid metabolism and promoted the development of various tumors^[103]7,[104]23,[105]24. Research has shown that ME1 expression is upregulated in liver cancer cells and human tissues, with patients exhibiting high expression levels having lower survival rates. Moreover, ME1 knockdown inhibits the migratory and invasive properties of HCC cells^[106]25. Similarly, PPT1 catabolized lipid-modified proteins^[107]26, which were significantly upregulated in HCC tissues and were significantly associated with a poor prognosis^[108]27. Research has shown that PPT1 levels are significantly upregulated in HCC tissues both in vivo and in vitro, and are strongly associated with poor prognosis. PPT1 knockdown enhances antitumor immunity^[109]28.LGMN promoted gastric cancer cell proliferation and angiogenesis through tumor-associated macrophages^[110]29. Knockdown of LGMN expression inhibits tumor growth and angiogenesis^[111]29. In line with these findings, our study further revealed the expressions of ME1 and LGMN were significantly increased in liver cancer. Furthermore, there were significant differences in the expression of these three key genes in different HCC cell lines. The expression of ME1 is low in normal tissues and elevated in liver cancer tissues. Next, GO analysis and KEGG analysis was performed to explore their potential biological mechanisms. The results suggested that dysregulation of immune regulation mediated the role of lipid metabolism in the occurrence and development of HCC. We also performed a GSEA analysis to further clarify the underlying mechanism. GSEA analysis directly clarify the expression trend of gene sets in different populations. In this study, GSEA results showed that the expression of immune cell differentiation was relatively low in the high-risk group. These results suggested that lipid metabolism disorders were related to immune disorders. The TIME is pivotal in determining patient prognosis, as tumor progression is often associated with changes in the surrounding stroma, with immune cells playing a key role in this stroma. Abnormal metabolic status of tumor cells also leads to alterations in TIME^[112]30. ESTIMATE algorithm is an innovative method for inferring the fraction of immune and stromal cells in tumors based on gene expression values. The immune score generated by the ESTIMATE algorithm quantitatively indicates the immune component in the tumor sample^[113]31. The results showed that patients with better prognosis had higher immune scores. In addition, two methods, TIMER and MCP Counter, were used to evaluate the immune status of these two molecular subpopulations. TIMER is a web-based tool that facilitates the quantification of six tumor infiltrating immune subsets^[114]32. TIMER analysis showed that the abundance of B cell, CD4^+ T cell, neutrophil, macrophage and dendritic cell in the low-risk group was significantly higher than the high-risk group. In addition, the MCP Counter algorithm showed that the low-risk group had higher T cells, monocytic cell and myeloid_ dendritic_ cells. These findings suggest that lower immune scores and suppressed immune status are associated with worse prognoses in high-risk patients. Subgroup analyses based on clinical features, including gender, age, and TNM stage, showed strong agreement between the risk scores and disease stages, further supporting the robustness of our prognostic model. These results emphasize the potential of our LMRG-based risk model in predicting HCC patient outcomes across diverse clinical subgroups. Despite the significance of our findings, several limitations exist in this study. First, the absence of detailed clinical data, such as tumor staging and comorbidities, may have impacted the accuracy of our model. For instance, patients in the high-risk group may have had more advanced disease stages, potentially confounding our results by skewing the observed survival outcomes. Second, although our findings are based on comprehensive bioinformatics analyses, experimental validation is necessary to confirm the functional roles of key LMRGs in HCC progression. Specifically, gene-knockdown studies and in vivo models will be essential to elucidate the biological roles of these genes in lipid metabolism and immune regulation. Furthermore, larger cohorts with complete clinical datasets should be used to verify the prognostic significance of these genes across different HCC stages. In summary, our study identified two molecular subtypes of HCC, categorized into low-risk and high-risk groups, with distinct lipid metabolism profiles. Patients in the high-risk group, characterized by immune suppression, exhibited poorer survival outcomes. Our LMRG-based risk model effectively predicts HCC prognosis and highlights the association between lipid metabolism and immune dysregulation. However, the lack of clinical data and experimental validation, including gene knockdown experiments, limits the current interpretation of our findings. Future studies should aim to address these limitations, incorporating experimental models and clinical datasets to explore potential therapeutic interventions targeting lipid metabolism and immune pathways in HCC. Electronic supplementary material Below is the link to the electronic supplementary material. [115]Supplementary Material 1^ (147.7KB, docx) [116]Supplementary Material 2^ (963.5KB, tif) [117]Supplementary Material 3^ (4.7MB, tif) [118]Supplementary Material 4^ (15.7KB, docx) [119]Supplementary Material 5^ (14.7KB, docx) Acknowledgements