A Novel and Robust Prognostic Model for Hepatocellular Carcinoma Based on Enhancer RNAs-Regulated Genes

Evidence has demonstrated that enhancer RNAs (eRNAs) play a vital role in the progression and prognosis of cancers, but few studies have focused on the prognostic ability of eRNA-regulated genes (eRGs) for hepatocellular carcinoma (HCC). Using gene expression profiles of HCC patients from the TCGA-LIHC and eRNA expression profiles from the enhancer RNA in cancers (eRic) data portal, we developed a novel and robust prognostic signature composed of 10 eRGs based on Lasso-penalized Cox regression analysis. According to the signature, HCC patients were stratified into high- and low-risk groups, which have been shown to have significant differences in tumor immune microenvironment, immune checkpoints, HLA-related genes, DNA damage repair-related genes, Gene-set variation analysis (GSVA), and the lower half-maximal inhibitory concentration (IC50) of Sorafenib. The prognostic nomogram combining the signature, age, and TNM stage had good predictive ability in the training set (TCGA-LIHC) with the concordance index (C-index) of 0.73 and the AUCs for 1-, 3-, and 5-year OS of 0.82, 0.77, 0.74, respectively. In external validation set (GSE14520), the nomogram also performed well with the C-index of 0.71 and the AUCs for 1-, 3-, and 5-year OS of 0.74, 0.77, 0.74, respectively. In addition, an important eRG (AKR1C3) was validated using two HCC cell lines (Huh7 and MHCC-LM3) in vitro, and the results demonstrated the overexpression of AKR1C3 is related to cell proliferation, migration, and invasion in HCC. Altogether, our eRGs signature and nomogram can predict prognosis accurately and conveniently, facilitate individualized treatment, and improve prognosis for HCC patients.

1 Introduction

Hepatocellular carcinoma (HCC) is the most common type (approximately 90%) of primary liver cancer, which has already become the sixth most commonly diagnosed cancer (906,000 new cases) and the third leading cause of cancer death (830,000 deaths) worldwide in 2020 (1). Chronic infection of hepatitis B virus (HBV) or hepatitis C virus (HCV), cirrhosis, excessive alcohol consumption, and type 2 diabetes are the main risk factors for HCC (2). HCC has a very poor prognosis due to its advanced stage, rapid progression, high recurrence rate, and limited treatment options (3). Traditionally, tumor stage is a widely used basis for predicting the prognosis of patients with HCC (4, 5). However, prognosis in HCC is complex and highly heterogeneous (6, 7), patients in the same tumor stage may present significantly different prognosis. A valuable, accurate strategy is still an urgent need to predict HCC prognosis.

Previous research has proposed that serum biochemical biomarkers, such as hepatic growth factor, osteopontin, BALAD scoring model composed of alpha-fetoprotein (AFP), alpha-fetoprotein lens culinaris agglutin-3 (AFP-L3), and Des-γ-carboxy prothrombin (DCP) could be used for predicting prognosis (8–10). With the development of high-throughput sequencing technology, genes were considered to be important factors in predicting survival of patients with HCC (11–14). Furthermore, gene signatures based on multigene expression, and nomograms including gene signature and clinical information have been constructed for patients’ prognosis in many published studies (15–17). However, it is still necessary to further mine omics data in combination with clinical characteristics to discover a novel and reliable prognostic model for HCC and guide patients’ individualized treatment.

Enhancer RNAs (eRNAs) are non-coding RNAs transcribed by enhancers that mediate the activation of target genes (18, 19). In human cancers, eRNAs can contribute to the activation of oncogenes or oncogenic signaling pathways and can be induced by oncogenes or tumor suppressors to directly participate in tumor promotion or inhibition process (20–22). Furthermore, eRNAs can also bind to DNA, proteins (e.g., transcription factors, cofactors, RNA-binding proteins) to regulate their function and activity (23). Super-enhancers, the clusters of enhancers, also play prominent roles in dysregulation of oncogenes expression, tumor suppressor genes expression, process of tumorigenesis, and proliferation of cancer cells (23–25). Although some data portals have already annotated a mass of enhancers, the direct interaction between eRNA and its target genes has not been elaborated until researchers integrated data from TCGA and other projects (26–28) and developed eRic data portal including a global eRNA-gene regulatory network across 31 cancer types in 2019 (29). This portal facilitates a deeper investigation of relationships between eRNAs and cancers for biomedical researchers. However, few studies have focused on eRNA-regulated genes (eRGs) to construct prognostic models of HCC patients for the management of individualized treatment.

This study attempted to construct a novel prognostic gene signature composed of eRGs and nomogram combining the signature with clinical characteristics by using univariate and Lasso-penalized Cox regression analyses in TCGA-LIHC dataset and validate it in GSE14520 dataset. We demonstrated the application potential of the signature and nomogram through bioinformatics methods. In addition, we validated the effect of important eRG in the signature on cell proliferation, migration, and invasion of HCC using two HCC cell lines (Huh7 and MHCC-LM3) in vitro.

2 Materials and Methods

2.1 Data Collection and Processing

Data used in this study were all publicly available. mRNA expression and clinical data of 374 LIHC patients were obtained from the Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/), and the eRNA expression profile across these TCGA samples as well as eRNA target genes were acquired from the enhancer RNA in cancers (eRic) data portal (https://hanlab.uth.edu/eRic/). The raw mRNA data in TCGA was processed to fragments per kilobase of transcript per million mapped reads (FPKM) and transformed based on log2 (log2FPKM). Patients with insufficient survival information or follow-up period less than 30 days were excluded, and 324 HCC samples were selected as the training set for subsequent analysis. Microarray dataset GSE14520 -GPL3921 includes 219 HCC patients with integral clinical information and survival time longer than 30 days were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) as the external validation set.

2.2 Construction and Validation of eRGs Signature

In this section, we identified eRNAs related to survival in HCC patients and the target genes regulated by them. Based on the target genes, an eRGs signature was constructed and validated for predicting prognosis.

Kaplan-Meier analysis and univariate Cox regression were used to screen eRNAs related to survival in TCGA dataset, and only the eRNAs with p < 0.05 in both above analyses were selected. Then the target genes of these eRNAs were identified by referring to the eRic database.

The target genes of survival-related eRNAs in TCGA-LIHC were firstly subjected to univariate Cox regression analysis and genes with p < 0.05 were considered as the candidate prognostic eRGs. Subsequently, Lasso-penalized Cox regression analysis was performed to build the prognostic eRGs signature, and 10-fold cross-validation was carried out to determine the optimal penalty parameter using “glmnet” package in R. Then the prognostic eRGs signature was established, and the risk score of each patient can be calculated based on the corresponding coefficients of genes from the Lasso Cox regression model (β) and their expression level:

$$riskScore=\sum _{i=1}^m{\beta }_i\times Ex{p}_i$$

where m is the number of the genes in signature, β_i is the coefficient from Lasso Cox analysis, Exp_i is the expression level of the eRGs in signature. The median risk score was used as the cutoff value to stratify the HCC patients into high risk and low risk groups. Thereafter, we performed a log-rank test to compare the survival rates between the two groups and plotted Kaplan-Meier survival curves using “survival” and “survminer” package in R. Time-dependent receiver operating characteristic (ROC) curves for 1-, 3-, and 5-year OS were also drawn based on the risk score using “timeROC” package in R to assess the prognostic performance of the signature. To validate the predictive capability and generalization of the eRGs signature, HCC patients in GSE14520 dataset with intact survival and clinical information were considered as the external validation set. Risk scores of these patients were calculated with the same model of prognostic eRGs signature, and Kaplan-Meier survival curves analysis and the ROC curves analysis were also conducted in this dataset.

2.3 Construction and Validation of Prognostic Nomogram

We performed univariate Cox regression and multivariable Cox regression analysis on the risk score of the gene signature and other clinical features (including age, TNM stage, AFP level, and BMI) in TCGA dataset to identify the independent prognostic factors, and factors with p < 0.05 were deemed statistically significant. The proportional hazard assumption of the model was tested by Schoenfeld residuals test. Then a nomogram based on these independent prognostic factors was developed using “rms” package in R to predict the overall survival time of HCC patients. Nomogram is a graphical representation of a complex mathematical formula (30). In medicine, a nomogram is usually used to graphically describe a statistical prognostic model that generates the probability of a particular individual’s clinical events (such as cancer recurrence or overall survival), which can be conductive to personalized medicine (31). Calibration curves were utilized to investigate the consistency between the nomogram-predicted probabilities and the actual survival rates, and decision curves analysis (DCA) was conducted to assess the clinical predictive value of the nomogram. Furthermore, time-ROC curves were also plotted in TCGA dataset to evaluate the prognostic performance of the nomogram. Based on the nomogram in the training set, calibration curves analysis, DCA analysis and time-ROC curves analysis were all performed in the external validation set GSE14520 to further validate the results.

2.4 Enrichment Analyses and Gene-Set Variation Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed with survival-related eRGs set, and terms with adjusted p < 0.05 were considered as significant enrichment. Gene Set Enrichment Analysis (GSEA) was also conducted on the genes in the survival-related signature to explore the potential pathways by using “clusterProfiler” package in R (32). After that, we applied gene-set variation analysis (GSVA) using “GSVA” package and “limma” package in R to identify the different pathways between the two different risk groups following the criteria of |log2FC| > 0.2 and p < 0.05 (33). The annotated gene set used in GSEA and GSVA was “c5.all.v6.2.symbols.gmt”, which can be downloaded from the Molecular Signatures Database (MSigDB).

2.5 Exploration of Immunotherapy Effect and Immune Landscape

Immune checkpoint inhibitors play an important role in the treatment of liver cancer (34), so we compared the differences in the expressions of some common immune checkpoints (CTLA4, PD-L1, TIGIT, and HAVCR2) between high risk patients and low risk patients using Wilcoxon rank-sum test. We also calculated the proportion of 22 immune-infiltrating cells in each patient using “Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT)” algorithm (35), and the patients with p < 0.05 were used for difference comparison in the two sub-groups. In addition, we utilized TIMER (Tumor Immune Estimation Resource) database to further explore the correlation between risk score and six immune cells including B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells (36). Human leukocyte antigen (HLA) is the major histocompatibility complex (MHC) in human, which goes together with the function of human immune system. Therefore, we also validated the differences in HLA-related genes between the two groups.

2.6 Evaluation of Drug Sensitivity

Half maximal inhibitory concentration (IC50) is a widely used indicator to evaluate the sensitivity of drug therapy (37). Sorafenib, a targeted therapy drug, has been considered the standard treatment for patients with hepatocellular carcinoma (HCC) since 2007 (38). Therefore, we estimated the IC50 of Sorafenib in high- and low-risk patients and determined whether there is a difference in the sensitivity of different patients to sorafenib. The drug IC50 was estimated using “pRRophetic” package in R based on Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) cell line expression data and TCGA-LIHC gene expression data (39, 40).

2.7 Identification of Transcription Factors Related to the Genes in Signature

Transcription factor (TF) list was downloaded from Cistrome (http://cistrome.org/), and the expression level of these TF were extracted from 324 TCGA-LIHC patients. Then we performed Spearman correlation test between the expression of eRGs in the signature and these TFs. The gene-TF pairs with absolute Spearman correlation coefficients > 0.4 and p < 0.05 were selected to further discussion.

2.8 In Vitro Experimental Validation of AKR1C3 in HCC Cell Lines

2.8.1 Cell Culture, Transfection and Sorafenib Treatment

Human hepatocellular carcinoma cell lines Huh7 and MHCC-LM3 were obtained from the Cell Bank of the Chinese Academy of Sciences (Beijing, China), both types of cells were cultured in DMEM-6429 (Sigma, MO, USA) with 10% fetal bovine serum (FBS, HyClone, Logan, UT, USA). The siRNA sequences targeting AKR1C3 used in this study were as follows: si-AKR1C3-1:(5’-CCAAACACCAGUGUGUAAATT-3’, 5’-UUUACACACUGGUGUUUGGTT-3’); and si-AKR1C3-2:(5’-GGAACUUUCACCAACAGAUTT-3’, 5’-AUCUGUUGGUGAAAGUUCCTT-3’) and negative controls (si-AKR1C3-NC); AKR1C3 overexpression plasmid:(5’-ATGGATTCCAAACACCAGTGT-3’, 5’-TTAATATTCATCTGAATATGG-3’) and empty plasmids were regarded as negative controls (NC) purchased from the Nantong Biomics Biotechnologies company. Transfection was performed using Lipofectamine 3000 (Thermo Fisher Scientific, Inc) according to the manufacturers’ instructions. The cells were harvested 48 hours after transfection. In addition, Sorafenib (8μM) was added to Huh7 and MHCC-LM3 cells prior the incubation.

2.8.2 Protein Extraction and Western Blot Analysis

Total protein was extracted using RIPA buffer and quantified with a BCA kit (Beyotime Biotechnology). Protein separation was performed using 10% SDS-PAGE and then transferred to PVDF membranes (Merck Millipore). Following being blocked with 5% non-fat milk for 2 h at 25°C, the PVDF membranes were incubated with AKR1C3(Abcam; 1:1000; ab209899) and GAPDH (Abcam; 1:5000; ab9485), and then incubated with a secondary antibody (Cell Signaling Technology) for 1 h. PVDF membranes were scanned by a chemiluminescence system.

2.8.3 Cell Viability and Colony Formation Assays

EdU cell proliferation assay was performed using a commercial EdU Kit (UE, China) according to the manufacturer’s protocol. Images obtained from a fluorescence microscope (Leika, Germany) were analyzed using Image J. The colony formation assay was used to evaluate the cell clonogenic ability. The transfected Huh7 and MHCC-LM3 cells were seeded in a 35mm-diameter petri dish and cultured for up to 14 days, respectively. Cell colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Beyotime) for 20 minutes, the colonies were counted under a light microscope.

2.8.4 Transwell Assay

Cell invasion was evaluated by performing the Chamber matrigel invasion 24-well units (Costar) according to the manufacturer’s instructions. The transfected cells were suspended in a serum-free medium and plated into the upper chamber of the transwell system with a pore size of 8 µm. The bottom chamber was filled with a medium containing 10% FBS. After incubation for 24 h, the migrated/invaded cells in the lower chamber (below the filter surface) were fixed in 4% paraformaldehyde, stained with crystal violet solution, and counted under a microscope.

2.8.5 Wound Scratch Assay

Wound scratch assays were used to assess the migratory ability of Huh7 and MHCC-LM3 cells in vitro. AKR1C3 downregulated Huh7 and MHCC-LM3 cells (including negative control cells) were planted in 3.5 cm dishes and grown until 80%–90% confluent. Then, a 100 μl yellow pipette tip was used to scratch the cell monolayers and the cells were maintained in DMEM-6429 medium. The area of the cell-free wound was measured with microscopy at 0 and 24 h.

2.8.6 CCK-8 Experiments

The CCK-8 assay (Dojindo, Japan) was performed to assess Huh7 and MHCC-LM3 cells proliferation. Cells were seeded at a density of 4×10³ cells/well in 96 wells plates, then added 20 µl of CCK-8 reagent to each well of a 96-well plate, and incubated the cells for 2 h at 37°C. At 6, 24, 48, 72, and 96 h, cell viability was detected by scanning with a microplate reader (Tecan, Switzerland) at 450 nm.

2.9 Statistical Analysis

R software 3.6.3 was used for all data management and analyses in present study. Wilcoxon rank-sum test was used to compare the differences of quantitative variables, and Spearman correlation test was used to explore the correlation between variables. Schoenfeld residuals test was performed to test the proportional hazard assumption of Cox regression model. All the statistical tests were two-sided, and p < 0.05 were considered to be statistically significant.

3 Results

3.1 Survival-Related eRNAs, eRGs and Prognostic Signature

The overall flow chart of our study is depicted in Figure 1. A total of 324 HCC patients from TCGA-LIHC were included as the training set, and 219 HCC patients from GSE14520-GPL3921 were used as the external validation set. The general clinical characteristics of the two datasets are exhibited in Supplementary Table S1.

FIGURE 1

Figure 1 Flowchart elaborating the scheme of construction and validation of prognostic signature and nomogram based on eRGs for HCC patients.

eRNA expression profile data of 324 TCGA-LIHC patients were downloaded from the eRic data portal, and 457 eRNAs were obtained for K-M analysis and univariate Cox regression. Based on the results of survival analysis, 46 eRNAs were identified that were significantly associated with overall survival (p < 0.05). Then 95 target genes of the survival-related eRNAs were obtained by referring to the eRNA target genes list from the eRic database. The survival-related eRNAs and their target genes list are shown in Supplementary Table S2.

We conducted univariate Cox regression analysis on survival-related eRNAs target genes in training set, and 22 eRGs were significantly associated with overall survival (p < 0.05). After that, Lasso-Cox regression analysis was carried out to discover eRGs related to survival and construct a prognostic signature, which composed of 10 genes, including SSRP1, SSB, IGFBP4, SUOX, RDH16, G6PC, AKR1C3, NUP205, ADAMTS5, and RRAGD, and the correlation between these genes and their corresponding eRNAs are all significant (correlation coefficients larger than 0.3, p < 0.0001). Part of the results are shown in Figure 2, and the horizontal coordinate indicates the expression of genes in the signature, the vertical coordinate indicates the expression of their corresponding eRNAs. Other results have been shown in Supplementary Figure S1.

FIGURE 2

Figure 2 The correlation between genes in the signature and their corresponding eRNA (part results).

The risk score of each patient can be calculated by using the following formula based on the corresponding coefficients of genes: riskScore  =  0.2209  ×  Exp_SSRP1 + 0.1789  ×  Exp_SSB  +  (−0.0131)  ×  Exp_IGFBP4  +  (−0.0522)  ×  Exp_SUOX  +  (−0.0053)  ×  Exp_RDH16  +  (−0.0356)  ×  Exp_G6PC  +  0.0563  ×  Exp_AKR1C3  +  0.0106  ×  Exp_NUP205  +  0.4235  ×  Exp_ADAMTS5  +  0.0469  ×  Exp_RRAGD . Subsequently, we obtained all the risk scores of patients in the TCGA dataset and considered the median as the cutoff point to divide the patients into high- and low-risk groups. Figure 3A presents the expression profiles of 10 survival-related eRGs in the signature, and the expressions of 10 genes were all significantly different between the two groups (Supplementary Figure S2). Among them, the expressions of IGFBP4, SUOX, RDH16, and G6PC were lower in the high-risk group than those in the low-risk group, while the expressions of SSRP1, SSB, NUP205, AKR1C3, ADAMTS5, and RRAGD were higher in high-risk group. Figure 3B shows risk scores and survival status in the two groups. The Kaplan-Meier survival curves in the training set are shown in Figure 3C, which revealed that patients in the high-risk group had a significantly worse prognosis than those in the low-risk group (p < 0.0001). The time-dependent ROCs for 1-, 3-, and 5-year OS are exhibited in Figure 3D, and their AUCs (area under the ROC curve) were 0.79, 0.73, and 0.68, respectively. The C-index (concordance index) of the signature was 0.70. All the above results demonstrated that the signature composed of 10 eRGs poses a good prognostic performance.

FIGURE 3

Figure 3 (A) Heatmap of the signature genes expression profiles in training set, and the colors represent centered and scaled log2FPKM value in the row direction (FPKM is defined as fragments per kilobase of transcript per million mapped reads); (B) Survival status distribution of patients in high risk and low risk group; (C) Kaplan-Meier curves of the signature for high risk and low risk group in training set; (D) Time-dependent ROC curves of the signature for 1-, 3-, 5-year overall survival in training set.

Furthermore, the GSE14520 dataset was considered as the external validation data to confirm the performance of the eRGs signature. The expression of genes, risk scores, and survival status are shown in Figures 4A, B. The Kaplan-Meier survival curves were also significantly different (p < 0.0001) between the two different risk groups (Figure 4C). The AUCs of time-dependent ROCs for 1-, 3-, and 5-year OS were 0.71, 0.73, 0.67, respectively (Figure 4D), and the C-index was 0.68. Taken together, these results suggested that the signature had a good capability of predicting survival.

FIGURE 4

Figure 4 (A) Heatmap of the signature genes expression profiles in validation set, and the colors represent centered and scaled log2FPKM value in the row direction (FPKM is defined as fragments per kilobase of transcript per million mapped reads); (B) Survival status distribution of patients in high risk and low risk group; (C) Kaplan-Meier curves of the signature for high risk and low risk group in validation set; (D) Time-dependent ROC curves of the signature for 1-, 3-, 5-year overall survival in validation set.

3.2 Prognostic Nomogram

Univariate Cox regression analysis and multivariate Cox regression analysis were used to identify the prognostic factors in TCGA training set, and the risk score of the eRGs signature was a crucial independent prognostic predictor (Figure 5A). The results of Schoenfeld residuals test can be found in Figure S3 (the global Schoenfeld test p = 0.29), which indicated that the Cox model satisfied the proportional hazard assumption. Then a nomogram containing age, TNM stage, and risk score was established to predict 1-, 3-, and 5-year survival based on the training set (Figure 5B). We conducted a series of internal validation on the performance of the nomogram. Calibration curves presented good concordance between the nomogram-predicted survival and the actual survival of 1-, 3-, and 5-year (Figure 5C). The C-index of the nomogram was 0.73. Time-dependent ROC curves at 1-, 3-, and 5-year were exhibited in Figure 5D, and AUCs of the nomogram at 1-, 3-, and 5-year were 0.82, 0.77, 0.74, respectively. Furthermore, the model with signature had significantly higher AUC values than the model without signature, which suggested that our eRGs signature and nomogram possessed an excellent prognostic performance. DCA curves of the nomogram for 1-, 3-, 5-year survival in training set are shown in Figure 5E, which demonstrates that the nomogram had a high net benefit.

FIGURE 5

Figure 5 (A) Forest plot for univariate Cox regression analysis and multivariate Cox regression analysis; (B) Nomogram predicting 1-, 3-, 5-year survival probability; (C) Calibration curves of the nomogram for 1-, 3-, 5-year overall survival in training set; (D) Time-dependent ROC curves of the nomogram for 1-, 3-, 5-year overall survival in training set; (E) DCA curves of the nomogram for 1-, 3-, 5-year overall survival in training set.

The prognostic performance of nomogram was also validated in GEO validation set. Kaplan-Meier survival curves of the nomogram in TCGA training set and GEO validation set were presented in Supplementary Figures S4A, B, and the curves also significantly different (p < 0.0001) between the high-risk and low-risk group. The C-index of the nomogram in validation set was 0.71. AUCs of the nomogram at 1-, 3-, and 5-year were 0.74, 0.77, 0.74, respectively, and significantly larger than that model without signature (Supplementary Figure S4D). In addition, calibration curves, ROC curves, and DCA curves (Supplementary Figures S4C, E) in GEO validation set were all further confirmation that the nomogram had a good predictive value and clinical application value.

3.3 Pathways and Mechanism Analyses

We performed GO and KEGG enrichment analysis on the survival-related eRGs and identified 28 GO terms and two KEGG pathways (Supplementary Table S3 and Figure S5). These survival-related eRGs were enriched in steroid metabolic process, lipid transport, carbohydrate catabolic process, coenzyme metabolic process, oxidoreductase activity, and other GO terms. The enriched KEGG pathways were ABC transporters and steroid hormone biosynthesis. The single gene GSEA results can be seen in Figures 6A, B, Supplementary Figure S6 and Supplementary Table S4 exhibit the different pathways between the high-risk and low-risk group. The up-regulated pathways in high-risk patients were pathogenic Escherichia coli infection, cell cycle, DNA replication, mismatch repair, spliceosome, and ribosome. In comparison, the down-regulated pathways in high-risk patients were fatty acid metabolism, drug metabolism cytochrome P450, steroid hormone biosynthesis, primary bile acid biosynthesis, PPAR signaling pathway, complement and coagulation cascades, amino acid metabolism, linoleic acid metabolism, etc.

FIGURE 6

Figure 6 Enrichment analyses. (A) Single gene GSEA results; (B) GSVA results: different pathways between high risk patients and low risk patients.

3.4 Immune Landscape

eRGs we identified were associated with immune function, and G6PC, SSRP1, NUP205, ADAMTS5, and RRAGD were all related to immune response or process, so we further explored the relationship between risk scores and immune landscape.

Macrophages and T cells had a large proportion in the immune infiltration of TCGA HCC patients (Figure 7A). The compositions of immune infiltration were significantly different between the two risk groups. The proportions of B cells naïve, Macrophages M2, Monocytes, NK cells resting, T cells gamma delta were lower in the high-risk group than those in the low-risk group, while dendritic cells resting, Macrophages M0, T cells follicular helper had higher proportions in the high-risk group (Figure 7B). The abundances of B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells were estimated based on Timer algorithm, and they were all significantly correlated with risk score (Spearman correlation test, p < 0.0001, Figure 7D). In addition, it can be seen from Figure 7C that, except for HLA-B, HLA-C, HLA-E, HLA-F, and HLA-G, the expression levels of HLA-related genes were significantly different between the two groups with a higher expression in the high-risk group.

FIGURE 7

Figure 7 (A) The composition of immune infiltrating cells in TCGA-LIHC patients; (B) Differences of immune infiltration cells between high and low risk group; (C) Differences of HLA-related genes between high and low risk group; (D) The correlation between risk score and the abundances of B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells. (*p < 0.05; **p < 0.01; ***p < 0.001; “ns” means “no significant difference”).

3.5 Drug Sensitivity, Immunotherapy Effect and Transcription Factors

The expressions of four common immune checkpoints (CTLA4, PD-L1, TIGIT, and HAVCR2) were all significantly lower in the low-risk group than high-risk group (p < 0.05), which indicated that low-risk patients may have better outcomes when treated with immune checkpoint inhibitors (Figure 8A). There was a significant difference in IC50 of Sorafenib between the two different risk groups, and low-risk patients have a lower IC50, which suggested that patients in the low-risk group may be more sensitive to Sorafenib (Figure 8B).

FIGURE 8

Figure 8 (A) Differences of immune checkpoints between high risk group and low risk group; (B) Difference of Sorafenib sensitivity (IC50) between high risk and low risk patients; (C) The network diagram of the correlation between TFs and genes in the prognostic signature. Green circles represent genes in the signature and purple polygons represent TFs. The red lines indicate positive correlation and the green lines indicate negative correlation; (D) Differences of DNA damage repair related genes between high and low risk group. (**p < 0.01; ****p < 0.0001).

We downloaded a list of 318 TFs from Cistrome and identified 14 TFs related to the genes in prognostic signature (Supplementary Table S5). The network diagram of the correlation between TFs and genes is shown in Figure 8C, and SSRP1, SSB, NUP205, and RDH16 have more co-expressed TFs with positive correlations. Furthermore, we also compared the differences of DNA damage repair related genes (PMS2, EPVAM, MLH1, and MSH2) between the high-risk group and low-risk group, and the expressions of these genes were all significantly higher in the high-risk group than low-risk group (Figure 8D).

3.6 Experimental Confirmation of AKR1C3 In Vitro on Two HCC Cell Lines

We utilized the UCSC Genome Browser (http://genome.ucsc.edu/index.html) visualized the location of AKR1C3 and relative enhancers (Figure 9A). The AKR1C3 located on chromosome 10:5094414-5107686, and the enhancers associated with AKR1C3 located on chromosome 10:5060919-5067935. After transfection with AKR1C3-specific siRNA in Huh7 and MHCC-LM3 cell lines, both of the two selected siRNAs could significantly decrease AKR1C3 expression compared with control cells according to the western blot analysis (Figure 9B). EdU staining and colony formation assays were applied to assess the effect of si-AKR1C3-transfection on proliferation. The results indicated that compared with the control (si-NC), the si-AKR1C3 significantly reduced cell viability (Figure 9C) and the number of colony formations (Figure 9D). We further explored the effects of AKR1C3 on migration and invasion capacity of Huh7 and MHCC-LM3 cells via Transwell chamber assays. The migration and invasion abilities were significantly inhibited in si-AKR1C3-1 and si-AKR1C3-2 groups compared to the si-NC group (Figures 9E, F). Consistently, the wound healing assay revealed that the si-NC group had a higher efficiency at closing the wound width than those in the AKR1C3 silencing group (Figure 9G). AKR1C3 overexpression promoted and silencing inhibited protein expression (Figure 10A) and cell proliferation (Figure 10B). In comparison with the control group and sorafenib groups, the si-AKR1C3-1 plus sorafenib groups showed significantly decreased cell proliferation rate, then it is of note that the AKR1C3 plus sorafenib groups further enhanced cell proliferation rate compared with the sorafenib groups (Figure 10C). Altogether, these experimental results suggested that the knockdown of AKR1C3 inhibited the cell proliferation, migration, and invasion in HCC cell lines, which indicated that AKR1C3 plays a key role in HCC cell proliferation and aggressiveness.

FIGURE 9

Figure 9 (A) The locations of AKR1C3 and enhancers associated with AKR1C3 on chromosome. (B) Western blot analysis to examine the efficiency of the AKR1C3 knockdown. (C) Proliferation ability in AKR1C3 knockdown Huh7 and MHCC-LM3 cells by EdU staining. (D) Colony-forming abilities in AKR1C3 knockdown Huh7 and MHCC-LM3 cells by clonogenic assays. (E, F) Transwell assays to detect the migration and invasive capacities in AKR1C3 knockdown Huh7 and MHCC-LM3 cells. (G) Wound-healing assay was performed to measure the migration ability of various cells as indicated. Magnification, × 200 (C, E, F), × 40 (G). Scale bar, 100 mm (C, E, F), 500 mm (G). Data were shown as mean ± SD of at least three independent experiments. (**p < 0.01; ***p < 0.001; ****p < 0.0001).

FIGURE 10

Figure 10 (A) Western blot analysis to examine the efficiency of the AKR1C3 knockdown and overexpression. (B, C) Proliferation curves were determined in AKR1C3 knockdown and overexpression Huh7 and MHCC-LM3 cells by cell counting kit-8 (CCK-8) assays in normal culture conditions and sorafenib sensitivity experiment. Data were shown as mean ± SD of at least three independent experiments. (***p < 0.001).

4 Discussion

Enhancer RNAs (eRNAs) regulate the expression of oncogenes or tumor suppressors and play a prominent role in the tumorigenesis, progression, and proliferation of cancers, which could be used as potential and valuable diagnostic, prognostic, and therapeutic markers for cancers. Considering this, we attempted to establish a novel and robust signature based on eRGs to provide a new perspective for prognostic prediction and optimization of personalized treatment in HCC patients for the first time.

In this study, we identified a signature including 10 survival-related eRGs (SSRP1, SSB, IGFBP4, SUOX, RDH16, G6PC, AKR1C3, NUP205, ADAMTS5, and RRAGD). The signature stratified HCC patients into high- and low-risk groups, and patients in the high-risk group significantly tended to have a poorer prognosis. Previous studies have already constructed different signatures for predicting the overall survival in HCC patients (41–43). The comparison of the AUCs for 1-, 3-, 5-year overall survival of our eRGs signature and those of previous signatures were presented in Table 1. Obviously, our signature has a better prognostic value than others, especially in external validation dataset, which indicates that our signature is more robust. Furthermore, we found that the signature was an independent prognostic factor, which can combine age and TNM stage to establish a nomogram. The calibration curves of the nomogram indicated that the predicted outcomes were in good agreement with the actual outcomes in training set and validation set. Furthermore, time-dependent ROC curves exhibited that the sensitivity and specificity of model with signature were significantly improve compared with those of the model with signature. According to the nomogram, clinicians can predict the prognosis of HCC patients and provide appropriate individualized treatment to improve their quality of life.

TABLE 1

Table 1 Comparison of AUCs for prognostic signatures in different studies.

Genes in the signature and their corresponding eRNAs were all significantly correlated. Among them, ENSR00000052553 is the cancer-type-specific eRNA of HCC, which regulates the expression of SUOX and RDH16, and can be considered as a potential target eRNA for further treatment studies of HCC. In addition, G6PC, AKR1C3, and ADAMTS5 were all regulated by two neighbor eRNAs, which indicated that the enhancers at those locations may cluster into super-enhancers. Studies have already demonstrated that super-enhancers play a significant role in oncogene activation, process of tumorigenesis, and tumor cell proliferation (23–25), so our findings may lay the groundwork for investigating the mechanism of super enhancers in HCC.

Compared with the low-risk group, four genes (IGFBP4, SUOX, RDH16, and G6PC) of the prognostic signature downregulated in the high-risk group, while the other six genes (SSRP1, SSB, NUP205, AKR1C3, ADAMTS5, and RRAGD) upregulated conversely. IGFBP4 is the smallest member of human insulin-like growth factors binding proteins (IGFBPs) (44), which is involved in the inhibition of oncogenic pathways and exerts a powerful tumor suppressor function in HCC cells. It also had been demonstrated that low expression of IGFBP4 is associated with poor prognosis in HCC patients (45). SUOX was found to decrease with the progression of HCC and considered as an independent prognostic factor for overall survival and time to recurrence in HCC patients (46). RDH16 is a tumor-suppressing gene, and it had been reported that downregulation of RDH16 occurs in approximately 90% of primary HCC patients with poor prognosis (47). In addition, RDH16 was also contained in a robust twelve-gene signature for predicting survival of HCC patients (48). The deficiency of G6PC can cause glycogen storage disease type Ia (GSD-Ia), which may lead to HCC (49). A previous study has identified G6PC as a potential prognostic target in clear cell renal cell carcinoma, and its low expression associated with poor survival and aggressive progression (50). SSRP1 can regulate the proliferation and metastasis of HCC, its aberrant overexpression is related to higher serum AFP level, larger tumor size, and higher T stage of HCC patients. It has been considered as a prognostic biomarker associated with CD8+ T cell infiltration in HCC, and patients with higher expression of SSRP1 have shorter overall survival and faster recurrence (51, 52). SSB plays a significant role in DNA replication (53), its overexpression may promote the proliferation of HCC cells. Xiong et al. suggested that the upregulation of NUP205 correlated with severe TNM stage and poor survival, which demonstrated that it can be seen as a biomarker for prognostic prediction in HCC patients (54). It was found that higher ADAMTS5 expression had a significant association with development and poorer survival of HCC, and its impact on prognosis was specific for HCC among other cancer types from TCGA project (55). Furthermore, ADAMTS5 was also identified as a prominent gene in a hypoxia-related and immune-associated prognosis signature for HCC (56). RRAGD can promote cell proliferation, invasion, migration, aerobic glycolysis, and Warburg effect (an important characteristic of cancer cell metabolism) of HCC. Upregulation of RRAGD is associated with poor prognosis (57).

AKR1C3 consisted of 323 amino acids with a predicted molecular weight of 36,853 Da Like other AKR enzymes it is a soluble monomeric NAD(P)(H) dependent oxidoreductase, the enzyme that converts carbonyl groups into secondary alcohols (58). Overexpression of AKR1C3 is usually associated with prostate cancer progression, aggressiveness, and resistance to AR-targeted therapies (59). According to a previous clinical study, upregulation of AKR1C3 is an indicator of poor prognosis in HPV16-associated and HPV-negative oropharyngeal squamous cell carcinoma (OPSCC) (60). The results from our bioinformatics analysis were confirmed by a series of experiments, and we found that the silence of AKR1C3 in Huh7 and MHCC-LM3 cells can significantly inhibit HCC cell viability, clone formation, migration, invasion ability, and wound closure potential. In addition, AKR1C3 may be regulated by super-enhancers based on our study. All the above indicated that AKR1C3 is an important eRG related to the progression and prognosis of HCC, which may be a potential biomarker for HCC treatment and intervention.

The enrichment results presented that survival-related eRGs were enriched in ATPase activity, lipid transporter activity, steroid hormone biosynthesis, and ABC transporters, etc. GSVA identified a number of significant pathways that may affect prognostic outcomes of HCC patients. The downregulation of fatty acid metabolism, drug metabolism cytochrome P450, steroid hormone biosynthesis, primary bile acid biosynthesis, complement and coagulation cascades, PPAR signaling pathway, and the upregulation of cell cycle, DNA replication may be the mechanism for poorer prognosis of HCC. Therefore, our enrichment results provide new insights into the deeper investigation of molecular mechanisms and the development of targeted drugs for HCC patients.

Immune microenvironment plays an important role in the occurrence and development of tumors, which has attracted attentions of researchers. According to our results, immune infiltration and HLA may be considered important factors in exploring the specific mechanism and improving the outcomes of HCC patients. Immune checkpoints have received a great deal of attention in cancer treatment in recent years. We checked the differences of expression of immune checkpoints (CTLA4, PD-L1, TIGIT, and HAVCR2) between the high-risk group and low-risk group, and found that they all had higher expression in the high-risk group, which was in line with previous study (61).

Although Sorafenib is the standard systemic therapeutic agent available in HCC patients, the mechanism of drug resistance and the existence of heterogeneity make patients have different drug treatment effects. So, we estimated IC50 of Sorafenib in the two groups and found that the low-risk patients performed better to Sorafenib, which may contribute to the efficient and rational medication of HCC patients. Co-expression analysis of genes and TFs showed that SSRP1 and SSB, NUP205 had positive correlations with their co-expression TFs, while G6PC and RDH16 were negatively correlated with their co-expression TFs. Meanwhile, expression of SSRP1, SSB, and NUP205 were all positively correlated with risk score, and expression of G6PC and RDH16 were negatively correlated with risk score. These findings suggested that eRNAs may bind to these TFs to regulate the expression of genes in the prognostic signature, which seemed to shed light on the investigation of HCC pathological mechanism and therapeutic targets.

Taken together, some advantages of our study deserve to be underscored. First, we constructed prognostic signature and nomogram for HCC based on eRGs for the first time, and the model performed a better predictive ability and provided a novel direction for the underlying pathological mechanism of HCC. Second, we identified survival-related eRNAs directly in HCC patients rather than select them from differential expressed eRNAs between patients and normal controls, so that more comprehensive information can be considered. Third, signature composed of specific genes is more economical and practical than whole-genome sequencing, and easy to routinely test. Fourth, the visualization of the nomogram is more convenient to assist clinicians in predicting patients’ prognosis and customizing individualized treatment plans. Last but not least, we confirmed the role of AKR1C3 in the progression and invasion of HCC through a series of in vitro experiments. Nonetheless, some limitations still exist in this study. Firstly, the external validation set was also from a public database, and a multicenter cohort study is needed to further prove the generalizability of the models. Secondly, only one eRG in the signature was experimentally verified in the present study, and other eRGs should be validated in the future studies. Thirdly, the regulatory relationships between eRNAs and their target genes from the eRic database used in our study were generated based on data-driven correlations, and their biological regulatory relationships need to be confirmed by rigorous experiments. Finally, multi-omics data such as methylation, long non-coding RNA, and proteomics should be integrated and analyzed to comprehensively elucidate biological regulatory networks.

5 Conclusion

Our study has been the first attempt to identify a prognostic signature composed of eRGs for HCC. We also constructed a nomogram incorporating the signature and clinical characteristics to predict overall survival of HCC patients accurately and robustly, which has been validated in external dataset. The signature and nomogram both performed good prognostic ability. AKR1C3 may be a potential biomarker for HCC treatment and intervention through a series of in vitro experiments. The findings of this study have significant practical implications in terms of providing a deeper insight into the investigation of pathogenesis of HCC, optimizing individualized treatment, and improving the prognosis of HCC patients.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/); enhancer RNA in cancers (eRic) (https://hanlab.uth.edu/eRic/); Gene Expression Omnibus, GSE14520 -GPL3921.

Author Contributions

Conceptualization, WZ and WT; methodology, WZ and KC; software, WZ; validation, QiZ; formal analysis, WZ and WT; investigation, LS and YW; resources, KC; data curation, WZ and KC; writing—original draft preparation, WZ and KC; writing—review and editing, QiuZ; visualization, WZ; supervision, QiuZ; project administration, LS and YW; funding acquisition, ML. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Science and Technology Major Project (2016ZX08011005), the National Natural Science Foundation of China (82073666), and the National Natural Science Foundation of China (82003556).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are very grateful to the staff who provided TCGA, GEO database and other online resources for providing us with a research foundation.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.849242/full#supplementary-material

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Singal AG, Elserag HB. Hepatocellular Carcinoma From Epidemiology to Prevention: Translating Knowledge Into Practice. Clin Gastroenterol Hepatol (2015) 13:2140–51. doi: 10.1016/j.cgh.2015.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Rastogi A. Changing Role of Histopathology in the Diagnosis and Management of Hepatocellular Carcinoma. World J Gastroenterol (2018) 24(35):4000–13. doi: 10.3748/wjg.v24.i35.4000

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Llovet J, Brú C, Bruix J. Prognosis of Hepatocellular Carcinoma: The BCLC Staging Classification. Semin Liv Dis (1999) 19(3):329-38. doi: 10.1055/s-2007-1007122

CrossRef Full Text | Google Scholar

5. Yau T, Tang VY, Yao TJ, Fan ST, Lo CM, Poon RT. Development of Hong Kong Liver Cancer Staging System With Treatment Stratification for Patients With Hepatocellular Carcinoma. Gastroenterology (2014) 146(7):16917–700.e3. doi: 10.1053/j.gastro.2014.02.032

CrossRef Full Text | Google Scholar

6. Li L, Wang H. Heterogeneity of Liver Cancer and Personalized Therapy. Cancer Lett (2016) 379(2):191–7. doi: 10.1016/j.canlet.2015.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Yang JD, Roberts LR. Epidemiology and Management of Hepatocellular Carcinoma. Infect Dis Clin North Am (2010) 24:899–919. doi: 10.1016/j.idc.2010.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Pan HW, Ou YH, Peng SY, Liu SH, Lai PL, Lee PH, et al. Overexpression of Osteopontin Is Associated With Intrahepatic Metastasis, Early Recurrence, and Poorer Prognosis of Surgically Resected Hepatocellular Carcinoma. Cancer (2003) 98(1):119–27. doi: 10.1002/cncr.11487

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Giordano S, Columbano A. Met as a Therapeutic Target in HCC: Facts and Hopes. J Hepatol (2014) 60:442–52. doi: 10.1016/j.jhep.2013.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Toyoda H, Kumada T, Kiriyama S, Sone Y, Tanikawa M, Hisanaga Y, et al. Prognostic Significance of Simultaneous Measurement of Three Tumor Markers in Patients With Hepatocellular Carcinoma. Clin Gastroenterol Hepatol (2006) 4(1):111–7. doi: 10.1016/S1542-3565(05)00855-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. He Q, Lin Z, Wang Z, Huang W, Tian D, Liu M, et al. SIX4 Promotes Hepatocellular Carcinoma Metastasis Through Upregulating YAP1 and C-MET. Oncogene (2020) 39(50):7279–95. doi: 10.1038/s41388-020-01500-y

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Huang W, Yu D, Wang M, Han Y, Lin J, Wei D, et al. ITGBL1 Promotes Cell Migration and Invasion Through Stimulating the TGF-β Signalling Pathway in Hepatocellular Carcinoma. Cell Prolif (2020) 53(7):e12836. doi: 10.1111/cpr.12836

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hong X, Luo H, Zhu G, Guan X, Jia Y, Yu H, et al. SSR2 Overexpression Associates With Tumorigenesis and Metastasis of Hepatocellular Carcinoma Through Modulating EMT. J Cancer (2020) 11(19):5578–87. doi: 10.7150/jca.44788

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Wang TH, Yu CC, Lin YS, Chen TC, Yeh CT, Liang KH, et al. Long Noncoding RNA CPS1-IT1 Suppresses the Metastasis of Hepatocellular Carcinoma by Regulating HIF-1α Activity and Inhibiting Epithelial-Mesenchymal Transition. Oncotarget (2016) 7(28):43588–603. doi: 10.18632/oncotarget.9635

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang Z, Weng W, Huang W, Wu B, Zhou Y, Zhang J, et al. A Novel Molecular-Clinicopathologic Nomogram to Improve Prognosis Prediction of Hepatocellular Carcinoma. Aging (2020) 12(13):12896–920. doi: 10.18632/aging.103350

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Liu GM, Xie WX, Zhang CY, Xu JW. Identification of a Four-Gene Metabolic Signature Predicting Overall Survival for Hepatocellular Carcinoma. J Cell Physiol (2020) 235:1624–36. doi: 10.1002/jcp.29081

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhao QJ, Zhang J, Xu L, Liu FF. Identification of a Five-Long Non-Coding RNA Signature to Improve the Prognosis Prediction for Patients With Hepatocellular Carcinoma. World J Gastroenterol (2018) 24:3426–39. doi: 10.3748/wjg.v24.i30.3426

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, Wu J, et al. Widespread Transcription at Neuronal Activity-Regulated Enhancers. Nature (2010) 465(7295):182–7. doi: 10.1038/nature09033

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Li W, Notani D, Rosenfeld MG. Enhancers as Non-Coding RNA Transcription Units: Recent Insights and Future Perspectives. Nat Rev Genet (2016) 17:207–7. doi: 10.1038/nrg.2016.4

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li W, Hu Y, Oh S, Ma Q, Merkurjev D, Song X, et al. Condensin I and II Complexes License Full Estrogen Receptor α-Dependent Enhancer Activation. Mol Cell (2015) 59(2):188–202. doi: 10.1016/j.molcel.2015.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hsieh CL, Fei T, Chen Y, Li T, Gao Y, Wang X, et al. Enhancer RNAs Participate in Androgen Receptor-Driven Looping That Selectively Enhances Gene Activation. Proc Natl Acad Sci USA (2014) 111(20):7319–24. doi: 10.1073/pnas.1324151111

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Melo CA, Drost J, Wijchers PJ, van de Werken H, de Wit E, Oude Vrielink JA, et al. eRNAs Are Required for P53-Dependent Enhancer Activity and Gene Transcription. Mol Cell (2013) 49(3):524–35. doi: 10.1016/j.molcel.2012.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lee JH, Xiong F, Li W. Enhancer RNAs in Cancer: Regulation, Mechanisms and Therapeutic Potential. RNA Biol (2020) 17:1550–9. doi: 10.1080/15476286.2020.1712895

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hnisz D, Schuijers J, Lin CY, Weintraub AS, Abraham BJ, Lee TI, et al. Convergence of Developmental and Oncogenic Signaling Pathways at Transcriptional Super-Enhancers. Mol Cell (2015) 58(2):362–70. doi: 10.1016/j.molcel.2015.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Lovén J, Hoke HA, Lin CY, Lau A, Orlando DA, Vakoc CR, et al. Selective Inhibition of Tumor Oncogenes by Disruption of Super-Enhancers. Cell (2013) 153(2):320–34. doi: 10.1016/j.cell.2013.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M, et al. An Atlas of Active Enhancers Across Human Cell Types and Tissues. Nature (2014) 507(7493):455–61. doi: 10.1038/nature12787

PubMed Abstract | CrossRef Full Text | Google Scholar

27. ENCODE Project Consortium. An Integrated Encyclopedia of DNA Elements in the Human Genome. Nature (2012) 489(7414):57–74. doi: 10.1038/nature11247

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic A, Meissner A, et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat Biotechnol (2010) 28(10):1045–8. doi: 10.1038/nbt1010-1045

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang Z, Lee JH, Ruan H, Ye Y, Krakowiak J, Hu Q, et al. Transcriptional Landscape and Clinical Utility of Enhancer RNAs for eRNA-Targeted Therapy in Cancer. Nat Commun (2019) 10(1):4562. doi: 10.1038/s41467-019-12543-5

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Grimes DA. The Nomogram Epidemic: Resurgence of a Medical Relic. Ann Internal Med (2008) 149:273–5. doi: 10.7326/0003-4819-149-4-200808190-00010

CrossRef Full Text | Google Scholar

31. Balachandran VP, Gonen M, Smith JJ, DeMatteo RP. Nomograms in Oncology: More Than Meets the Eye. Lancet Oncol (2015) 16:e173–80. doi: 10.1016/s1470-2045(14)71116-7

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hnzelmann S, Castelo R, Guinney J. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7–7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

34. Liu X, Qin S. Immune Checkpoint Inhibitors in Hepatocellular Carcinoma: Opportunities and Challenges. Oncology (2019) 24:S3. doi: 10.1634/theoncologist.2019-IO-S1-s01

CrossRef Full Text | Google Scholar

35. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–e110. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Aykul S, Martinez-Hackert E. Determination of Half-Maximal Inhibitory Concentration Using Biosensor-Based Protein Interaction Analysis. Analytic Biochem (2016) 508:97–103. doi: 10.1016/j.ab.2016.06.025

CrossRef Full Text | Google Scholar

38. Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc JF, et al. Sorafenib in Advanced Hepatocellular Carcinoma. N Engl J Med (2008) 359(4):378–90. doi: 10.1056/NEJMoa0708857

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): A Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res (2013) 41(Database issue):D955–61. doi: 10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zhang BH, Yang J, Jiang L, Lyu T, Kong LX, Tan YF, et al. Development and Validation of a 14-Gene Signature for Prognosis Prediction in Hepatocellular Carcinoma. Genomics (2020) 112(4):2763–71. doi: 10.1016/j.ygeno.2020.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Li B, Feng W, Luo O, Xu T, Cao Y, Wu H, et al. Development and Validation of a Three-Gene Prognostic Signature for Patients With Hepatocellular Carcinoma. Sci Rep (2017) 7(1):5517. doi: 10.1038/s41598-017-04811-5

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhang FP, Huang YP, Luo WX, Deng WY, Liu CQ, Xu LB, et al. Construction of a Risk Score Prognosis Model Based on Hepatocellular Carcinoma Microenvironment. World J Gastroenterol (2020) 26(2):134–53. doi: 10.3748/wjg.v26.i2.134

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Chelius D, Baldwin MA, Lu X, Spencer EM. Expression, Purification and Characterization of the Structure and Disulfide Linkages of Insulin-Like Growth Factor Binding Protein-4. J Endocrinol (2001) 168:283–96. doi: 10.1677/joe.0.1680283

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Lee YY, Mok MT, Kang W, Yang W, Tang W, Wu F, et al. Loss of Tumor Suppressor IGFBP4 Drives Epigenetic Reprogramming in Hepatic Carcinogenesis. Nucleic Acids Res (2018) 46(17):8832–47. doi: 10.1093/nar/gky589

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Jin GZ, Yu WL, Dong H, Zhou WP, Gu YJ, Yu H, et al. SUOX Is a Promising Diagnostic and Prognostic Biomarker for Hepatocellular Carcinoma. J Hepatol (2013) 59(3):510–7. doi: 10.1016/j.jhep.2013.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhu YH, Li JB, Wu RY, Yu Y, Li X, Li ZL, et al. Clinical Significance and Function of RDH16 as a Tumor-Suppressing Gene in Hepatocellular Carcinoma. Hepatol Res Off J Jap Soc Hepatol (2020) 50(1):110–20. doi: 10.1111/hepr.13432

CrossRef Full Text | Google Scholar

48. Ouyang G, Yi B, Pan G, Chen X. A Robust Twelve-Gene Signature for Prognosis Prediction of Hepatocellular Carcinoma. Cancer Cell Int (2020) 20:207. doi: 10.1186/s12935-020-01294-9

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Cho JH, Lee YM, Bae SH, Chou JY. Activation of Tumor-Promoting Pathways Implicated in Hepatocellular Adenoma/Carcinoma, A Long-Term Complication of Glycogen Storage Disease Type Ia. Biochem Biophys Res Commun (2020) 522:1–7. doi: 10.1016/j.bbrc.2019.11.061

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Xu WH, Xu Y, Tian X, Anwaier A, Liu WR, Wang J, et al. Large-Scale Transcriptome Profiles Reveal Robust 20-Signatures Metabolic Prediction Models and Novel Role of G6PC in Clear Cell Renal Cell Carcinoma. J Cell Mol Med (2020) 24(16):9012–27. doi: 10.1111/jcmm.15536

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Luo G, Xu J, Xia Z, Liu S, Liu H, He K, et al. SSRP1 Is a Prognostic Biomarker Correlated With CD8(+) T Cell Infiltration in Hepatocellular Carcinoma (HCC). BioMed Res Int (2021) 2021:9409836. doi: 10.1155/2021/9409836

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Ding Q, He K, Luo T, Deng Y, Wang H, Liu H, et al. SSRP1 Contributes to the Malignancy of Hepatocellular Carcinoma and Is Negatively Regulated by miR-497. Mol Ther J Am Soc Gene Ther (2016) 24(5):903–14. doi: 10.1038/mt.2016.9

CrossRef Full Text | Google Scholar

53. Pandita RK, Chow TT, Udayakumar D, Bain AL, Cubeddu L, Hunt CR, et al. Single-Strand DNA-Binding Protein SSB1 Facilitates TERT Recruitment to Telomeres and Maintains Telomere G-Overhangs. Cancer Res (2015) 75(5):858–69. doi: 10.1158/0008-5472.Can-14-2289

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Xiong DD, Feng ZB, Lai ZF, Qin Y, Liu LM, Fu HX, et al. High Throughput circRNA Sequencing Analysis Reveals Novel Insights Into the Mechanism of Nitidine Chloride Against Hepatocellular Carcinoma. Cell Death Dis (2019) 10(9):658. doi: 10.1038/s41419-019-1890-9

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhu Z, Xu J, Wu X, Lin S, Li L, Ye W, et al. In Silico Identification of Contradictory Role of ADAMTS5 in Hepatocellular Carcinoma. Technol Cancer Res Treat (2021) 20:1533033820986826. doi: 10.1177/1533033820986826

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Hu B, Yang XB, Sang XT. Development and Verification of the Hypoxia-Related and Immune-Associated Prognosis Signature for Hepatocellular Carcinoma. J Hepatocell Carcinom (2020) 7:315–30. doi: 10.2147/jhc.S272109

CrossRef Full Text | Google Scholar

57. Ding L, Liang X. Ras Related GTP Binding D Promotes Aerobic Glycolysis of Hepatocellular Carcinoma. Ann Hepatol (2021) 23:100307. doi: 10.1016/j.aohep.2021.100307

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Jin Y, Penning TM. Aldo-Keto Reductases and Bioactivation/Detoxication. Annu Rev Pharmacol Toxicol (2007) 47:263–92. doi: 10.1146/annurev.pharmtox.47.120505.105337

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Liu C, Yang JC, Armstrong CM, Lou W, Liu L, Qiu X, et al. AKR1C3 Promotes AR-V7 Protein Stabilization and Confers Resistance to AR-Targeted Therapies in Advanced Prostate Cancer. Mol Cancer Ther (2019) 18(10):1875–86. doi: 10.1158/1535-7163.Mct-18-1322

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Huebbers CU, Verhees F, Poluschkin L, Olthof NC, Kolligs J, Siefer OG, et al. Upregulation of AKR1C1 and AKR1C3 Expression in OPSCC With Integrated HPV16 and HPV-Negative Tumors Is an Indicator of Poor Prognosis. Int J Cancer (2019) 144(10):2465–77. doi: 10.1002/ijc.31954

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Liu Z, Jiao D, Liu L, Zhou X, Yao Y, Li Z, et al. Development and Validation of a Robust Immune-Related Risk Signature for Hepatocellular Carcinoma. Medicine (2021) 100(10):e24683. doi: 10.1097/md.0000000000024683

PubMed Abstract | CrossRef Full Text | Google Scholar