The HCC genomic data GSE46408 were obtained from the Gene Expression Omnibus database (ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE46408; Fig. 1) (19). mRNA microarray data were retrieved and the limma package (version 3.50.0) in R was used to identify differential expression of genes between HCC and normal liver tissue (23). For clinical differential expression studies, clinical data from The Cancer Genome Atlas (TCGA; portal.gdc.cancer.gov/) and Genotype-Tissue Expression (GTEx; gtexportal.org/) databases included 365 liver cancer and 160 normal tissue samples which were obtained using the University of California Santa Cruz Xena platform (xena.ucsc.edu/) datasets cohorts ‘GDC TCGA Liver Cancer (LIHC)’ (ID: TCGA-LIHC.htseq_counts.tsv) and ‘GTEX’ (ID: gtex_RSEM_Hugo_norm_count) (24–26).

Kaplan-Meier survival analysis was used to identify the association between PGK1 expression levels and the overall survival (OS) of patients with HCC. The R package TSHRC (version 0.1-6; cran.r-project.org/web/packages/TSHRC/index.html), was used to analyze data using the two-stage test with default settings (27). HCC samples were divided into high and low PGK1 expression groups and the mean levels of PGK1 in all samples were set as the cut-off value. A total of 365 samples from patients with HCC were obtained from the TCGA database, which included high expression of PGK1 in 187 cases and low expression levels of PGK1 in 178 cases. For the clinical characteristics analyses, samples from 365 patients with HCC were evaluated.

The clusterProfiler package (version 4.2.2) in R was used to perform GO functional annotation and KEGG pathway analysis (28). Using clusterProfiler, the GO terms were obtained with a setting of P≤0.01 and q≤0.05 considered to indicate a significant enrichment. The ggplot2 package (version 3.4.0) in R was used to plot the bar and bubble charts of enrichment results (29). STRING (string-db.org/) (version 11.5) was used to screen the protein-protein interaction (PPI) (30). Cytoscape software (version 3.8.0) was used to visualize the interaction network (31).

Human liver cancer cell lines (SNU182, SNU449, JHH5, HuH7, HepG2 and HCCLM3) were obtained from The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences. All cells were cultured in DMEM (HyClone; Cytiva) with 10% FBS (Gibco) in a 5% CO₂ incubator at 37°C.

The cell lines were authenticated using short tandem repeat (STR) profiling at the beginning and end of the study. STR profiling of SNU182 and HCCLM3 cell lines was performed Suzhou Jianda Biotechnology Co., Ltd. using the ABI Prism^® 3130 XL Genetic Analyzer (Applied Biosystems; Thermo Fisher Scientific, Inc.). STR profiling of HuH7 and JHH5 cell lines was performed by Shanghai Biowing Applied Biotechnology Co., Ltd., STR profiling of SNU449 cell line was performed by Guangzhou Jennio Biotechnology Co., Ltd. and STR profiling of HepG2 cell line was performed by Genesky Bio-Tech Co., Ltd. using the ABI Prism^® 3730 XL Genetic Analyzer (Applied Biosystems; Thermo Fisher Scientific, Inc.). STR profiling was performed according to the American Type Culture Collection (ATCC) recommended procedure (32) and the STR profiles matched the known ATCC STR profile of the cells (atcc.org/search-str-database).

Cells were digested and lysed on ice in RIPA buffer (Shanxi ZHHC Biotechnology Co., Ltd.; PL001-2A) with a protease inhibitor PMSF (Shanxi ZHHC Biotechnology Co., Ltd.; PL012-1), then quantified by BCA assay. Subsequently, 4–12% SDS-PAGE was used to separate 10 µg/lane protein lysate, which was then transferred to a PVDF membrane (MilliporeSigma). The membrane was blocked with 5% fat-free milk at RT for 1 h, then incubated with PGK1 (1:1,000; ABclonal Biotech Co., Ltd.; A12686) or GAPDH (1:50,000; ABclonal Biotech Co., Ltd.; A19059) antibodies at 4°C overnight. The membranes were incubated with anti-rabbit IgG, HRP-linked secondary antibody (1:5,000; CST; 7074S; www.cellsignal.cn) at RT for 1 h. Next, the membrane was washed using TBS (Guangzhou Roles-BIO Co., Ltd.; RBG6-1; www.rolesbio.com) with 0.05% Tween, visualized with Ultra-sensitive ECL chemiluminescent substrate (Biosharp Life Sciences; BL523B) and the protein expression levels analyzed by ImageJ software (version 1.52v; National Institutes of Health).

DEGs identification based on the GEO database was performed using the limma package in R and expression values were subjected to log2 transformation. Changes in gene expression levels >1.5-fold or <-1.5-fold and with P<0.05 were regarded as statistically significant, and the P-value was adjusted by the false discovery rate method. Expression analysis of clinical data based on TCGA and GTEx database was performed using an unpaired independent samples t test for clinical differential expression studies. Clinical and pathological classification data were analyzed by non-parametric Kruskal-Wallis followed by post hoc Dunn's test. The Dunn's test was calculated by the FSA package (https://fishr-core-team.github.io/FSA/) (version 0.9.5) in R, and the P-value was adjusted by the Bonferroni method. OS times were calculated by Kaplan-Meier method and log-rank test was used to assess statistical significance. Where late-stage crossover in Kaplan-Meier curves occurred, the two-stage test was used. Clinical association variables analysis based on TCGA database was evaluated by univariate Cox regression, χ² or Fisher's exact test, ‘Unknown’ groups were not included in the analysis. The mean level of PGK1 was set as the cut-off value in OS, Cox regression model, and χ² or Fisher's exact test analyses. Statistical analysis was computed using SPSS software (version 26.0; IBM Corp.) and R (version 4.0.4; RStudio, Inc.) software. The association between PGK1 mRNA expression levels and clinical and pathological stage was computed by the Kaplan-Meier plotter tool, with the condition set as default. The cut-off value was the best-performing threshold value, follow up threshold were all included and all other influences were contained all to remain constant except for the variables (33). PPI network was analyzed and screened by Cytoscape (version 3.8.0). Enrichment analysis of GO and KEGG annotations was performed using clusterProfiler in R software, and the P-value was adjusted by the Benjamini-Hochberg method. P<0.05 was considered to indicate a statistically significant difference. All experiments were performed with three replicates for each protein and data were expressed as the mean ± standard deviation. The results were analyzed by ImageJ software (version 1.52v; National Institutes of Health) and displayed by GraphPad Prism software (version 8.0; Dotmatics).

HCC-associated DEGs were identified based on the GSE46408 dataset, which included six pairs of HCC and their corresponding normal tissue samples. A total of 4,841 DEGs were identified; 2,606 were significantly up- and 2,235 were significantly downregulated in HCC compared with normal tissue (Fig. 2B; Table SI). Regarding genes involved in biological processes (BPs), the terms that appeared most frequently were those associated with metabolic processes such as ‘organic acid catabolic process’ (GO:0016054; 178 genes), ‘carboxylic acid catabolic process’ (GO:0046395; 178 genes) and ‘small molecule catabolic process’ (GO:0044282; 236 genes; Fig. 2B; Table SII). The cellular component (CC) genes were mainly in ‘chromosomal region’ (GO:0098687; 143 genes) and ‘mitochondrial matrix’ (GO:0005759; 173 genes; Fig. 2C) and the main DEGs involved in molecular function (MF) were small molecule binding and catalytic activity such as ‘flavin adenine dinucleotide binding’ (GO:0050660; 48 genes) and ‘vitamin binding’ (GO:0019842; 68 genes; Fig. 2D). The KEGG pathway results were mainly enriched in ‘Coronavirus disease-COVID-19’ (hsa05171; 89 genes), ‘cell cycle’ (hsa04110; 66 genes), ‘carbon metabolism’ (hsa01200; 51 genes), ‘complement and coagulation cascades’ (hsa04610; 46 genes) and ‘peroxisome’ (hsa04146; 41 genes; Fig. 2E; Table SIII). Taken together, these results demonstrated that DEGs were mainly enriched in the metabolism and energy metabolism pathways, which were selected for further research.

A total of 208 DEGs were associated with carbohydrate metabolism and derived from carbohydrate-associated KEGG entries after de-duplication, including 36 significantly up- and 172 significantly downregulated genes (Fig. 3A and B; Tables SIV and SV). In terms of the metabolic processes and pathways of HCC, these 208 DEGs based on carbohydrate metabolism were selected through GO functional annotation. For BP terms, the metabolic processes were mostly enriched in the ‘small molecule catabolic process’ (GO:0044282; 105 genes), ‘organic acid catabolic process’ (GO:0016054; 82 genes) and ‘carboxylic acid catabolic process’ (GO:0046395; 82 genes; Fig. 4A). The CC terms were mainly clustered in the ‘mitochondrial matrix’ (GO:0005759; 70 genes; Fig. 4B) and MF terms were largely observed in small molecule binding and catalytic activity such as ‘vitamin binding’ (GO:0019842; 34 genes) and ‘oxidoreductase activity, acting on CH-OH group of donors’ (GO:0016614; 29 genes; Fig. 4C; Table SVI). The upregulated genes were further investigated (Fig. 4D. Next, a protein interaction analysis was conducted to determine the functions of DEGs through the PPI network (Figs. 4E and S1; Tables SVII and SVIII).

To determine levels of PGK1 in the tissues of patients with HCC, analysis of the expression of PGK1 in clinical HCC tissues and normal tissue samples based on TCGA and GTEx databases was conducted. The patient sample IDs from TCGA and GTEx databases are listed in Table SIX. The expression levels of PGK1 were significantly increased in HCC compared with the normal tissue (Fig. 5A). To explore the association between expression of PGK1 and the prognosis of patients with HCC, OS time in the high and the low PGK1 group was compared using the Kaplan-Meier method based on TCGA database. Patients with high expression of PGK1 had a median survival of 3.1 years, which was significantly lower compared with the 5.1-year median survival time in the low PGK1 expression group (Fig. 5B).

The effect of treatment on OS time was analyzed using the Kaplan-Meier method (Fig. 5C; Table SX). A significant difference in OS between treatment groups was observed. The ablation embolization therapy group had the highest survival rate, with a mean survival time of 5.8 years. Patients who underwent radiation, combination, chemotherapy and targeted therapy had a median survival time of 4.4, 3.8 and 2.8 years, respectively. The surgical group of patients had the lowest median survival time of 1.1 years. Pairwise comparisons demonstrated that ablation embolization was significantly superior to chemo-targeted therapy and surgical monotherapy and chemo-targeted therapy was significantly superior to surgical monotherapy. However, other intergroup comparisons showed no statistically significant differences.

A nomogram was used to predict disease prognosis by integrating relevant variables. According to univariate Cox analysis, PGK1 and other independent prognostic factors (age, sex, ethnicity, clinical and TNM stage and histological grade) were selected for the construction of the nomogram. Staging was determined according to the 7th edition of the American Joint Committee on Cancer (AJCC) staging manual (34). PGK1 expression [hazard ratio (HR)=1.524; 95% CI, 1.079-2.154), clinical (HR=2.408; 95% CI, 1.663-3.486) and T (HR=2.501; 95% CI, 1.760-3.554) and M stages (HR=3.912; 95% CI, 1.230-12.441), were significantly associated with poor OS (Fig. 5D). To investigate the association between these factors affecting OS probability, PGK1 mRNA expression in 365 tumor samples taken from TCGA at different clinical pathological stages was compared. These results demonstrated a positive association between PGK1 expression and HCC clinical stage progression (from stage I to IV, H=14.642; Fig. 5E and F). There was significantly higher PGK1 expression in stage III tumors compared with stage I tumors. Similarly, PGK1 was significantly upregulated in the T stage, which increased from T1 to T4 based on status of the primary tumor (H=14.595). The expression levels of PGK1 increased significantly in T2 and T3 stages compared with the T1 stage. Moreover, to investigate its expression in HCC tumors in detail, the characteristics of the patients were explored. These data demonstrated that PGK1 expression was significantly associated with advanced tumor nodes and clinical stage. It was also demonstrated that PGK1 expression levels were associated with HCC pathological stage progression (Table I). The association between PGK1 expression and the specific clinical or pathological stage in prognosis was also demonstrated. Expression of PGK1 was significantly associated with the early stages (stage I and T1) of HCC (OS: HR=2.81 and 3.04, respectively; Table II). These results suggested that PGK1 expression levels were positively associated with progression of liver cancer and were negatively associated with survival prognosis of patients, which may be a risk factor in the malignant progression of HCC.

To verify the aforementioned predictions, validation tests were performed on HCC and paracancer tissue in addition to six liver cancer cell lines. Compared with paracancerous tissue, there was a markedly increased distribution of PGK1 in HCC tissue (Fig. 6A). Western blotting was performed to detect the protein expression of PGK1 in six liver cancer cell lines, SNU449, SNU182, HuH7, JHH5, HepG2 and HCCLM3 (Fig. 6B). HuH7 demonstrated the lowest protein expression levels of PGK1. By contrast, PGK1 was overexpressed in HepG2, JHH5, SNU449, SNU182 and HCCLM3 cells.