This article is mentioned in:

Introduction

Liver cancer is the 6th most common cancer worldwide and causes 830,000 deaths every year, with ~906,000 new cases diagnosed annually (1–5). As the most common form of liver cancer, hepatocellular carcinoma (HCC) is increasing in incidence worldwide and it has been predicted that HCC will cause at least 1 million deaths annually by 2030 (6,7). Although there are etiological agents responsible for HCC, such as hepatitis B and C virus infections, the molecular pathogenesis of HCC is still unclear (8). Nonetheless, cirrhosis and alcoholism, as well as metabolic syndrome and diabetes mellitus, are major etiologies of HCC globally (9). Therefore, unlike other types of solid malignancy, other complications or comorbidities should also be considered regarding HCC to identify biomarkers involved in its metabolic processes, which is important for diagnosis and increased understanding of tumor biology.

Due to lack of specific symptoms during the early stages of HCC, >60% of patients are diagnosed with advanced disease following metastasis, making the overall 5-year survival rate <20% (9–11). If patients are diagnosed with early-stage disease, the survival rate increases to >70% (12,13). Therefore, early detection of HCC is key to the success of tumor therapy administered when tumors are small to increase the survival rate of patients. Clinically, surveillance imaging combined with α-fetoprotein (AFP) and des-gamma-carboxyprothrombin (DCP) measurements are applied for early detection of HCC in patients (14). However, the roles of these detection methods in surveillance are controversial due to their high sensitivity for early diagnosis (15).

Over the past decade, the search for biomarkers has been part of a new era of ‘omics’ and has been promoted by a number of novel technologies, such as next-generation sequencing and microarray, which make tens of thousands of molecular targets analyzable and operable (16). Although numerous circulating biomarkers have been identified, few biomarkers have been applied clinically because of their low predictive accuracy and high cost (17,18).

The web-based opening genomic databases can not only provide gene expression profiles but can also provide paired normal tissue to screen for potential therapeutic targets for HCC (19,20). Glycolysis is key for metabolism of cells and is considered to be a good target for cancer therapy (21). Moreover, analyzing the glycolysis process in HCC has aided in understanding drug resistance mechanisms, pathogenesis and potential treatment paths for HCC (22). Therefore, the present study aimed to assess whether glycolysis signaling and the associated proteins in paired normal tissues adjacent to HCC can be a potential molecular prognostic marker.

Materials and methods

Expression analysis

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).

Figure 1. Flow chart of the experimental design. The GSE46408 HCC dataset (n=6 pairs of samples) was selected for analysis. TCGA served as an external validation dataset (n=365). PGK1, involved in carbohydrate metabolism, was identified in the TCGA training dataset by differential expression analysis, survival analysis and clinical characteristics studies. The risk-stratification ability of PGK1 and its association with tumor tissues were explored. Good predictive performance of the signature was confirmed in expression and clinical characteristics analyses. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; DEG, differentially expressed gene; GO, Gene Ontology; GEO, Gene Expression Omnibus; PPI, protein-protein interaction; GTEx, Genotype-Tissue Expression; PGK1, phosphoglycerate kinase 1.

Survival analysis and clinical characteristics studies

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.

Bioinformatics analysis

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).

Cell culture

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% CO2 incubator at 37°C.

Cell line authentication

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).

Western blotting

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).

Statistical analysis

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, χ2 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 χ2 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).

Results

Identification of DEGs and functional enrichment

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.

Figure 2. DEG selections and enrichment analysis of the GSE46408 dataset. (A) Volcano plot of DEGs. Red, upregulated genes (>1.5-fold); green, downregulated genes (<-1.5-fold); black, no significant differential expression. GO functional annotation of the DEGs. The DEGs (>|1.5-fold| and P<0.05) in hepatocellular carcinoma tissue compared with normal tissue involved in (B) biological processes, (C) cellular components and (D) molecular functions. The top 20 enriched GO terms of significant DEGs measured by adjusted P-value are shown. (E) Top 25 enriched KEGG terms of significant DEGs. DEG, differentially expressed gene; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of metabolic-associated DEGs in HCC

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).

Figure 3. DEGs associated with carbohydrate metabolism were identified in the GSE46408 dataset. The top 15 carbohydrate metabolism-associated KEGG terms and the associated DEGs were visualized by (A) pie chart and (B) bubble diagram. KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change; DEG, differentially expressed gene.

Figure 4. Enrichment analysis of DEGs associated with carbohydrate metabolism in the GSE46408 dataset. DEGs were measured by GO annotation analysis to extract terms associated with (A) biological process, (B) cellular components and (C) molecular function. (D) Heatmap showed upregulated DEGs of associated terms. (E) PPI network analysis demonstrated a number of associated terms. The interaction network of proteins was analyzed using the STRING database under the condition of high confidence score (0.700), which included 201 nodes and 1,132 edges with a PPI enriched P-value <1.0×10−16. The degree rank gradient from high to low was represented by the color gradient from blue to yellow. DEG, differentially expressed gene; GO, Gene Ontology; PPI, protein-protein interaction; STRING, search tool for the retrieval of interacting genes/proteins.

Analysis of PGK1 expression in patients with HCC

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).

Figure 5. PGK1 expression levels based on clinical data from TCGA and GTEx databases. (A) PGK1 expression levels based on HCC (n=365) and normal tissues (n=160) from TCGA and GTEx databases. Statistical significance was tested by independent samples t test. (B) Kaplan-Meier survival analysis of OS time of patients with HCC based on PGK1 expression. Statistical significance was determined by the two-stage method. OS analysis was considered from date of diagnosis to date of death. (C) Kaplan-Meier curves of OS time of patients with HCC based on treatment types. Survival differences were tested by the log-rank test and late-stage crossover accounted for using the two-stage method. (D) Univariate Cox regression analysis for OS probability plotted on a forest diagram. Analysis of PGK1 expression levels in tissue samples from patients with different (E) clinical and (F) pathological stages of HCC. Data were analyzed by the non-parametric Kruskal-Wallis with multiple comparisons Dunn's test. P-value was adjusted by the Bonferroni method. TCGA, The Genome Cancer Atlas; GTEx, Genotype-Tissue Expression; PGK1, phosphoglycerate kinase 1; HR, hazard ratio; OS, overall survival; HCC, hepatocellular carcinoma.

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.

Table I. Patient clinical characteristics according to phosphoglycerate kinase 1 expression.

Table I.

Patient clinical characteristics according to phosphoglycerate kinase 1 expression.

Characteristic	Low expression (n=178)	High expression (n=187)	P-value
Age, years			0.113a
<60	88	77
≥60	90	110
T stage			0.027a*
T1	101	79
T2	38	53
T3	32	46
T4	4	9
Unknown	3	0
N stage			0.364b
N0	127	121
N1	1	3
Unknown	50	63
M stage			1.000b
M0	133	130
M1	2	1
Unknown	43	56
Clinical stage			0.025b
I	98	72
II	37	47
III	33	50
IV	2	2
Unknown	8	16
Histological grade			0.312a
G1	27	28
G2	82	93
G3	57	61
G4	9	3
Unknown	3	2
Sex			0.183a
Male	114	132
Female	64	55
Family cancer history			0.184a
No	107	97
Yes	50	62
Unknown	21	28

{ label (or @symbol) needed for fn[@id='tfn1-ol-27-3-14242'] } Data analyzed using

a Pearson's χ² or

b Fisher's exact test. ‘Unknown’ groups were not included in the analysis.

Table II. Association between phosphoglycerate kinase 1 mRNA expression and clinicopathological stage of hepatocellular carcinoma for prognosis using a Kaplan-Meier plotter.

Table II.

Association between phosphoglycerate kinase 1 mRNA expression and clinicopathological stage of hepatocellular carcinoma for prognosis using a Kaplan-Meier plotter.

		Overall survival (n=364)		Recurrence-free survival (n=316)		Progression-free survival (n=370)		Disease-specific survival (n=362)
Characteristic	Patients (n)	HR	P-value	HR	P-value	HR	P-value	HR	P-value
Clinical stage
I	171	3.04	<0.001	2.23	0.004	2.11	0.003	3.53	0.005
		(1.62-5.68)		(1.28-3.87)		(1.27-3.52)		(1.39-8.96)
II	86	1.28	0.550	0.78	0.450	0.66	0.170	0.69	0.510
		(0.57-2.88)		(0.40-1.51)		(0.36-1.21)		(0.23-2.11)
III	85	0.56	0.065	0.54	0.088	0.52	0.040	0.57	0.140
		(0.30-1.05)		(0.27-1.11)		(0.28-0.98)		(0.27-1.21)
IV	5	-	-	-	-	-	-	-	-
T stage
T1	181	2.81	<0.001	2.17	0.004	2.05	0.004	2.95	0.008
		(1.55-5.11)		(1.27-3.72)		(1.25-3.38)		(1.27-6.83)
T2	94	1.43	0.370	0.83	0.550	0.67	0.160	1.82	0.250
		(0.65-3.12)		(0.44-1.55)		(0.38-1.18)		(0.64-5.17)
T3	80	0.60	0.120	0.65	0.240	0.59	0.067	0.58	0.170
		(0.32-1.15)		(0.31-1.34)		(0.34-1.04)		(0.26-1.27)
T4	13	-	-	-	-	-	-	-	-

[i] HR, hazard ratio.

Verification of increased PGK1 expression in 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.

Figure 6. Verification of increased expression of PGK1 in HCC tissue and liver cancer cell lines. (A) PGK1 expression in normal liver tissue and HCC specimens. Images were obtained from Human Protein Atlas online database (images.proteinatlas.org/73644/164848_B_8_5.jpg; http://images.proteinatlas.org/10065/25633_B_8_8.jpg). Compared with normal liver tissue, PGK1 was overexpressed in HCC cells. Scale bar, 35 µm. (B) Endogenous expression of PGK1 in six liver cancer cell lines was analyzed by western blotting. GAPDH was used as the internal control (P<0.01). HCC, hepatocellular carcinoma; PGK1, phosphoglycerate kinase 1.

Discussion

HCC is a complex disease caused by numerous risk factors. Therefore, it is difficult to diagnose HCC using the small number of currently available biomarkers. The discovery of biomarkers from different energy-producing mechanisms in HCC is key as it provides an improved understanding of these mechanisms, which may provide new strategies to treat HCC (35). To investigate whether glycolysis-associated proteins can serve as the potential molecular prognostic marker of HCC, a gene signature of HCC tissue was constructed and the mRNA profile obtained from the GSE46408 HCC dataset was assessed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. The key glycolytic enzyme involved in metabolic genes and pathways, PGK1, was selected for further analysis (36,37). Carbohydrate metabolism, which serves an important role in metabolism of energy and biosynthesis of macromolecules and occurs primarily in the liver, was the primary focus of this analysis (38,39). PGK1, which is one of the critical enzymes in the aerobic glycolysis process, was identified as serving an important role in a variety of biological functions to address the high metabolism requirements of HCC (40). In the present study, PGK1 was identified as a good predictive prognostic marker for the early stage of HCC. Moreover, PGK1 levels were positively associated with the progression of liver cancer, which was confirmed via immunohistochemical analysis of human tumor tissue and upregulated protein levels found in several tumor cell lines. PGK1 was identified as a potential marker for early HCC diagnosis and diagnostic features of PGK1 expression were identified that may potentially be used to distinguish between HCC and adjacent non-cancerous tissue.

The present study demonstrated that higher PGK1 expression was associated with a worse prognosis of patients with HCC. Consistent with the results of the present study, significant upregulation of PGK1 in HCC tissues has a negative association with survival of patients with HCC (41,42). Of note, high PGK1 expression was associated with the T stage of TNM stage progression in HCC and significantly different in early stage T1 compared with other stages. Nevertheless, there was no significant association between PGK1 expression and age, sex, histological grade, ethnicity and family cancer history. The upregulated PGK1 expression was demonstrated in HCC tissues and verified through immunohistochemical analyses and western blotting. A previous study identified that increased expression levels of PGK1 in the liver cancer cell lines SNU449 and HCCLM3 accelerate proliferation and metastasis through the activation of the Warburg effect (37). In the present study, upregulated expression of PGK1 in the liver cancer cell lines was observed, with the exception of the HuH7 cell line. The lowest protein expression levels of PGK1 in HuH7 consistent with previous studies (41,43), which suggested that HuH7 could potentially be used for PGK1 overexpression assays and the other five cell lines could be used for PGK1 knockdown in subsequent future experiments.

A number of studies on both cell and animal models have focused on the mechanism of PGK1 in liver cancer, but analysis of clinical data has not yet been studied, to the best of our knowledge (41,44–46). A previous study reported upregulation of PGK1 in some HCC cases, but did not compare this data with other HCC cases regarding PGK1 (47). Web-based genomic databases not only provide gene expression profiles, but also contain detailed clinical data from each patient. The correlation analysis between key genes and progression of disease is more systematic and scientific compared with traditional offline analytical approaches utilizing small localized genomic datasets (48). A previous study demonstrated that upregulation of PGK1 promotes progression of HCC; to the best of our knowledge, however, no strategy has been proposed for the early diagnosis of HCC using PGK1 as a biomarker or in combination with other biomarkers (49). Biomarkers are a key component of clinical management of cancer because they improve survival rate (50). Despite current efforts toward the discovery of novel prognostic or predictive biomarkers in solid tumors, <1% of these biomarkers are estimated to enter clinical practice because of the lack of external validation and clinical applicability of prognostic studies (51,52). Clinically, combination of AFP, AFP-L3 and DCP as non-invasive biomarkers is currently the most widely used combination, except in early diagnosis (53,54). In previous years, due to the key role of PGK1 in tumors, PGK1 has been investigated as a potential biomarker for tumor therapy (55,56).

There are limitations to the present study. First, cut-off values for the high and low PGK1 expression groups varied in different analyses. For example, the mean level was set as the cut-off value in the OS, Cox regression and χ2 or Fisher's exact test analyses, whereas the best-performing threshold value was used as the cut-off in the Kaplan-Meier plotter tool. The Kaplan-Meier plotter provides the best cut-off value by evaluating all possible cut-off points within the range of observed gene expression values in the dataset to avoid missing associations due to the use of a specific cut-off. This data-driven approach objectively identifies the cut-off with the strongest association with outcome rather than relying on arbitrary values. It helps maximize the power to detect true relationships and predict prognosis. The use of a single cut-off value in OS, Cox regression and χ2 or Fisher's exact test analyses rather than the best cut-off value is a limitation. Second, not all of the tissue samples contained complete clinical information, resulting in not all samples being included in all analyses. Third, stage information was determined according to the 7th edition of the AJCC, which is now outdated. Staging criteria from other countries and regions were also not considered. In future, a large number of population-based cohort studies are needed to determine whether PGK1 detection can improve early diagnosis of HCC by combining current effective methods for HCC diagnosis with PGK1 detection.

In the present study, survival and Kaplan-Meier curve analyses suggested that chemotherapy and targeted therapy provide a survival advantage over time compared with other types of treatment. The development of targeted drugs continues to expand cancer treatment options and improve patient survival. It also demonstrated the ability of PGK1 expression to predict prognosis, particularly in the early stages of disease. Univariate Cox regression analysis demonstrated that the risk score calculated by the model formula was an independent risk factor. In addition, PGK1 levels were positively associated with HCC progression, which was confirmed by immunohistochemistry analysis.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the General Program of the National Natural Science Foundation of China (grant nos. 82074305 and 32273036), Basic and Applied Basic Research Fund of Guangdong Province (grant nos. 2021A1515010822 and 2019A1515011994) and the Project of Administration of Traditional Chinese Medicine of Guangdong Province of China (grant no. 20182022).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

JH and HX were responsible for study conception and design. JY, XL, WH, WY and YQ were responsible for data analysis and interpretation. JY, JH and HX wrote the manuscript. JH and HX confirm the authenticity of all the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References