Comprehensive computerized searches of three publicly available datasets were conducted to procure the messenger RNA (mRNA) expression profiles. The TCGA LIHC cohort (28) included 370 patients and the ICGC LIRI-JP cohort included 232 samples (29). The microarray datasets, including 225 samples of GSE14520, were downloaded from the GEO. A total of 870 CRs were retrieved from previous research (4). The “limma” R package was utilized to select the CRs related to differentially expressed genes (DEGs) between nontumor and tumor tissues in the TCGA LIHC cohort, with a P-value of<0.05 and |log2FC| of ≥0.2.

WGCNA was applied for pinpointing the HCC clinical characteristic-specific module by running the R package “WGCNA” (30, 31). The expression profiles of CR-related DEGs were utilized as an import for the WGCNA, and clinical characteristics were analyzed and defined as the sample phenotype. The power of β = 10 and scale-free R2 of 0.95 was instituted as the soft threshold parameters to ensure a signed scale-free co-expression gene network. Correlations were calculated between the module eigengenes and clinical information based on the eigengenes function. Several hub genes were considered functionally significant because they were markedly interconnected with nodes in a module. Our study selected an attractive module and identified hub genes by clinical trait significance and module connectivity.

A co-expression network based on the selected module was constructed by the exportNetworkToCytoscape function in the WGCNA R package and visualized in Cytoscape software to obtain hub nodes (32). Hub genes were calculated by applying the cytoHubba plugin based on the maximal clique centrality (MCC) algorithm (33).

A consensus clustering algorithm was applied based on the hub genes to confirm the number of clusters in the TCGA cohort and further validated in the ICGC-LIRI-JP cohort and GSE14520. This step was run and repeated 1,000 times in R using the package ConsensusClusterPlus to guarantee classification stability.

A single sample gene set enrichment analysis (ssGSEA) was used to disclose the relative amount of infiltration of 28 immune cells in the TME according to the TCGA-HCC dataset (34). The marker gene sets for TME infiltration of immune cell types were procured from Charoentong et al. (35). The content of immune cells in individual samples in the ssGSEA was estimated by utilizing differentially expressed marker genes. Each enrichment score was denoted by the relative content of each immune cell type. Furthermore, the Kruskal–Wallis test was used to analyze the distinctions in immune cell abundance between CR clusters to better comprehend the associations between CR clusters and immune cell infiltration in HCC.

Gene-annotation enrichment analyses were utilized to explore the differences in biological processes between distinct CR-related patterns through the package clusterProfiler in R (36). The gene sets of h.all.v7.5.1.symbols were procured from the Molecular Signatures Database v5.1 (MSigDB) (http://www.broad.mit.edu/gsea/msigdb/). Herein, a distinct energy metabolic scoring system was defined based on Dr. Yu et al.’s energy metabolism classifier for breast cancer (36) and the gene sets of HALLMARK_GLYCOLYSIS.v7.5.1, which was acquired from the MSigDB.

Additionally, 10 cuproptosis-related genes were retrieved from the literature and divided into activated and inhibitor groups ( Supplement Table 3 ) (26). Moreover, these analyses were run by implementing the ComplexHeatmap and gene set variation analysis (GSVA) (37) packages in R to quantify the heterogeneity in different CR-related HCC patterns.

A scoring system, named CRsscore, was constructed to assess the epigenetic regulation pattern of individuals with HCC as follows. DEGs were identified between CR clusters via the package limma in R. The significance criteria for determining DEGs were a P-value of<0.001 and |log2 fold change (FC)| of >2.0 (38). DEG intersection from different CR clusters in the TCGA-HCC cohort and the genes involved in the ICGC LIRI-JP cohort were regarded as the ultimate DEGs. The prognostic genes in the TCGA cohort were screened using univariate Cox regression analysis on the premise of correlating (P< 0.01) with the overall survival (OS) of patients. All 152 genes were further incorporated into a least absolute shrinkage and selection operator (LASSO) analysis for dimension reduction in the “glmnet” R package. Next, a multivariate Cox analysis further screened five genes based on the lowest Akaike information criterion value. The CRsscore of our model for each sample was determined by the relative expression of each CR phenotype-related gene and its associated Cox coefficient.

The ICGC LIRI-JP cohort was used to validate the prediction effect of the model.

In order to explore the difference in the responses to chemotherapeutic drugs between the two sets, the semi-inhibitory concentration (IC50) values of the chemotherapeutic drugs which are usually utilized to treat LIHC was calculated by using the “pRRophetic” package.

The tissue microarray (TMA), including 30 paired liver cancer tissues and para-carcinoma tissues, was derived from 2019 to 2021 at the First Affiliated Hospital of Wenzhou University. This study was supported by the hospital’s ethics committee, and all the patients provided informed consent.

TMA sections (4-μm thick) were deparaffinized and hydrated, and 0.3% hydrogen peroxide, and incubated with primary antibody overnight at 4°C and with secondary biotinylated goat anti-rabbit antibody successively; the sections were then stained using SignalStain^® DAB (Cell Signaling Technology, Danvers, MA) and counterstained with hematoxylin QS (Vector Laboratories). The intensity of staining (0, 1, 2, 3) and the proportion of positive cells (0%–100%) were semi-quantified, and scored from 0 (no stained cells) to 3 (all cells intensely stained). The detail information of tissue microarrays is available in Supplementary Table 1 .

Pearson correlation was used to analyze the correlations between variables, and a t-test was used to explore the continuous variables that conformed to a normal distribution between binary groups. The Kruskal–Wallis test was utilized to distinguish the differences for comparison of the three clusters. The cutoff values of each dataset were evaluated with the survival outcome and CRsscore in each dataset using the R package survminer. The Kaplan–Meier method was applied to depict survival curves for the subsets in each cohort, and the log-rank test was used to statistically identify significant differences. Significance was defined at P< 0.05 in the premise. All statistical analyses were performed using R, version 4.1.0.

A flow chart of the data processing and course in this study is presented in Figure 1 .

The expression matrix was obtained from the 370 samples in the TCGA cohort after data processing and quality assessment. A total of 549 DEGs (267 upregulated and 282 downregulated) were derived for subsequent analysis under the threshold of a P-value of<0.05 and |log2FC| of ≥0.2 ( Supplement Table 2 ).

The samples of the TCGA LIHC cohort were clustered using the average linkage and Pearson’s correlation methods ( Figure 2A ), and the co-expression network was developed by implementing the co-expression analysis. The power of β = 10 (scale-free R2 = 0.95) in this study was screened as the soft-thresholding parameter to guarantee a scale-free network ( Figures 2B, C ). In total, two modules were found by the average linkage hierarchical clustering. The turquoise module had the most significant relationship with the tumor stage ( Figures 2D, E ), and this module was screened as the crucial clinical module for further exploration. Finally, the association between different modules was illustrated by an eigengene adjacency heatmap ( Figure 2F ).

A total of 30 hub genes were calculated by applying the cytoHubba plugin based on the MCC algorithm (33)( Supplement Figure 1A ). The intersection between the 30 hub genes and the gene involved in GSE14520 was taken to apply this classifier to multiple datasets. Finally, an RF-related classifier involving 16 hub genes was customized ( Supplement Table 3 ).

Three unique CR clusters ( Figure 3A ) were identified by unsupervised clustering in the TCGA cohort according to 16 hub genes. Importantly, analysis from the ICGC LIRI-JP cohort as well as GSE14520 externally verified the stability of our clustering results ( Supplement Figures 2A, B ). The Kaplan–Meier analysis illustrated that cases in cluster 3 correlated with more adverse prognoses ( Figure 3B ). Similar results were obtained for OS in the ICGC LIRI-JP cohort and GSE14520 ( Figures 3C, D ). Hence, the classifier robustness was well validated. The cluster-related gene expression distribution and specific clinical characteristics between the subgroups were revealed ( Figure 3E ). Our CR-based classification revealed the CR-related gene expression levels, which were abundant in cluster 3. Similar results were obtained in the ICGC LIRI-JP cohort as well as GSE14520 ( Supplement Figures 3A, B ). These results draw the same conclusion that patients with HCC with poor outcomes are abundant with the 16 hub gene expressions.

The relative amount of infiltrating immune cell constitution in the TME of HCC between the CRsclusters was calculated via the ssGSEA algorithm to find the correlation between TME and CR-related subtypes. Moreover, a heatmap was utilized to visualize the significant differential levels of infiltrating immune cells, which were defined with a strict cutoff of p< 0.05 by the Kruskal–Wallis test ( Figure 4A ). CRscluster1 had features of high TME immune cell infiltration with conspicuous surges in the infiltration of natural killer cells, B cells, pDCs, Th2 cells, mast cells, DCs, and neutrophils as well as cytolytic activity, tumor necrosis factor (TNF) II interferon (IFN) response, TNF I IFN response, and CCR. Figure 4B reveals the surprising negative correlation between the 16 hub gene expression and most immune signatures. This is consistent with our conclusion.

The biological process among CR-related clusters was further explored via performing ssGSEA for hallmark gene sets. Figure 4A shows that CRscluster1 was strongly activated in stromal and metabolism pathways, such as glutathione, fatty acid, and phenylalanine metabolism pathways. CRscluster3 was markedly related to carcinogenic and DNA damage repair-associated pathways. CRscluster2 was the intermediate state of the other two clusters. These phenomena are similar to the GSEA analysis results ( Supplement Figure 3A ).

Subsequently, markers that represent DNA damage repair signaling pathways were screened ( Supplement Table 4 ) and determined among different clusters by ssGSEA. The phenomenon ( Figure 4C ) illustrated that CRscluster3 is associated with better DNA damage repair than the other CRsclusters. The ICGC LIRI-JP cohort and GSE14520 analysis drew the same conclusion ( Figures 4D, E ), considering the survival analysis results ( Figures 3C, D ).

The customized energy metabolic scoring system containing four central metabolic pathways, including glycolysis, glutaminolysis, fatty acid oxidation (FAO), and the pentose phosphate pathway (PPP) (39)( Figure 5A , Supplement Table 4 ), was used to further explore the metabolic heterogeneity among different clusters. The abundance of the activities of the four metabolic pathways was then evaluated by ssGSEA among different clusters. As expected, CRscluster1 was more dependent on FAO and glutaminolysis. CRscluster3 was more dependent on glycolysis and PPP. CRscluster2 was the intermediate state of the other two clusters ( Figure 6A ). Surprisingly, the ICGC LIRI-JP cohort and GSE14520 analysis externally verified the robustness of our results ( Figure 6B,C ). Moreover, the relationship between the four metabolic pathway activities and outcomes in this study followed breast cancer (39). We then analyzed the activity of cuproptosis-related genes by ssGSEA among different clusters. CRscluster1 was abundant with cuproptosis activity ( Figure 6D ). The ICGC-LIRI-JP cohort and GSE14520 show similar results in combination with the finding of the survival analysis ( Figures 6E, 6F ).

The enrichment score of 10 classical oncogenic pathways was evaluated via referred signatures ( Figure 6G , Supplement Table 4 ). Oncogenic pathways, such as hippo-related signaling and phosphatidylinositol 3-kinase (PI3K) signaling, had higher scores in CRscluster3. These results are consistent with those of previous studies associated with the glycolytic cancer tendency (39).

CRs exert a profound effect on shaping different TME landscapes, but the CR-related pattern in individuals cannot be conveniently predicted. Hence, we tried to develop a set of CRsscores to quantify the CR-related pattern of individuals with HCC. We first identified 902 DEGs across clusters 1–3 ( Supplement Table 5 ). Additionally, we applied a GO and KEGG analysis to explore the biological pathways associated with the DEGs. DEGs between diverse CR phenotype-related patterns were found to be enriched in metabolism and epigenetic-related biological processes ( Supplement Figures 4A, B ).

In order for the CRsscore to be well validation by other datasets, the final 429 DEG2, which intersected between the 902 DEGs and the gene involved in the ICGC LIRI-JP cohort, were further analyzed as the candidates ( Supplement Figure 5A ). The univariate Cox analysis selected 152 CR phenotype-related genes ( Supplement Table 6 ) that were incorporated in the LASSO and multivariate Cox analyses. Eventually, five genes (CDCA8, NEIL3, ANXA10, PON1, and CYP26B1) were identified as independent HCC prognosis indicators. Consequently, we developed the CRsscore using the following formula: CRsscore = (0.088144 ´ expression of CDCA8) + (0.18774´ expression of NEIL3) + (-0.023476 ´ expression of ANXA10) + (-0.002680 ´ expression of PON1) + (0.081717 ´ expression of CYP26B1).

Samples with HCC in the high CRsscore set had a shorter OS according to the K-M survival analysis (P< 0.001, Figure 7A ). Moreover, the area under the curve (AUC) of the CRsscore was 0.823, which shows a more accurate predictive ability than that of the traditional clinicopathological features ( Figure 7F ). The predictive value of the AUC of the CRsscore regarding the 1-, 2-, and 3-year survival rates was 0.823, 0.736, and 0.731, respectively ( Figure 7B ).

The hazard ratio and 95% CI of the CRsscore in the univariate (P< 0.001) and multivariate Cox regression analyses (P< 0.001) respectively elucidated the CRsscore as a cancer indicator ( Supplement Figures 5B, C ) and independent prognosis index of OS in patients with HCC.

The heatmap of the association between the clinicopathological features and the CRsscore is also presented ( Figure 7G ). A hybrid nomogram (c-index = 0.735) encompassing the CRsscore and clinicopathological features is shown in Supplement Figure 5D . The practical and predicted 1-, 2-, and 3-year survival rates following the reference curve via the calibration curve analysis are depicted in Supplement Figure 5F . These findings suggest that the nomogram was precise and steady; therefore, its implementation in the clinical services of patients with HCC is appropriate.

The CRsscore of each patient in the ICGC LIRI-JP cohort was also calculated, and the cohort was divided into two groups based on the median value. A survival analysis illustrated a better outcome in the low-risk group (log-rank test; p< 0.001; Figure 7D ). An analysis of the 1-, 2-, and 3-year prognostic prediction classification efficiencies suggested that the CRsscore still had relatively high AUC values ( Figure 7E ), indicating that the CRsscore had a prominent ability to predict HCC prognosis.

The energy metabolism level, cuproptosis level, and DNA damage repair pathway were further explored among CRsscore-high and -low groups by applying ssGSEA ( Figure 8A , Supplement Figure 4A ). CRsscore-high, which was associated with adverse outcomes, showed more abundant DNA damage repair pathway, glycolysis, and PPP. CRsscore-low, which was related to a benign prognosis, showed more abundant cuproptosis, FAO, and glutaminolysis. These results were identical to the results of the CR cluster analysis. Moreover, Supplement Figures 7A–E reveal that the levels of cuproptosis activity, FAO, and glutaminolysis in patients with HCC significantly decreased as the CRSsscore increased. In contrast, the glycolysis and PPP levels increased. Lastly, we further explored the correlation between the CRsscore and cuproptosis activity from a genomics perspective ( Supplement Figure 7F ). Surprisingly, the CRsscore had a significant negative correlation with FDX1 expression, which was the most crucial for cuproptosis regulation (26).

The IC50 values which can reflect the sensitivity to chemotherapeutic drugs of usual chemotherapeutic drugs were predicted and compared. Patients in the low-CRsscore group were more susceptible to Sorafenib and Gefitinib, whereas patients in the high-CRsscore set were more responsive to Cisplatinl, Mitomycin.C and Doxorubicin( Supplement Figures 8A–E ).

Immunohistochemical images of CDCA8, NEIL3, ANXA10, PON1, and CYP26B1 were obtained from the First Affiliated Hospital of Wenzhou University cohort( Figure 9 ). These results have indicated that higher expressions of PON1, ANXA10 genes in para cancerous tissues. Meanwhile, higher expressions of CDCA8, NEIL3 genes in liver cancer and the expression of CYP26B1 was no significant difference between para-cancer and cancer. Most important of all, there were many studies involved in the relative expression of PON1 (40–43), ANXA10 (44, 45), CDCA8 (46–48) and NEIL3 (49, 50) in HCC consisting of my experiment. This further revealed validity of the CRsscore.