Six hundred forty-three HCC patients from TCGA, GEO, and ICGC were enrolled. A total of 370 HCC patients with comprehensive transcriptional and clinical data from the TCGA database (https://portal.gdc.cancer.gov/) were used as a training cohort. In this study, two external cohorts of HCC were used as validation groups. In validation group 1, 243 HCC patients were selected from the ICGC database (https://dcc.icgc.org/projects/). The validation group 2 consists of 30 HCC patients from the GSE107943 database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE107943).

Reference genes were chosen from the Molecular Signatures Database (MSigDB) (http://www.gsea-msigdb.org/gsea/msigdb/search.jsp) (Table 1).

We calculated GSEA in the training group using the R package “GSVA” based on the above-mentioned gene set (Lee et al., 2008). R package for WGCNA was used to evaluate TCGA database mRNA matrix data (Langfelder and Horvath, 2008). We generated an adjacency matrix and then converted it into a topological overlap matrix (TOM) to assess the correlation strength between the nodes. By performing hierarchical clustering, we ensured that each module contained at least 50 genes. Finally, we integrated similar modules to find module genes having a good correlation with ssGSEA E2F target scores. By combining these genes with the E2F target gene set in HALLMARKS, we got a new gene set named gene set A.

Univariate Cox regression was used to identify gene set B’s prognosis-related genes from gene set A. We used R package “limma” to determine DEGs between TCGA-LIHC tumors and healthy controls with a false discovery rate [FDR] < 0.05 and a log|fold change [FC]|> 1) (Ritchie et al., 2015). E2F-related DEGs and gene set B were intersected to obtain candidate genes to construct the gene signature. To avoid overfitting, LASSO Cox regression (R package “glmnet”) was conducted to exclude collinear genes (T and ibshirani, 1997). Finally, the three best model genes were selected, and their coefficients were recorded to create an E2F target gene signature.

According to the contents of the “hist_hepato_carc_fact” record in the clinical information provided in the UCSC database, liver cancer samples were divided into HBV, HCV, ALD, and NASH. Finally, 104 samples with Hepatitis B, 56 samples with Hepatitis C, and 20 samples with Non-Alcoholic Fatty Liver Disease were retrieved. Alcoholic liver disease (ALD) related content was not obtained. Therefore, we divided the remaining 200 samples into others. Further, we compared the distribution of samples of different subgroups in the high-low-risk group, and no significant difference was observed (Table 2, p = 0.533).

We use one Class Linear Regression (OCLR) to quantify the stemness of tumor samples. Two stemness indices were constructed from stem cell transcriptome, methylation group, and multi-platform data: mRNA expression-based stemness index (mRNAsi) represents gene expression, and epigenetically regulated-mRNAsi (DNAsi) measures epigenetics characteristics of stem cells. A stemness index (si) is a measure of stemness ranging from low (0) to high (1).

By developing regression models from cell line and gene expression profiles of Genomics of Drug Sensitivity in Cancer (GDSC), the pRRophetic algorithm predicted drug maximum 50% inhibitory concentration (IC50): (www.cancerrxgene.org/).

All analyses and graphs were constructed using R v4.1.1. We calculated the gene signature’s inter-gene correlations using the Pearson correlation test. The training set (TCGA cohort) and two test sets (ICGC and GEO cohorts) were divided into high- and low-risk groups based on the median value of the E2F target-related risk score. Group differences were calculated using the t-test and Chi-square test. The survival distribution and gene expression patterns of HCC patients were visualized using scatter plots and heat maps. Kaplan-Meier survival analysis determined whether high- and low-risk groups had distinct prognoses. A time-dependent ROC curve was calculated with R package time ROC. To test the independent predictability of the gene signature, the R package “survival” was used to run multivariate Cox stepwise regression models.

We used GSEA v4.1.0 to examine whether the E2F target pathway was activated in different risk categories. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to explore the functions of the DEGs in high- and low-risk groups (Subramanian et al., 2005).

The CIBERSORT algorithm (https://cibersortx.stanford.edu/) and the immunophenoscore (IPS) were used to evaluate immune cells infiltrated in high- and low-risk groups (Newman et al., 2015; Charoentong et al., 2017). Finally, we used CellMiner to find new medications and biological targets based on gene-drug sensitivity (https://discover.nci.nih.gov/cellminer/home.do) (Reinhold et al., 2012).

Figure 1A illustrates the overall workflow. Using univariate Cox regression, we calculated 50 pathways ssGSEA scores in 370 TCGA HCC samples. The top 20 pathways with statistical significance are listed in Figure 1B, among which MYC_TARGET_V1, G2M_CHECKPOINT, and E2F_TARGETs are the first three items (p < 0.0001). Next, we compared the expression scores of these pathways in tumor tissue and paracancer tissue. E2F target showed the most significant score difference between tumor and paracancer tissue (Figure 1C). The score of G2M_CHECKPOINT in the tumor was also higher than that in para-cancer tissue, and the difference between tumor and para-cancer was the second largest (Figure 1C). The score of MYC_TARGET_V1 in tumors was significantly lower than in adjacent tissues, as shown in Figure 1C. By applying Decision Curve Analysis (DCA), we found that MYC_TARGET_V1 performed slightly better than the other two pathways, and the evaluation effects of the E2F target and G2M_CHECKPOINT were very similar, as shown in Figure 1D. Then, we performed ROC curve analysis, and the results showed that G2M_CHECKPOINT had the highest AUC value, but the AUC values of the three channels were all above 0.71 (Figure 1E). There was no significant difference in ROC diagnostic performance among the three channels. Based on the above results, we finally choose the E2F target as the main content of this study. According to the median value (0.4072497) of the E2F target pathway, we classified HCC patients as low- and high-score groups. Kaplan–Meier survival curves showed that high-score patients have a shorter overall survival time than low-score individuals (p < 0.001, Figure 1F).

WGCNA co-expression algorithm identified co-expressed coding genes and modules of 370 HCC patients from the TCGA database. As a first step, we clustered the samples using hierarchical clustering and calculated the distance between each gene using the Pearson correlation coefficient. We used WGCNA to build a weight co-expression network with a soft threshold 10 to identify co-expression modules. To ensure that the network is scale-free, β= 10 is selected (Supplementary Figure S1). The expression matrix was converted into an adjacency matrix, which was then converted into a topological matrix. A hierarchical clustering method based on average linkage was employed to cluster genes, and each gene network module had a minimum of 50 genes per mixed dynamic shear tree standard. Cluster analysis was performed on the gene modules after the eigengenes of each module were calculated. Modules with relatively close distances were grouped as new modules.

When setting minModuleSize = 50, deepSplit = 3, and height = 0.25, we got a sum of 10 modules (Figure 2A). Grey modules are unaggregatable gene sets. Moreover, we analyzed the correlation between each module and each subtype and identified 609 genes with high correlation coefficients (above 0.5) with E2F, including the red (84), pink (82), and brown modules (443). These 609 genes and MSigDB E2F target pathway genes yielded 809 E2F target pathway genes. Using univariate Cox regression analysis, 576 E2F-related prognostic genes were identified (p < 0.05). Meanwhile, edgR and Deseq2 were used to analyze the differences between liver and para-cancer samples. A total of 1,400 DEGs (1,222 upregulated and 178 downregulated in the tumor) were identified by edgR, while 1,313 DEGs (1,100 up-expressed and 213 down-expressed) were obtained by Deseq2. These differential genes and 576 E2F-related prognostic genes revealed 52 differential prognostic genes (Figure 2C). Figure 2D exhibited the univariate Cox regression results of these 52 genes. Figure 2E shows the expression heatmap of 52 genes in HCC and adjacent liver tissue. Using the R software package glmnet, we conducted a lasso cox regression analysis. In Figures 2F,G, we analyzed the confidence interval for each lambda. Figure 2H showed that when lambda = 0.07227, there are three genes remained in the model, namely, growth hormone receptor (GHR), thyroid hormone receptor interactor 13 (TRIP13), and cell division cycle associated 8 (CDCA8). The final model formula is: RiskScore = −   0.0302 * exp ⁡ ⁡ G H R + 0.0694 * exp ⁡ ⁡ T R I P 13 + 0.1701 * exp ⁡ ⁡ C D C A 8

The formula suggested that increased GHR expression is a protective factor in HCC and related to low risk, while TRIP13 and CDCA8 were associated with high RiskScore and poor outcomes. TRIP13 accelerates the mitotic process and leads to chromosome instability via playing a role in spindle assembly checkpoint and DNA repair pathways (Lu et al., 2019). High TRIP13 gene expression is found in HCC tissues, which promotes cell growth and metastasis via activating AKT/mTOR and silencing TGF-β1/smad3 pathway (Yao et al., 2018; Zhu et al., 2019). Many researchers have identified CDCA8 as a novel oncogene, which predicts a poor prognosis in HCC as well as promotes tumor proliferation via MEK/ERK, AKT/β-Catenin, and CDK1/cyclin B1 signaling (Jeon et al., 2021; Cui and Jiang, 2023). In the terms of the GHR gene, some researchers have observed its upregulation promoted HCC development (Garcia-Caballero et al., 2000; Haque et al., 2022). On the other hand, its downregulation has been noted in HCV-induced HCC and related to an unfavorable outcome (Lin et al., 2021; Abu El-Makarem et al., 2022).

We searched the clinical information of the TCGA HCC cohort (n = 380) from the UCSC database. As a result, 104 samples with Hepatitis B, 56 samples with Hepatitis C, and 20 samples with Non-Alcoholic Fatty Liver Disease were retrieved. Alcoholic liver disease (ALD) related content was not obtained. We divided the remaining 200 samples into others. Next, we analyzed the prognostic KM curves of the three genes in the analysis model in different subpopulations, and the results showed that GHR showed low expression and poor prognosis in HCV samples and total samples, which was consistent with the trend presented in our model (Supplementary Figures S2A–D). TRIP13 and CDCA8 only showed high expression and poor prognosis in total samples, which was also consistent with the trend shown in our model (Supplementary Figures S2A–D). Finally, we observed the expression of the three genes in different subgroups, and the results showed that the expression of GHR in different subgroups was significantly lower than that in adjacent tissues. TRIP13 and CDCA8 were significantly higher than paracancer tissues, as shown in Supplementary Figure S2E.

We generated the risk score for each training group sample and presented the RiskScore distribution (Figure 3A). Patients with a high RiskScore have a worse prognosis than those with a low RiskScore for HCC (Figure 3A). TRIP13 and CDCA8 were risk factors whose expression altered with increasing risk value, whereas high GHR expression was a protective factor linked with low risk (Figure 3B). Using GSEA, it was confirmed that genes associated with the E2F target pathway were significantly enriched in high-risk individuals (Figure 3C). We divided TCGA HCC samples into high- and low-risk groups based on the median risk score (0.528115). The KM curve displayed that the OS of the high-risk group was significantly lower than the low-risk group [Figure 3D, log-rank p < 0.0001, HR = 2.044 (1.435–2.912)]. Using R software package timeROC, we assessed prediction power over one, two, three, and 5 years. This risk model has a very high AUC area below the line (AUC >0.69, Figure 3E).

Two published papers have constructed E2F-related models in liver cancer, including E2F target gene characteristics of five genes constructed by Hu et al. (2022), and two gene models constructed by Wang et al. (2023), respectively. To figure out whether our new E2F target signature has an advantage, we calculated the risk score of each sample according to the formulas provided in the paper. 371 TCGA HCC samples and corresponding prognostic follow-up information provided in TCGA were used for the KM curve and AUC analysis. The results showed that the Kaplan-Meier curves of all models showed significant differences, as shown in Supplementary Figure S3A. The AUC value of the model constructed by Hu et al. (2022) was above 0.658. The AUC value of the model constructed by Wang et al. (2023) is above 0.658, as shown in Supplementary Figure S3B; Compared with the above results, we found that the AUC values of the two models are lower than ours. Also, through c-index analysis, the c-index value of our model is 0.68 (Cl 95%:0.63–0.73), and the c-index value of the model constructed by Hu et al. (2022) is 0.65 (Cl 95%: 0.63–0.71). The c-index value of the model Wang et al. (2023) constructed was 0.64 (Cl 95%:0.59–0.69). Finally, through the DCA decision curve analysis, we found that the return rate of the E2F model we built was also higher than that of the other two models, as shown in Supplementary Figure S3C.

ICGC and GEO HCC cohorts were used as external validation sets to evaluate this risk model’s predictive value. First, the risk score of each sample was calculated in two datasets, and the RiskScore distribution of the samples was plotted as shown in Figures 4A, B. Next, we separated patients in two validation cohorts into high- and low-risk groups according to the median risk score (ICGC cohort: 0.3101617, GEO cohort: 0.7349518). Kaplan-Meier plotter revealed that patients in high-risk groups have poor prognoses within 6 years in the ICGC cohort and 8 years in the GEO cohort (Figures 4C, D). The prediction power was evaluated for one, two, three, and 5 years using the R software package timeROC. Figures 4E, F showed that the AUC area below the line is very high for the risk model, with 0.76 for the ICGC cohort and 0.657 for the GEO cohort.

Pearson correlation analysis assessed whether risk groups and pathological characteristics were correlated. According to our findings, advanced T Stage, advanced stage, and high grade were significantly correlated with high-risk groups (p < 0.05, as shown in Figure 5).

We calculated the DEGs between high- and low-risk groups using the “limma” algorithm, and the threshold was set as p < 0.05. Figure 6A shows 473 upregulated and 115 downregulated genes in the high-risk group. To analyze functional enrichment in DEGs, we used R’s clusterProfiler package and parameterized the threshold to an FDR of 0.05. According to KEGG analysis, the high-risk group showed activation of the DNA replication pathway, cell cycle pathway, and p53 signaling pathway (Figures 6B, C). Based on the GO analysis, E2F targets are primarily involved in DNA replication and cell cycle checkpoints (Figure 6D). The association between risk ratings and tumour stem cell features were also evaluated. Risk score has no relationship with DNA stem index (mDNAsi) but is significantly correlated with RNA stem index (mRNAsi), as shown in Figures 6E, F.

CIBERSORT was utilized as an online tool to deconvolute expression matrices of human immune cell subtypes using linear support vector regression (Newman et al., 2015). We characterized the differences of 22 immune cell infiltration between high- and low-risk groups. In the high-risk group, there were more infiltrating B cell plasma, T cell follicular helper, T cell regulatory (Tregs), macrophage M0, neutrophils, and T cell CD4⁺ memory activated within tumor tissues than in the low-risk group (Figure 7A). In contrast, the high-risk group had fewer T cell CD4⁺ memory resting, monocyte, NK cell resting, mast cell activated, B cell naïve infiltrated compared to the low-risk group (Figure 7A). According to the CIBERSORT algorithm, the high-risk group displayed immunosuppressive infiltration and a lack of immunoactive cells. Moreover, we analyzed the expression levels of 41 immune suppressor genes from the tumor immunophenotype database (TIP, http://biocc.hrbmu.edu.cn/TIP/index.jsp). Results revealed that 31 genes out of 41 (75.6%) were significantly upregulated in the high-risk group, including programmed cell death 1 (PDCD1), hepatitis A virus cellular receptor 1/2 (HAVCR1/2), T cell immunoreceptor with Ig and ITIM domains (TIGIT), indoleamine 2,3-dioxygenase 1 (IDO1), enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), and DNA methyltransferase 1 (DNMT1) (Figure 7B).

ESTIMATE predicts tumor purity and infiltrated matrix/immune cells in tumor tissue based on ssGSEA. Estimate score (tumor purity), stromal score (matrix in tumor tissue), and immune score (infiltration of immune cells in tumor tissue) are the three main scores derived from ESTIMATE. Figure 7C shows the proportion of stromal scores is significantly lower in a high-risk group than in the low-risk group. Immune and ESTIMATE Scores had no obvious difference between the two groups (Figures 7D, E). The results suggest that the stromal component of TME is more suitable for distinguishing high-risk from low-risk patients.

The TCGA data was used by Thorsson et al. (2018) to identify six immune subtypes: C1 (Wound Healing), C2 (IFN-γ Dominant), C3 (Inflammatory), C4 Lymphocyte Depleted), C5 (Immunologically Quiet), and C6 (TGF-β Dominant). The risk score of C1, C2, and C4 types with poor prognoses was higher, and C3 types with good prognoses were primarily found among patients with low-risk scores (Figure 7F). This phenomenon indicated that the immune subtype and risk score are related.

Drug resistance often emerges during cancer treatment, leading to poor efficacy and unfavorable outcome of HCC. To test E2F risk models in chemotherapy, we predicted the IC50 of 138 medicines in high- and low-risk patients using the pRRophetic algorithm. A total of 65 drugs showed significant variations between high-risk and low-risk patients, with Tipifarnib, Camptothecin, Salubrinal, and Nilotinib being more sensitive to high-risk groups. (Figure 8A, p < 0.05). Moreover, we found that most 65 drugs showed a significant correlation between their IC50s and risk scores (Figure 8B). The above results suggest that the E2F risk model can help HCC patients choose chemotherapy medicines based on clinical drug effects.