Data regarding TCGA-LIHC cohort from TCGA (https://portal.gdc.cancer.gov/) was downloaded using the TCGAbiolinks package (Colaprico et al., 2016). The RNA-seq data (normal: 50, tumor: 374) of HCC patients were obtained. The expression profiles were standardized in fragments per kilobase per million (FPKM) format and log2-transformed, while corresponding clinical and survival information of the HCC patients were collected from the UCSC Xena browser (http://genome.ucsc.edu) (Goldman et al., 2020). In this study, the TCGA-LIHC dataset was severed as the training set.

Additionally, we downloaded HCC related microarray data of GSE25097 (Ivanovska et al., 2011) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE25097) and GSE54236 (Zubiete-Franco et al., 2019) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54236) datasets from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database using the GEOquery (Davis and Meltzer, 2007) package. The microarray platform for the GSE25097 dataset was GPL10687. This datasets included transcriptome profiles of 268 HCC tumor tissue samples. The microarray platform for the GSE54236 dataset was GPL6480. This datasets included gene expression of 81 HCC tumor tissue samples. The GSE25097 and GSE54236 datasets were served as the verification datasets for analysis.

The LAMRGs were collected from the GeneCards (https://www.genecards.org/) (Stelzer et al., 2016) and MSigDB database (https://www.genecards.org/) (Liberzon et al., 2015). We used the keyword “lactic acid metabolism” as the search keyword to obtain 1,986 LAMRGs from the GeneCards database. The same keyword was used to obtain 194 genes LAMRGs from the MSigDB database. In addition, 22 LAMRGs were obtained from literature (Zhao et al., 2022). Finally, we intersected the LAMRGs from the three sources to obtain 2139 genes related to lactic acid metabolism (Supplementary Table S1).

The limma package was used to screen the DEGs from expression profile of TCGA-LIHC dataset. The threshold was set as | logFC | > 2 and p. adjust <0.05. To obtain the LAMRDEGs, all DEGs from TCGA-LIHC datasets and LAMRGs were intersected and displayed using a Venndiagram. Volcano plots of the differences were generated using the “ggplot2”package. The heatmap regarding the LAMRDEGs was drawn using R package “pheatmap.”

The R package “clusterProfiler” (Yu et al., 2012) was used to conduct Gene Ontology (GO) analysis (Mi et al., 2019) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) enrichment analysis regarding the LAMRDEGs. The selection criteria were set as p. adjust< 0.05 and false discovery rate (FDR) < 0.25. The Benjamini-Hochberg method was used for p correction.

The “c2. cp.all.v2022.1. Hs.symbols.gmt” gene sets were obtained from the MSigDB database (Liberzon et al., 2015). The “clusterprofiler” package was used for the GSEA (Subramanian et al., 2005) analysis of the DEGs obtained from TCGA-LIHC dataset. Furthermore, the significant enrichment standards were p. adjust< 0.05 and FDR <0.25.

We identified 14 key prognostic-related genes from the LAMRDEGs using Cox regression model. The expression of the 14 key genes was then included in the univariate Cox regression model and a forestplot was drawn. These 14 key genes were included in the multivariate Cox regression model, and a nomogram model was established to predict the 1-year survival of patients with HCC. Finally, a calibration curve was constructed to evaluate the prediction ability of the model. The “rms” package was used to construct the nomogram and calibration curves. We used the R package ggDCA (Tataranni and Piccoli, 2019) to conduct decision curve analysis (DCA) to evaluate the effect of the nomogram model in predicting 1-, 3-, and 5-year survival outcomes in patients with HCC. The formula for calculating the prognostic risk scores was: risk S c o r e = ∑ i C o e f f i c i e n t g e n e i ∗ m R N A E x p r e s s i o n g e n e i

We obtained the gene set “H.A.v7.4. Symbols.gmt” from the MSigDB database (Liberzon et al., 2015), and conducted GSVA (Hänzelmann et al., 2013) regardingthe HCC samples in TCGA-LIHC dataset. We divided the HCC samples into high- and low-risk groups using a Cox model. The functional differences in the enrichment pathways between the high- and low-risk groups were analyzed.

To obtain the model of key genes in TCGA-LIHC dataset, we used the glmnet package (Yang et al., 2022) based on 14 key genes to conduct a LASSO regression (Cai and van der Laan, 2020). Through the LASSO model, we screened six key genes, visualized the results of the LASSO regression and presented the expression levels of each gene between the two groups in the diagnostic model.

The STRING database (Szklarczyk et al., 2015) was used to construct a protein-protein interaction (PPI) network using hub genes. In addition, we used the Comparative Toxicogenomics Database (CTD) (Davis et al., 2021) (http://ctdbase.org/) to predict potential drugs or small molecule compounds that could interact with hub genes. Cytoscape software was used to visualize the mRNA-drug interaction network.

We excluded samples with no prognostic information and obtained expression profile data from TCGA-LIHC dataset. The information on the high- and low-risk groups in TCGA-LIHC dataset was obtained using a risk score model. The gene expression matrix of the dataset was uploaded to CIBERSORTx (https://cibersortx.stanford.edu) (Steen et al., 2020). Combined with the immune cell marker matrix, an immune cell infiltrating matrix was obtained (filter standard, p < 0.05). We then filtered out the data with an enrichment score of immune cells >0 and obtained the immune cell infiltration matrix. Finally, we acquired immune cells with significant differences in infiltration between the high- and low-risk groups using the LASSO model and plotted a heatmap and accumulation histogram of the correlation between the six key genes and these immune cells.

All data processing and analyses were performed using the R software (version 4.1.2). To compare two groups of continuous variables, the independent Student’s t test was used to estimate the statistical significance of normally distributed variables. The Mann-Whitney U test (Wilcoxon rank test) was used to analyze the differences between non-normally distributed variables. The Kruskal-Wallis test was used to compare three or more groups. Furthermore, the Chi-square test or Fisher’s exact test was used to compare the differences between categorical variables. Correlations between the different genes were assessed using Spearman’s correlation analysis. Kaplan-Meier (K-M) analysis was used to evaluate the survival value of the key genes. The level of statistical significance was set at a p-value <0.05.

The flow chart of this study is presented in Figure 1. The data from TCGA-LIHC dataset were divided into HCC and normal groups. Ninety-four DEGs between the two groups were identified using the limma package. The filtering threshold was set at |logFC| > 2 and p. adjust< 0.05. There were 15 upregulated and 79 downregulated genes (Figures 2A, B). The LAMRGs were obtained from the GeneCards database and literature. To obtain LAMRDEGs, we intersected genes from the LAMRGs and the above DEGs. We then identified five upregulated (CCNB1, E2F1, NQO1,LCN2,SPP1) and 23 downregulated (ADAMTS13, ALPL, BCHE, CFP, CXCL12, CYP2A6, CYP2C9, DBH, FBP1, FOS, GPD1, GYS2, HBB, HP, IGF2, IGFBP3, KDM8, LCAT, PCK1, RCAN1, SERPINE1, SHBG, and TDO2) LAMRDEGs.

The venndiagram and heatmap were constructed to illustrate the LAMRDEGs between the HCC and normal tissues (Figures 2C–E).

GO and KEGG enrichment analyses of the LAMRDEGs were conducted. The results are illustrated in Figures 3A, B. Subsequently, we conducted GO and KEGG enrichment analyses of the LAMRDEGs combined with the logFC values. The results are presented in Figures 3C, D. As shown in Figure 3A-D, the LAMRDEGs were mainly enriched in biological processes, including response to xenobiotic stimulus (GO, 0009410), response to glucocorticoid (GO, 0051384), response to corticosteroid (GO, 0031960). The LAMRDEGs were significantly enriched in the secretory granule lumen (GO, 0034774), cytoplasmic vesicle lumen (GO, 0060205), and vesicle lumen (GO, 0031983) of GO cellular components (CC). The LAMRDEGs were also enriched in heme binding (GO, 0020037), tetrapyrrole binding (GO, 0046906), integrin binding (GO, 0005178) of molecular functions (MF). In addition, the KEGG pathway analysis showed that the pathways were enriched in cellular senescence (hsa04218), and the p53 signaling pathway (hsa04115) (Table 1).

To further elucidate the potential mechanisms of the DEGs in HCC, we used GSEA to explore prominent KEGG pathways in the expression of DEGs from TCGA-LIHC dataset. The enrichment results are presented as a ridgelineplot (Figure 4A). The results showed that the DEGs were significantly enriched in different pathways (Table 2), including oxidative stress induced senescence (Figure 4B), glycolysis gluconeogenesis (Figure 4C), fatty acid metabolism (Figure 4D), and the defective intrinsic pathway for apoptosis (Figure 4E).

The univariate Cox regression analysis was used to assess the prognostic value of LAMRDEGs. The screening revealed that 14 key genes (E2F1, SERPINE1, GYS2, SPP1, PCK1, CCNB1, CYP2C9, IGFBP3, KDM8, RCAN1, ALPL, FBP1, NQO1, and LCAT) were associated with the over survival (Table 3; Figure 5A). We constructed a nomogram to assess the predictive value of the Cox model using these genes (Figure 5B). In addition, we constructed 1-year, 3-year, and 5-year (Figure 5C) calibration curves regarding nomogram. The results indicated that the prediction effect of the model in the 5-year survival analysis was better than that in the 1-year and 3-year survival analyses. Moreover, we used the DCA to evaluate the clinical utility of the prognostic model in terms of the 1-year (Figure 5D), 3-year (Figure 5E), and 5-year (Figure 5F). The results suggested that the prediction effect of the model in the 5-year survival analysis was better than that in 1- and 3-year survival analyses.

In addition, we divided TCGA-LIHC and GSE25097 samples into high- and low-risk groups according to the risk scores and compared expression of 14 key genes between the high- and low-risk groups from TCGA-LIHC dataset (Figure 5G). In the GSE25097 dataset, only the expression levels of 12 key genes (E2F1, SERPINE1, GYS2, SPP1, PCK1, CCNB1, CYP2C9, RCAN1, ALPL, FBP1, NQO1, and LCAT) differed significantly between the high- and low-risk groups (Figure 5H).

To explore the differences in the hallmark gene set between the high- and low-risk groups, GSVA was performed on TCGA-LIHC dataset. The results showed that different genes between the high- and low-risk groups in the LIHC samples were significantly enriched in 17 pathways, of which 13 were hallmark gene sets (angiogenesis, coagulation, complement, epithelial mesenchymal transition, G2M checkpoint, IL2/STAT5 signaling, inflammatory response, KRAS signaling up, myogenesis, PI3K/AKT/mTOR signaling, reactive oxygen species pathway, TGF-β signaling, and UV response DN) (Figure 6A; Table 4). A comparison of the hallmark gene sets is presented in Figure 6C. A correlation heatmap of the 14 key genes is presented in Figure 6B.

We screened key prognosis-related genes using LASSO regression (Figures 7A, B) and constructed a predictive model for prognosis. The prognostic model comprised six key genes: SERPINE1, SPP1, CCNB1, CYP2C9, NQO1, and LCAT (Table 5). A risk scoreplot for the six key genes is illustrated in Figure 7C. Finally, according to the risk score of the predictive model, the samples from TCGA-LIHC and GSE54236 datasets were divided into high- and low-risk groups. The expression levels of the six key genes between the high-and low-risk groups from TCGA-LIHC (Figure 7D) and GSE54236 (Figure 7E) were compared.

The Kaplan-Meier survival analysis was performed to examine the prognostic value of the six key genes (SERPINE1, SPP1, CCNB1, CYP2C9, NQO1, and LCAT) in the high expression and low-expression groups (Figures 8A–F). The results showed that the high expression of SERPINE1 (Figure 8A), SPP1 (Figure 8B), CCNB1 (Figure 8C) and NQO1 (Figure 8E) and low expression of CYP2C9 (Figure 8D) and LCAT (Figure 8F) were associated with poor prognosis. In addition, 1-, 3- and 5-year time-dependent receiver operating characteristic curves were plotted to analyze the prognostic value of the six key genes in TCGA-LIHC dataset. The results showed that the expression levels of SERPINE1 (Figure 8G), SPP1 (Figure 8H), CCNB1 (Figure 8I), and NQO1 (Figure 8K) had prognostic predictive value in HCC patients, while CYP2C9 (Figure 8J) and LCAT (Figure 8L) had no significant prognostic value.

Since genes that regulate the same biological processes have close relationships, we further analyzed the interaction between the key genes. The STRING database was used to construct a PPI network of the six key genes (SERPINE1, SPP1, CCNB1, CYP2C9, NQO1, and LCAT) (Figure 9A). We found that NQO1 had the strongest positive correlation with SPP1, CYP2C9 had the strongest negative correlation with CCNB1, and SERPINE1 had the weakest correlation with the other genes. In addition, we constructed and predicted the functionally similar gene interaction network of six key genes using GeneMANIA and found that 20 genes were co-expressed with the six key genes (Figure 9B).As shown in the mRNA-drug interaction network constructed using the CTD database (Figure 9C), we identified 11 potential drugs or compounds for the six key genes.

Finally, chromosomal mapping was performed regarding the six key genes. The results showed that SPP1 was distributed on chromosome 4, CCNB1 on chromosome 5, SERPINE1 on chromosome 7, and CYP2C9 on chromosome 10. LCAT and NQO1 were distributed on chromosome 16 (Figure 9D).

The expression profile data and grouping information were collated and uploaded to the CIBERSORTx. The correlation between the 22 types of immune cells and the expression profiles of the high- and low-risk group samples were assessed. According to the analysis results, the distribution of immune cell infiltration in TCGA-LIHC dataset is illustrated in Figure 10A . The high-risk group had a higher abundance of macrophage, CD4⁺ memory active T cells and gamma-delta T cells (γδT), and had a lower abundance of mast cells. In addition, the correlation between immune cell infiltration abundance and the six key genes was analyzed, and the correlation heatmap is illustrated in Figure 10B.The positive correlation between CCNB1 and macrophages M0 was the top1 (p < 0.001; r = 0.316; Figure 10C). The negative correlation between CCNB1 and mast cells resting was the top1 (p < 0.001; r = −0.326; Figure 10D).