Development and validation of a prognostic classifier based on HIF-1 signaling for hepatocellular carcinoma

Expression profile of HIF-1 related genes in HCC

Sixteen HIF-1 related genes were identified in previous study and KEGG pathways, including Solute Carrier Family 2 Member 1 (SLC2A1), Solute Carrier Family 2 Member 3 (SLC2A3), Hexokinase 1 (HK1), Hexokinase 1 (HK2), Aldolase, Fructose-Bisphosphate A (ALDOA), Aldolase, Fructose-Bisphosphate C (ALDOC), Phosphoglycerate Kinase 1 (PGK1), Enolase 1 (ENO1), Pyruvate Kinase M1/2 (PKM), Lactate Dehydrogenase A (LDHA), Pyruvate Dehydrogenase Kinase 1 (PDK1), BCL2 Interacting Protein 3 (BNIP3), Phosphofructokinase, Liver Type (PFKL), Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH), 6-Phosphofructo-2-Kinase/Fructose-2,6-Biphosphatase 3 (PFKFB3), Hypoxia Inducible Factor 1 Subunit Alpha(HIF1A) [17]. The study design was displayed in Figure 1. 365 patients with HCC in The Cancer Genome Atlas (TCGA) database were enrolled as the derivation cohort. HIF-1 related genes were significantly upregulated in tumors compared with adjacent normal tissues (Figure 2A, 2B). Moreover, correlation analysis validated the association among these genes in HIF-1 related signaling (Figure 2C).

Prognostic relevance of HIF-1 related genes in HCC

Kaplan-Meier analysis and log-rank test were adopted to evaluate the prognostic value of sixteen HIF-1 related genes in HCC (Figure 3). SLC2A1, ENO1, LDHA, ALDOA, GAPDH, HK2, PKM, PGK1, HIF1A, PFKFB3, and PDK1 forecasted poor overall survival (P<0.05). The univariate Cox regression model was used to further identify the prognostic relevant genes for overall survival in patients with HCC. Eleven genes, including SLC2A1, ENO1, LDHA, ALDOA, GAPDH, HK2, PKM, PGK1, HIF1A, PFKFB3, PDK1 were significantly correlated with clinical prognosis (P<0.05) (Figure 4A).

Development of a HIF-1 related gene-based risk model and prognostic score

Furthermore, the identified prognostic features were included in the LASSO Cox regression model for prognostic model and risk score construction (Figure 4B, 4C). Four genes (SLC2A1, ENO1, LDHA, and GAPDH) were selected to establish a risk scoring system due to their integrated prognostic relevance. As SLC2A1, ENO1, LDHA, and GAPDH were recognized as genes which promote anaerobic metabolism in human cells, we further developed an oxygen-sensing specific risk score for patients with HCC. The coefficients for corresponding features were generated according to the partial likelihood deviance and determined with its lowest value at a log λ=-3.45. Four features were included in a formula based on gene expression and the risk score were calculated as follows:

$\text{risk}\text{score}={\text{e}}^{0.111\ast \text{SLC}2\text{A}1+0.186\ast \text{ENO}1+0.278\ast \text{LDHA}+0.042\ast \text{GAPDH}}$

The patients were further assigned to the low-risk and high-risk groups by median (risk score=23). PCA and t-SNE analysis confirmed the clustering ability of this four-gene based risk score (Figure 4D, 4E). Functional analysis showed different expression profile between the low-risk and high-risk group (Supplementary Figure 1). GO and KEGG analysis indicated that HIF-1 and hypoxia-related metabolic signaling were enriched in high-risk group (Supplementary Figure 2). Furthermore, PPI network analysis revealed the interactions of these differently expressed genes (Supplementary Figure 3).

Prognostic value of four-gene based risk score in patients with HCC

Based on the scoring system, the prognostic risk of each patient was individualized. The time-dependent receiver operating characteristic (ROC) analysis showed that the area under the curve (AUC) value of four-gene based model was 0.753 (95% CI, 0.673-0.834), 0.710 (95% CI, 0.619-0.800) and 0.669 (95% CI, 0.557-0.781) at 1 year, 3 years and 5 years after diagnosis in the TCGA cohort, respectively (Figure 5A–5C). Moreover, the risk score presented with a better predictive ability compared with the TNM stage during the short-term follow-up (Supplementary Figure 4). The risk score seemed to be more precise in the short-term follow-up of male patients, while more precise in 5-year follow-up of female patients (Supplementary Figure 5). Besides, patients with continuous risk scores harbored various clinical outcomes in different groups (Figure 5D, 5E).

Kaplan-Meier plot and log-rank test showed that patients at high-risk exhibited a significantly poorer overall survival compared to those at low-risk (Figure 6A). Restricted cubic spline (RCS) revealed that accompanied by the score increases, the risk raised. And the four-gene based risk scoring system was equipped with a continuous predictive ability (Figure 6B). Univariate Cox hazard model in patients with HCC indicated that TNM stage (HR: 1.369, 95% CI: 1.204−1.557, P<0.001), T stage (HR: 1.572, 95% CI: 1.331−1.857, P<0.001), N stage (HR: 1.175, 95% CI: 0.973−1.419, P=0.094), M stage (HR: 1.237, 95% CI: 1.024−1.495, P=0.028) and risk score (HR: 2.220, 95% CI: 1.546−3.187, P<0.001) could predict the overall survival of these patients (Figure 6C). Multivariate Cox hazard model revealed that T stage (HR: 1.527, 95% CI: 1.195−1.951, P<0.001), M stage (HR: 1.344, 95% CI: 1.039−1.738, P=0.024) and the risk score (HR: 2.238, 95% CI: 1.535−3.263, P<0.001) could be independent predictors for the above patients (Figure 6D). Furthermore, Chi-square analysis indicated that the risk score was significantly correlated with TNM stage and T stage in patients with HCC (Figure 6E).

Infiltrating immune cells profile in different groups

To further identify the difference in tumor-infiltrating immune cells profile between patients with low risk and high risk, we calculated the immune cells fraction using CIBERSORT algorithm (Figure 7A). Patients with high risk exhibited a lower level of CD8+ T cells (P=0.072), activated NK cells (P=0.004), macrophages M1 (P=0.011) and resting mast cells (P<0.002). On the other hand, patients with high risk showed a higher level of regulatory T cells (P=0.077), macrophages M0 (P<0.001) and neutrophils (P=0.002). Moreover, an ESTIMATE algorithm was used for tumor microenvironment scoring, by which the immune score and stromal score of HCC were estimated. Patients with high risk showed a significantly higher immune score compared with those at low risk, but there was no significant difference in terms of stromal score (Figure 7B, 7C). These results further confirmed that the four-gene based risk score could serve as a reliable predictor for clinical outcomes and might be an indicator of immune cells infiltration profile in HCC.

Further validation of the clinical relevance of risk score in GEO dataset

To further validate the risk score system, a GEO dataset, GSE14520, was included in our study. A total of 445 HCC samples (GPL3921 platform) were analyzed. Results suggested that there was a significant correlation among HIF-1 related genes (Figure 8A). Among the four selected genes, ENO1, GAPDH, and SLC2A1 were significantly upregulated in tumor tissues compared with normal tissues (Figure 8B). The stratification ability of the risk score system was also confirmed by PCA and t-SNE analysis in tumor tissues (N=225) (Figure 8C, 8D). Time-dependent ROC analysis showed an AUC value of 0.671 (95% CI: 0.571-0.772) at 5 years (Figure 8E).

Survival analysis was performed to confirm the clinical relevance of the risk score system. Patients in the high-risk group exhibited a higher rate of mortal events (Figure 9A). Kaplan-Meier analysis and log-rank test showed a poorer overall survival in the high-risk group (P<0.01) (Figure 9B). Univariate Cox hazard model indicated that TNM stage (HR: 2.088, 95% CI: 1.619−2.693, P<0.001) and risk score (HR: 1.958, 95% CI: 1.265−3.031, P=0.003) were reliable for survival prediction (Figure 9C). Moreover, the multivariate Cox hazard model confirmed that TNM stage (HR: 1.958, 95% CI: 1.509−2.542, P<0.001) and the risk score (HR: 1.604, 95% CI: 1.023−2.514, P=0.039) were independent prognostic indicators (Figure 9D).

Study subjects

For TCGA (The Cancer Genome Atlas) and GEO (Gene Expression Omnibus) datasets, the inclusion criteria were as follow: 1) available for OS data and mRNA expression; 2) histologically diagnosed as malignant hepatocellular carcinoma; 3) available for active follow-up data. Patients without detailed or active follow-up were excluded. Due to the unknown identity of patients from the above two datasets, written informed consents were waived by the Ethics Committee of The Affiliated Foshan Hospital, Sun Yat-sen University.

The TCGA dataset was obtained from TCGA using gdc-client. The mRNA expression matrices of GSE14520 were downloaded from the GEO dataset [28]. 424 samples from TCGA and 445 samples from GSE14520 were finally enrolled in our study as the derivation and validation cohort, respectively.

Development of risk score based on HIF-1 related genes

The data of TCGA was downloaded using gdc-client and the expression data was further normalized using edgeR package in R software (Version 3.5.3). The data of GSE14520 was downloaded from the GEO dataset and the limma package was applied to normalization. Univariate Cox hazard model was adopted to identify prognostic relevant signatures from HIF-1 related genes. Eleven signatures were subjected to the LASSO Cox regression model using the glmnet package [29]. Risk scores were calculated using the generated coefficients and corresponding expression. The patients were further divided into low-risk and high-risk groups according to risk scores.

Assessment of risk score system

To further confirm the relevance of the aforementioned system, the clustering ability of risk scores was further examined by PCA and t-SNE analysis. Time-dependent ROC analysis was performed and AUC at different time points were calculated to assess the diagnostic value of risk scores [30]. Survival analyses, including Kaplan-Meier plots, log-rank test, univariate and multivariate Cox hazard model were also adopted. GO and KEGG analyses were conducted using clusterProfiler package in R [31]. GO and PPI networks were constructed by Metascape [32].

Evaluation of infiltrating immune cells, immune score and stromal score

The tumor-infiltrating immune cells fraction was evaluated by uploading the expression matrices to the CIBERSORT, which calculated the fractions according to the LM22 signature with 1000 permutations [33]. To further quantify the tumor microenvironment of patients with HCC, the immune score and stromal score were calculated using the estimate package in R software (Version 3.5.3) [34].

Statistical analysis

The Chi-square test and the Mann-Whitney U test were used to compare the difference between subgroups. Student's t-test was performed to compare the expression level in different groups when appropriate. Spearman analysis was used for evaluating the correlation between gene expressions. Kaplan-Meier method and Cox regression model were conducted for survival analysis. The Lasso Cox regression model was used to determine the key relevant prognostic phenotypes and the corresponding ideal coefficient. The partial likelihood deviance was calculated by the glmnet package in R software for ideal coefficient selection. Statistical analysis was performed on SAS 9.4 (SAS Institute, Inc), Graphpad Prism 7 (Graphpad) and R 3.5.3 software (https://www.r-project.org). A P value <0.05 was determined as statistically significant and all P values were two-tailed.