A glycolysis-related gene pairs signature predicts prognosis in patients with hepatocellular carcinoma

Data source

The expression profiles and clinicopathological data of HCC and normal tissues were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and International Cancer Genome Consortium Japan (ICGC, https://dcc.icgc.org/releases/current/Projects/LIRI-JP). The clinical information inclusion criteria were set as follows: (1) patients had completely detailed clinical information; (2) The follow-up time of samples exceeded 30 days. GSE10140, GSE10141, GSE10143, GSE15654, GSE14520, GSE76427, GSE45114 expression profile was derived from the Gene Expression Omnibus (GEO) database, including 1,469 samples. The TCGA was used as training dataset. Other databases were used as validation datasets.

Gene set enrichment analysis

To explore whether the specific gene sets were significant different between tumor group and normal group, we performed gene sets enrichment analysis (GSEA) (http://www.broadinstitute.org/gsea/index.jsp). The mRNAs expression levels between tumor and non-tumor groups were analyzed to confirm which function could be available for further study. Normalized P-value ≤ 0.05 were considered to be statistically significant.

Construction and evaluation of glycolysis‑related gene pairs signature

First of all, the glycolysis‑related mRNAs level in the same sample was pairwise compared to generate a score for each glycolysis‑related gene pair (GRGP). If the expression level of gene 1 was greater than gene 2, the GRGP score was 1, otherwise it was 0 (Li et al., 2017). GRGPs with high variation in TCGA dataset (median absolute deviation >0.05) were included in subsequent model construction. The prognostic value of GRGPs was identified by univariate Cox regression. Then, GRGPs with P ≤ 0.05 in univariate analysis were incorporated into Lasso regression model in order to establish a GRGPs signature. A risk score was established according to the following formula: risk score = ${\sum }_i\mathrm{C}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\left(\mathrm{G}\mathrm{R}\mathrm{G}{\mathrm{P}}_i\right)\ast \mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(\mathrm{G}\mathrm{R}\mathrm{G}{\mathrm{P}}_i\right)$. To classify patients into low-risk group and high-risk group, the optimal cut off of the GRGPs signature was defined using time-dependent receiver operating characteristic (ROC) curve analysis at 1 year in the TCGA dataset for overall survival (OS). The OS between two groups was compared utilizing Kaplan–Meier and Log-rank test. Risk score and other clinicopathological characteristics were included in the model so as to confirm whether risk score was an independent factor to predict the progress of the patients. Further, the clinical value of the GRGPs signature was verified by comparing the risk scores of patients with different ages, gender, grade and stage.

The OS of the patients with HCC at 1 year, 3 years and 5 years were predicted using a nomogram. Index of concordance (C-index) and Calibration curves were applied to explore the accuracy of the model established from TCGA dataset.

Validation of the GRGPs signature

Seven GEO databases (GSE10140, GSE10141, GSE10143, GSE15654, GSE14520, GSE76427, and GSE45114) and ICGC database were enrolled for subsequent verification. According to human race, all the databases were divided into European and American dataset and Asian dataset. European and American dataset (765 samples) included GSE10140, GSE10141, GSE10143, and GSE15654 database; Asian dataset (947 samples) included GSE14520, GSE76427, GSE45114 and ICGC database. Due to the lack of detailed clinical information (such as age, gender, stage) in European and American dataset, thus it was only used for survival analysis. Asian dataset was used for subsequent survival analysis and independent prognostic analysis. The GRGPs signature was calculated, and then Kaplan–Meier, Log-rank test and Cox regression were used to verify whether the GRGPs signature was significantly associated with OS. The ROC curve, C-index and Calibration curves were constructed to examine whether the model established by TCGA dataset could effectively predict survival in patients with HCC.

Statistical analysis

Cox regression was utilized to evaluate the associations between the expression levels of the glycolytic-related mRNAs and OS. Moreover, univariate and multivariate Cox regression were applied to determine the prognostic values of the GRGPs signature and various clinicopathological characteristics. The prediction accuracy of the risk score for 1-year, 3-years and 5-years survival was evaluated using ROC curve analysis. Statistical tests were two-sided, and P values below 0.05 were considered to be statistically significant.

Preliminary selection of genes using GSEA

According to inclusion criteria, 349 patients with HCC were obtained. Expression signatures of marker gene sets were obtained by condensing multiple gene sets from the Molecular Signature Database (MSigDB). Each expression signature involved 50 specific gene sets that stand for clearly defined biological statuses and processes (Zhang, Zhang & Yu, 2019). GSEA was applied to investigate whether the identified gene sets revealed significant differences between the tumor and normal groups. Forty-four gene sets were upregulated in HCC, and 20 gene sets were highly enriched, with normalized P < 0.05 among the 50 gene sets (Fig. 1; Table S1). As can be seen in the Table S1, the biological process of glycolysis (P < 0.05) contained 199 genes, which was the second largest in size.

Identification of prognostic GRGPs

Based on 199 GRGs, 19,701 GRGPs were established. After removing GRGPs with comparatively low variation (MAD > 0.05), 1,102 GRGPs were reserved and selected as initial candidate GRGPs. Cox regression was utilized to select prognostic GRGPs. According to the results of univariate Cox, 35 GRGPs had prognostic values for the patients with HCC (P < 0.05, Table S2). Subsequently, lasso regression model revealed that 17 GRGPs as prognostic factors (Table 1; Figs. 2A–2B). As shown in Table 1, five GRGPs (HMMR_PFKFB1, CHST1_GYS2, MERTK_GYS2, GPC1_GYS2 and LDHA_GOT2) were found to be harmful prognostic factors and twelve GRGPs (IDUA_GNPDA1, IDUA_ME2, IDUA_G6PD, IDUA_GPC1, MPI_GPC1, SDC2_LDHA, PRPS1_PLOD2, GALK1_IER3, MET_PLOD2, GUSB_IGFBP3, IL13RA1_IGFBP3 and CYB5A_IGFBP3) were verified to be beneficial prognostic factors.

Construction of GRGPs signature

These seventeen GRGPs were utilized to establish a GRGPs signature. Risk score of GRGPs signature for each patient was calculated utilizing the following formula (formula 1): risk score = (−0.30794 * IDUA_GNPDA1) − (0.15299 * IDUA_ME2) − (0.16389 * IDUA_G6PD) − (0.35599 * IDUA_GPC1) + (0.04846 * HMMR_PFKFB1) − (0.35632 * MPI_GPC1) − (0.29752 * SDC2_LDHA) − (0.09077 * PRPS1_PLOD2) − (0.06137 * GALK1_IER3) + (0.02511 * CHST1_GYS2) − (0.26287 * MET_PLOD2) − (0.00305 * GUSB_IGFBP3) + (0.34302 * MERTK_GYS2) + (0.20608 * GPC1_GYS2) − (0.31484 * IL13RA1_IGFBP3) + (0.17629 * LDHA_GOT2) − (0.10962 * CYB5A_IGFBP3). The cutoff point of risk score was set at −0.698 utilizing ROC curve analysis, which classified the patients into high-risk group or low-risk group (Fig. 2C). Risk score was significantly associated with OS of the patients with HCC, in which OS in low-risk group was longer than that in high-risk group (P = 2.898e−14, Fig. 3A). The survival time of patients with HCC decreased with risk score increasing (Fig. 4).

Association between risk score and clinicopathologic factors

Risk score increased with age, stage, and survival status, demonstrating that the GRGPs signature might be relevant to the progression of HCC. Risk score of patients with advanced-stage and advanced-age were significantly higher than those with early-stage and early-age (Figs. S1A and S1B). Risk scores of dead patients were higher than those of living patients (Fig. S1C). The heat map illustrated that the high expression of these seventeen GRGPs were significantly related to female, lower survival status of patients, higher stage and higher grade (Fig. S1D).

The GRGPs signature was an independent prognostic factor

In TCGA dataset, univariate Cox regression revealed that GRGPs signature was associated with OS and its HR was 3.508 (95% CI [2.608–4.720], P < 0.001, Table S3; Fig. 5A). After controlling clinical features such as gender, age, tumor stage and grade, GRGPs signature remained to be an independent prognostic indicator (HR = 3.204, 95% CI [2.293–4.476], P < 0.001, Table 2; Fig. 5B). GRGPs signature and TNM stage were independent prognostic factors based on the TCGA. Thus, these factors were included in nomogram. GRGPs signature was the largest contributor to 1-year, 3-year and 5-year OS (Fig. 6A). The C-index of the nomogram was 0.764 (95% CI [0.715–0.813]). The areas under the ROC curve (AUC) corresponding to the survival at 1 year, 3 years and 5 years were 0.803, 0.777 and 0.774, respectively (P < 0.05). The C-index, ROC curve and Calibration explained that the GRGPs signature had better accuracy (Figs. 6B–6D).

Validation of the GRGPs signature

In validation datasets, the risk score of GRGPs signature was calculated according to formula 1. The risk score was also significantly correlated with OS of patients with HCC (European and American dataset: P = 1.143e−02, Fig. 3B; Asian dataset: P = 6.342e−08, Fig. 3C). Univariate independent prognostic analysis indicated that GRGPs signature were independent prognostic factors in Asian dataset (HR of risk score = 2.661, 95% CI [1.862–3.803], P < 0.001, Table S4; Fig. 5C). After controlling stage and gender, GRGPs signature remained an independent prognostic indicator in multivariate analysis (HR = 2.567, 95% CI [1.714–3.844], P < 0.001, Table S5; Fig. 5D). According to results of independent prognostic analysis, GRGPs signature, stage and gender were included in nomogram based on Asian dataset. GRGPs signature and age were the largest contributor to 1-year, 3-year and 5-year OS in Asian dataset nomogram (Fig. 6E). The C-index of nomogram based on Asian dataset was 0.705 (95% CI [0.660–0.750]). The areas under the ROC curve at 1 year, 3 years and 5 years in Asian dataset were 0.694, 0.664, 0.536, respectively (Fig. 6H). The 3-year and 5-year calibration curves also proved that the GRGPs signature had great accuracy and robustness (Figs. 6F and 6G).