The transcriptome profiles and the corresponding clinicopathological data of patients with liver cancer were obtained from TCGA database (https://portal.gdc.cancer.gov/) as the training cohort. In addition, RNA-sequencing (RNA-seq) data and clinical trait information from patients with liver cancer were downloaded from the ICGC database (LIRI-JP dataset; http://icgc.org/) and Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/) for validation. Specifically, 12,518 CRGs in the GSE117623 dataset were downloaded from the GEO database (24). Transcriptomic and matched clinical data from the IMvigor210 cohort of patients treated with anti-PD-L1 were collected (research-pub.gene.com/IMvigor210CoreBiologies) to explore the value of model genes in assessing response to immunotherapy (25).

To acquire the differentially expressed genes (DEGs) associated with CTCs/CTMs, the limma R package (version 2.7, bioinf.wehi.edu.au/limma) was used to process the RNA-seq data using a false discovery rate (FDR) <0.05 and |log2 (fold change)|>2 as the cutoff criteria. A Venn diagram was then constructed using a Venn webtool (http://bioinformatics.psb.ugent.be/webtools/Venn/) to illustrate the intersection among TCGA-DEGs, ICGC-DEGs and genes from the GSE117623 dataset. These intersected genes were considered to be the CRGs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted to explore the functional roles and pathways associated with the CRGs using the clusterProfiler R package (version 3.19) (26). The cut-offs were set as P<0.05 and FDR <0.05. Gene set enrichment analysis (GSEA) was performed to investigate the common biological pathways (27) using cp.kegg.v7.1.symbols.gmt as a reference gene set with a threshold of P<0.05, to screen for key enriched pathways in different risk groups. In addition, interactions among the CRGs were illustrated by the construction of a PPI network using the STRING database (https://string-db.org/), with an interaction score >0.7 being considered significant. Moreover, Cytoscape software was utilized to visually represent the PPI network. Specifically, the Cytoscape plug-in Molecular Complex Detection (MCODE) (version 2.0.3) was utilized to identify the highly interconnected modules of the PPI network with the following criteria: Degree cut-off, 2; node score cut-off, 0.2; k-core, 2; and max. depth, 100 (28). In addition, another Cytoscape plug-in, cytoHubba (version 0.1), was used to rank the nodes in the network according to their network functionality (29). The gene set variation analysis (GSVA) package (version 3.19) was used to explore the signaling pathways between high- and low-risk groups (30).

Univariate Cox regression analysis was performed to determine the prognostic CRGs and the CRGs associated with survival time, with P<0.01 considered to be statistically significant. Then, least absolute shrinkage and selection operator (LASSO) penalized Cox regression analysis was performed to further filter prognostic CRGs associated with the overall survival (OS) of patients with liver cancer (31). Subsequently, a risk signature was developed via stepwise multivariate Cox proportional hazards regression analysis. Prognostic gene signatures were constructed based on linear combinations of regression coefficients derived by multiplying the LASSO Cox regression model coefficients by their mRNA expression levels (32): Risk score=∑ (βmRNA × mRNA)n, where β represents the regression coefficient for the mRNA, mRNA represents the expression level of the mRNA, and n represents the specific gene. Receiver operating characteristic (ROC) curves and Kaplan-Meier curves were constructed to evaluate the predictive performance of the prognostic model in TCGA cohort. pheatmap R package (version 1.0.12; cran.r-project.org/web/packages/pheatmap/index.html) to plot images describing gene expression heatmaps, risk scores and OS for high and low risk groups. Data from the ICGC and GEO databases were used as external validation data to test the predictive capability of the model.

A non-negative matrix factorization (NMF) clustering algorithm was utilized to analyze the five signature genes in the risk score model, and determine the subtypes of CRGs in liver cancer using the NMF R package (version 0.27) (33). Using conformal, scatter and silhouette features, the optimal number of clusters with n=2 was determined.

Nomograms are widely used as a tools for the prognostic analysis of patients with tumors (34). A simplified liver cancer nomogram was constructed for each dataset based on the CRG model and its predictive performance was evaluated by plotting calibration curves.

The association between the low- and high-risk groups and clinical characteristics were explored using Chi-square tests, and the results were displayed as a heatmap. In addition, the associations between the signature genes and immune cell infiltration were analyzed. Six algorithms, namely CIBERSORT-ABS (35), TIMER (36) (https://cistrome.shinyapps.io/timer/), QUANTISEQ (37), MCPCOUNTER (38), XCELL (39) and EPIC (40,41), were used to evaluate the differences in the immune microenvironment between the two risk groups. Tumor-associated immune comprehensive score was assessed via ImmunoPhenoScore in R package IOBR (42) (version 0.99.9, http://github.com/IOBR/IOBR). Waterfall plots for the two risk groups were produced using the maftools (github.com/PoisonAlien/maftools) R package (version 3.19). Differences in the expression of major histocompatibility complex (MHC) molecules, human leukocyte antigen (HLA) signature, chemokines and potential immune checkpoints were also compared between the two groups. To investigate the association between signature genes and immune subtypes, ‘Subtypes’ module of the TISIDB database (http://cis.hku.hk/TISIDB/index.php). Pearson correlation coefficients of the signature genes expression with the immune checkpoints (PD-1, PD-L1 and CTLA4) were calculated using R language to assess the correlation. In this study, the OCLR algorithm and the Primary Cell Biology Consortium (PCBC, http://progenitorcells.org/) stemness score model were used to calculate the mRNAsi of cells in the TCGA-LIHC dataset and to assess the correlation between the stemness index and the risk score (43). Pearson correlation coefficients of the signature genes expression with the immune checkpoints were calculated using the R language to assess the correlation.

To identify small molecule compounds that may be suitable for the treatment of liver cancer, the pRRophetic (genemed.uchicago.edu/~pgeeleher/pRRophetic/) R package (version 3.19) was used to calculate the half-maximal inhibitory concentration (IC₅₀) based on data from the Genomics of Drug Sensitivity in Cancer database (44).

HepG2 and MHCC97H human liver cancer cells (cat. nos. CTCC-001-0014 and CTCC-400-0192, respectively) were obtained from the Meisen Chinese Tissue Culture Collections. The cell lines were authenticated by short tandem repeat testing. Both cell lines were cultivated in high-glucose Dulbecco's modified Eagle's medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) at 37°C in 5% CO₂. Two small interfering RNA (siRNAs) targeting tumor protein p53 inducible protein 3 (TP53I3), namely si-TP53I3-1 and si-TP53I3-2, and an siRNA negative control were synthesized by and purchased from Sangon Biotech Co., Ltd. The sequences of siRNAs are listed in Table SI. Cell transfection was conducted in 6-well plates when cell confluence was 60–70%, with a final siRNA concentration of 50 nM per well. Transfection of the liver cancer cells was performed using Lipofectamine^® 2000 reagent (Invitrogen; Thermo Fisher Scientific, Inc.) following the manufacturer's instructions. Transfection was performed for 6–8 h at 37°C in 5% CO₂. The cells were harvested at 24 h post-transfection for reverse transcription-quantitative PCR (RT-qPCR) analysis and at 48 h post-transfection for western blot and in vitro functional assessment.

Cells were lysed on ice with RIPA buffer (Wuhan Boster Biological Technology, Ltd.) containing protease inhibitor cocktail (MedChemExpress) for 20 min. The protein contents of the cell lysates were quantified using a BCA protein assay kit (Beyotime Institute of Biotechnology). Then, 30 µg protein/lane was separated by 10% SDS-PAGE (Boster Biological Technology) and transferred to PVDF membranes (EMD Millipore). The membranes were blocked with 5% defatted milk at room temperature for 2 h, then incubated with anti-TP53I3 (#14828-1-AP; 1:1,000; Proteintech Group, Inc.) and anti-ACTB (#AC006; 1:3,000; ABclonal Biotech Co., Ltd.). primary antibodies at 4°C for 12–16 h, followed by HRP-conjugated Affinipure goat anti-rabbit IgG (H+L) (SA00001-2; 1:5,000; Proteintech) secondary antibodies at room temperature for 2 h, and the signal was detected using Pierce^® ECL Western Blotting Substrate (Thermo Fisher Scientific, Inc.). Finally, the bands were detected and analyzed using ChemiDoc^™ XRS+ with Image Lab^™ software (version 6.0, Bio-Rad Laboratories, Inc.).

Total RNA was extracted from cells using FreeZol reagent (Vazyme Biotech Co., Ltd.) and synthesized into cDNA using PrimeScript^™ RT Master Mix (Takara Bio, Inc.), according to the manufacturer's instructions. qPCR was then carried out using the CFX96 Real-Time PCR System (Bio-Rad Laboratories, Inc.) with the SYBR Green PCR kit (Thermo Fisher Scientific, Inc.) according to the standard protocol. The thermocycling conditions used were as follows: 95°C for 30 sec pre-cycling, and then 40 cycles of 95°C for 10 sec and 60°C for 30 sec. The primer pairs were synthesized by Sangon Biotech Co., Ltd. and their sequences are presented in Table SII. The relative expression of TP53I3 was calculated using the formula 2^−ΔΔCq with GAPDH as the reference gene (34).

Cell Counting Kit 8 (CCK-8) assay (ABclonal Biotech Co., Ltd.) was utilized to assess the proliferation ability of the cells. Cells (3,000/well) were plated in a 96-well plate and incubated overnight at 37°C to allow adhesion. At 24, 48 and 74 h, 100 µl 10% CCK-8 solution was added to each well and the cells were cultured in a cell incubator for 2 h, after which absorbance was measured at 450 nm using a microplate reader (Thermo Fisher Scientific, Inc.).

Cells were seeded into 6-well plates at a concentration of 1,000 cells/well. The cells were cultured at 37°C with 5% CO₂ in fresh medium and allowed to grow for 14 days. The colonies were then fixed for 15 min at room temperature in 4% paraformaldehyde (Wuhan Servicebio Technology Co., Ltd.), and stained with crystal violet (0.5% wt./vol.) at room temperature for 15 min. Finally photographs of the plates were taken and the colonies were quantified using an inverted microscope (Guangzhou Micro-shot Technology Co., Ltd.). The number of colonies was counted manually. Each independently counted colony refers to a cell cluster of ≥50 cells. The experiment was repeated three times.

The BeyoClick^™ EdU-555 Cell Proliferation Kit (Beyotime Institute of Biotechnology) was employed to investigate the proliferation rate of the human liver cancer cells according to the manufacturer's protocols. Briefly, after incubation with 1X EdU (10 µM) solution for 2 h at 37°C, cells were fixed with paraformaldehyde (4%) for 30 min at room temperature, then permeabilized with 0.3% Triton X-100 for 15 min and finally stained with Hoechst 33342 and 4′,6-diamidino-2-phenylindole in the absence of light for 30 min at room temperature. Finally, the cells were imaged by fluorescence microscopy.

Bioinformatics analysis and mapping were accomplished using R software. Survival rates were compared using Kaplan-Meier analysis with the calculation of P-values using log-rank tests, or the 2-stage test in the plot with late-stage crossover (cran.r-project.org/web/packages/TSHRC/TSHRC.pdf). In addition, the Chi-square test was used for comparisons between categorical variables, and unpaired Student's t-test was utilized to evaluate the discrepancies between the two risk groups. Correlations between variables were assessed using Spearman's correlation test. The cell groups were compared by one-way ANOVA followed by Dunnett's post hoc tests. For each statistical analysis, P<0.05 was considered to indicate a statistically significant result.

The liver cancer (liver hepatocellular carcinoma) gene expression profiles were downloaded from TCGA and ICGC portals and 1,622 and 628 DEGs, respectively, were screened out using the limma R package. The DEGs from TCGA and ICGC databases are shown as volcano plots in Fig. 1A and B, respectively. A Venn diagram was then constructed to filter out the differentially expressed CRGs (Fig. 1C). The intersection of the DEGs from TCGA and ICGC databases with the 12,518 CRGs from GSE117623 yielded a total of 258 differentially expressed CRGs (Table SIII).

To investigate the interrelationship of the differentially expressed CRGs and identify hub genes, a PPI network was constructed and module analysis performed to determine co-expression networks. Firstly, the 258 differentially expressed CRGs were uploaded to the STRING database, and the minimum required interaction score was set to 0.7, which indicates a strong interaction between the CRGs. The STRING interactions were then analyzed using Cytoscape and the resulting co-expression network, which contained 155 nodes and 2,731 edges, is shown in Fig. 2A. In addition, modules with >50 genes were identified using the MCODE plug-in and 10 hub genes in that module, namely topoisomerase IIa, cyclin B2, cell division cycle associated 8 (CDCA8), BIRC5, aurora kinase B, cyclin B1, BUB1 mitotic checkpoint serine/threonine kinase (BUB1), BUB1B, kinesin family member 20A and TTK protein kinase, were characterized using the cytoHubba plug-in (Fig. 2B). This included 57 nodes and 1,497 edges. These potential hub genes may be instrumental in the biological progression of liver cancer.

To explore the biological categories and biological processes associated with the differentially expressed CRGs, GO and KEGG enrichment analyses were conducted using R software, and the enrichment results are shown in bubble charts (Fig. 2C and D). The GO enrichment analysis revealed that the differentially expressed CRGs were principally concentrated in the biological process terms ‘nuclear division’, ‘organelle fission’, ‘chromosome segregation’ and ‘mitotic nuclear division’. In addition, the main cellular component terms included ‘chromosomal region’, ‘chromosome, centromeric region’, ‘spindle’ and ‘kinetochore’. Moreover, the molecular function terms associated with the CRGs were ‘tubulin binding’, ‘microtubule-binding’, ‘ATPase activity’ and ‘motor activity’ (Fig. 2C). Regarding the KEGG analysis, the primary terms are shown in Fig. 2D, which reveals that the differentially expressed CRGs were particularly enriched in ‘cell cycle’, ‘microRNAs in cancer’, ‘p53 signaling pathway’ and ‘cellular senescence’.

Univariate Cox regression analysis demonstrated that 88 CRGs were strongly associated with survival in patients with liver cancer (P<0.01), all of which were prognostic risk factors (Fig. 3A). Then, the 88 CRGs were regression penalized using LASSO Cox regression to exclude relatively insignificant parameters (Fig. 3B and C). Stepwise multivariate Cox regression was subsequently employed to construct a predictive signature for patients with liver cancer in TCGA cohort (Fig. 3D). The five genes in the signature were CDCA8, TP53I3, hepatitis A virus cellular receptor 1 (HAVCR1), MYCN proto-oncogene (MYCN) and thioredoxin reductase 1 (TXNRD1). The formula for risk score calculation was as follows: Risk score=(0.0826 × expression level of CDCA8) + (0.0112 × expression level of TP53I3) + (0.0824 × expression level of MYCN) + (0.0376 × expression level of HAVCR1) + (0.0120 × expression level of TXNRD1) (Table SIV). Patients in TCGA cohort were classified into high- and low-risk groups using the median predictive index as the cut-off point. As Fig. 4A shows, the low-risk group was significantly associated with improved survival (P<0.05). To evaluate the predictive ability of this prognostic signature, an ROC analysis of the risk score was conducted. The area under the curve (AUC) values predicted from the ROC curves for 1-, 3- and 5-year OS were 0.804, 0.736 and 0.707, respectively (Fig. 4B). In Fig. 4C, the upper panel shows the expression heat map of the five prognostic model genes in the high- and low subgroups, the middle panel reveals that the risk of patients with liver cancer increases as risk score increases, and the lower panel demonstrates the poor OS of the patients in the high-risk group compared with those in the low-risk group. To validate the predictive power of the signature, the same formula was used to analyze the risk score of each patient in the ICGC dataset, for independent external validation. The Kaplan-Meier curves also displayed a poor prognosis of patients in the high-risk group in this dataset (P<0.05; Fig. 4D). In addition, the ROC curve showed the strong predictive ability of the risk-score signature for prognosis, with AUCs for the prediction of 1-, 3- and 5-year OS of 0.714, 0.730 and 0.726, respectively (Fig. 4E). Also, the expression of the five CRGs and the mortality of the patients increased as the risk scores increased (Fig. 4F).

To investigate whether the CTC/CTM-associated 5-gene signature could be an independent prognostic factor for patients with liver cancer, the prognostic value of this signature was compared with that of several clinicopathological factors, including age, sex, grade and American Joint Committee on Cancer (AJCC) stage in both cohorts using univariate and multivariate Cox regression analyses. For TCGA cohort, 365 valid patients were included, 182 in the high-risk group and 183 in the low-risk group. Univariate Cox analysis indicated that risk score [P<0.001; hazard ratio (HR), 1.258; 95% confidence interval (95% CI), 1.183–1.338)], AJCC stage (P<0.001; HR, 1.658; 95% CI, 1.340–2.053) and T status (P<0.001; HR, 1.633; 95% CI, 1.332–2.003) were candidate factors. Further multivariate Cox regression analysis emphasized that risk score was an independent risk factor for patients with liver cancer (P<0.001; HR, 1.227; 95%CI, 1.142–1.318) (Fig. 4G). For the ICGC cohort, 229 patients were included, 59 in the high-risk group and 170 in the low-risk group. Univariate Cox regression analysis demonstrated that sex (P=0.033; HR, 0.505; 95% CI, 0.270–0.946), risk score (P=0.010; HR, 1.076; 95% CI, 1.017–1.137) and tumor stage (P<0.001; HR, 2.252; 95% CI, 1.543–3.287) were potential risk factors. Multivariate Cox regression analysis also indicated that sex (P=0.004; HR, 0.389; 95% CI, 0.204–0.740), risk score (P=0.030; HR, 1.070; 95% CI, 1.006–1.138 and tumor stage (P<0.001; HR, 2.330; 95% CI, 1.601–3.391) were independent predictors for patients in the ICGC cohort (Fig. 4H). In conclusion, these findings indicated that the 5-CRG risk signature was closely associated with the clinical characteristics of patients with liver cancer, had a fine predictive capacity and has potential as a prognostic indicator for these patients.

To further verify the predictive power of the prognostic signature, the GSE14520 dataset was analyzed. In this dataset, 221 valid patients were included, with 100 in the high-risk group and 121 in the low-risk group. As the results in Fig. 5A illustrate, patients in the low-risk group had improved survival outcomes compared with those in the high-risk group (P<0.05). In the GEO cohort, the AUC for 5-year OS was 0.684 (Fig. 5B). The expression of model genes, risk score distribution and survival status for each patient in this validation cohort are shown in Fig. 5C. Following univariate Cox regression analysis (Fig. 5D), the results of independent prognostic analysis revealed that risk score (P<0.001; HR=1.280; 95% CI, 1.145–1.432), as well as AJCC stage (P<0.001; HR, 2.017; 95% CI, 1.518–2.681) (Fig. 5E) were independent risk factors in this cohort. These findings indicate that the prognostic model is promising as a predictive signature.

Patients were clustered into different subtypes based on the expression levels of the five prognostic signature genes using the NMF algorithm. To ensure the robustness of the clustering results, the coefficient of correlation was used to determine the optimal number of clusters, and when the number of clusters was 2, clear boundaries were observed for both subtypes. This indicated the stable and reliable clustering of the liver cancer samples (Fig. 6A). The OS of patients in cluster 1 (C1) was significantly improved compared with that of C2 (P=0.002; Fig. 6B). Most immune checkpoints were upregulated in the C2 group compared with the C1 group (Fig. 6C). In addition, it was also found that the level of immune infiltration in the tumor microenvironment was also distinct in the two groups, with immune score, stromal score and ESTIMATE score of the C1 group being significantly lower than those of the C2 group (P<0.05; Fig. 6D). This suggests that C1 molecular subtype tends to present ‘cold tumors’, whereas the C2 molecular subtype tends to present ‘hot tumors’. It was also noted that the level of neutrophil infiltration was higher in the C1 group than in the C2 group, whereas the levels of monocytic lineage and fibroblast infiltration were lower in the C1 group than in the C2 group (P<0.05; Fig. 6E).

The sensitivity of the high- and low-risk groups to various chemotherapeutic agents was compared to evaluate drugs for potential use in liver cancer. The findings indicate that the low-risk group was associated with a higher IC₅₀ for chemotherapeutic compounds including ABT.888 (veliparib), AS601245 (an ATP-competitive JNK inhibitor), AG.014699 (rucaparib), A.443654 (a pan-Akt inhibitor), ATRA (tretinoin) and AUY922 (luminespib). By contrast, axitinib, A.770041 (an LCK inhibitor), AZD.0530 (saracatinib), AMG.706 (motesanib), AKT.inhibitor.VIII and AICAR (acadesine) had a higher IC₅₀ in the high-risk group, indicating that patients in the low-risk group may benefit more from treatment with these compounds (P<0.05; Fig. 7). The sensitivity of the two liver cancer subtypes to various chemotherapeutic drugs was also evaluated. The results suggested that patients with the C1 subtype might be more sensitive to metformin, lapatinib, elesclomol, docetaxel, camptothecin, bosutinib, axitinib and vinblastine, while patients in group C2 would likely benefit by treatment with cisplatin, bortezomib, bleomycin, bicalutamide, mitomycin C, imatinib, etoposide and gemcitabine (Fig. S1).

To further explore the molecular mechanism associated with the signature genes and the prognostic module, GSEA was performed in TCGA liver cancer cohort. Fig. 8A-E reveals the KEGG pathways of the five signature genes, namely CDCA8, HAVCR1, MYCN, TP53I3 and TXNRD1, showing the five most upregulated and downregulated pathways for each gene. The signature genes are mainly concentrated in KEGG pathways including ‘cell cycle’, ‘p53 signaling pathway’, ‘complement and coagulation cascades’ and ‘drug metabolism cytochrome p450’. In addition, GSEA was used to compare the high- and low-risk groups based on the risk scores. The KEGG pathways enriched in the high and low risk groups are shown in Fig. 8F.

GSVA was also utilized to analyze the differences in biological behavior between the high- and low-risk groups. The results demonstrated that pathways associated with tumor progression, such as ‘cell cycle’, ‘DNA replication’, ‘RNA degradation’, ‘mTOR signaling pathway’ and ‘P53 signaling pathway’, were mainly concentrated in the high-risk group. By contrast, metabolism-related pathways, including ‘fatty acid metabolism’, ‘propanoate metabolism’, ‘butanoate metabolism’ and ‘tyrosine metabolism’, were mainly present in the low-risk group of patients (Fig. 8G).

In view of the important role of immune checkpoints in tumor immunotherapy, the differential expression of immune checkpoint genes was analyzed between risk subgroups. The results revealed that common immune checkpoint genes, including cytotoxic T-lymphocyte associated protein 4 (CTLA4), CD274, programmed cell death 1, and T-cell immunoreceptor with Ig and ITIM domains were upregulated in the high-risk group compared with the low-risk group (P<0.05; Fig. 9A). This suggests that the poor prognosis of high-risk patients with liver cancer may at least partially be attributed to an immunosuppressive microenvironment. Chemokines and their receptors are necessary for the targeted migration of immune cells and the initiation and execution of the immune response (45,46). Therefore, the differential expression of chemokines and their receptors was analyzed in the two risk subgroups, which revealed higher levels of expression for the majority of these chemokines and receptors in patients in the high-risk group (P<0.05; Fig. 9B and C). An association between risk score and HLA-associated gene expression was also observed. As shown in Fig. 9D, the abundance of HLA-related genes was higher in patients at high risk than those in the low-risk group (P<0.05). The results of algorithms were visualized using heat maps, including assessment of immune cell infiltration in the two risk subgroups, and the results suggest that the high-risk group has more abundant immune cell infiltration (Fig. 9E). In addition, further exploration of the association between risk score and immune pathway activity revealed that cytolytic activity, type I IFN response and type II IFN response scores were higher in the low-risk group, and conversely, the MHC class I score was higher in the high-risk group (P<0.05; Fig. 10F). These results demonstrate that patients in the high-risk group are more likely to benefit from immunotherapy.

To investigate the prognostic value of the CRG signature in patients with different clinical features, a heat map was drawn to reveal whether there was a potential association with clinicopathological features in the high- and low-risk subgroups (Fig. 10A). The expression levels of CDCA8 and TXNRD1 were higher in the high-risk group than in the low-risk group. In addition, the results revealed that the high-risk score was closely associated with a higher T stage (P<0.001), higher grade (P<0.001), higher tumor stage (P<0.001) and poor patient survival status (P<0.001).

Since TMB and tumor immune dysfunction and exclusion (TIDE) are good indicators of the response to immunotherapy, sample scores were calculated for each patient with liver cancer and variability between the high- and low-risk subgroups was assessed. The results revealed that the high-risk group had a higher TMB and lower TIDE index, which further demonstrates that patients in the high-risk group should be more responsive to immunotherapy (P<0.05; Fig. 10B and C). The mRNA expression-based stemness score revealed a positive correlation between liver cancer tumor stemness and the risk score, indicating that tumors in the high-risk group are more likely to undergo malignant progression and thus lose their differentiated phenotype (P<0.05; Fig. 10D). In addition, the maftools R package was used to visualize the differences in somatic mutation distribution between the high- and low-risk groups. The results demonstrated that the high-risk group had a higher mutation frequency compared with the low-risk group (91.38 vs. 77.84%, respectively). The most mutated gene in the low-risk group was catenin b1 (23%) and the most mutated gene in the high-risk group was TP53 (43%) (Fig. 10E).

Immune checkpoints have an important role in immune regulation, and immune checkpoint inhibitors are used in cancer therapy. Therefore, the associations between the signature genes and the expression of immune checkpoint genes, namely programmed cell death protein 1 (PD-1), programmed death-ligand 1 (PD-L1) and CTLA4, were investigated. The results in Fig. 11A indicate that the expression of CDCA8 was positively correlated with that of the three immune checkpoints, PD-1 (R=0.3; P=4.9×10⁻⁹), PD-L1 (R=0.32; P=3.8×10⁻¹⁰) and CTLA4 (R=0.32; P=3.9×10⁻¹⁰). In addition, the expression of HAVCR1 was positively correlated with CTLA4 expression (R=0.32; P=1.3×10⁻¹⁰). TISIDB portal was used to analyze the expression of signature genes in different immune subtypes, specifically: C1, wound healing; C2, IFN-g dominant; C3, inflammatory; C4, lymphocyte depleted; C5, immunologically quiet; and C6, TGF-b dominant (47). The results indicated that the roles of these five genes differ among the different immune subtypes, with CACA8, MYCN and TXNRD1 being differentially expressed among the immune subtypes. Specifically, the TISIDB analysis revealed that CDCA8 was highly expressed in the C1 and C2 types, MYCN was highly expressed in the C1 type, and TXNRD1 was mainly expressed in the C2 and C4 types (Fig. 11B). In addition, the IMvigor dataset was used to predict the responsiveness of the five signature genes to atelelizumab treatment. Notably, consistent with the previous findings, the analysis suggested that patients with high expression of CDCA8 and TXNRD1 may obtain improved treatment outcomes (Fig. 11C). The correlations between tumor immune infiltration by CD4⁺ T cells, CD8⁺ T cells, B cells, neutrophils, macrophages and dendritic cells, and the expression of the five signature genes were also investigated (Fig. S2). In this analysis, correlation coefficients >0.3 and P<0.05 were considered as distinctive; partial.cor denotes partial correlation, indicating the correlation of gene expression with immune cell infiltration in the TIMER database. The results show that CDCA8 expression is positively correlated with the infiltration of six types of immune cells: B cells (partial.cor, 0.441; P=9.08×10⁻¹⁸), CD8⁺ T cells (partial.cor, 0.303; P=1.03×10⁻⁸), CD4⁺ T cells (partial.cor, 0.359; P=6.74×10⁻¹²), macrophages (partial.cor, 0.439; P=1.70×10⁻¹⁷), neutrophils (partial.cor, 0.368; P=1.63×10⁻¹²) and dendritic cells (partial.cor, 0.465; P=1.22×10⁻¹⁹). Similarly, HAVCR1 expression was found to be positively correlated with the infiltration of B cells (partial.cor, 0.302; P=1.14×10⁻⁸), macrophages (partial.cor, 0.302; P=1.34×10⁻⁸), neutrophils (partial.cor, 0.392; P=4.00×10⁻¹⁴) and dendritic cells (partial.cor, 0.317; P=2.18×10⁻⁹), and TXNRD1 expression was positively associated with neutrophil infiltration (partial.cor, 0.322; P=8.67×10⁻¹⁰).

To better assess the predictive efficacy of the CRG prognostic model, the risk signature was compared with six published liver cancer prognostic models. The signature of Du et al (48) was a m6A-based gene signature; the signature of Fu and Song (49) was a pyroptosis-related gene signature; the signature of Guo et al (50) was a signature containing nine genes; Lei et al (51) devised a starvation-based nine-mRNA signature; Tian et al (52) proposed a five-gene prognostic signature for liver cancer; and the signature of Zheng et al (53) comprised five pyroptosis-related genes. When the accuracy of these models and the current model were compared, it was found that the C-index and restricted mean survival of the CRG signature were higher than those of the other six models, which indicates that the present model is optimal (Fig. 12A and B). Additionally, the AUCs of the CRG model for 1-, 3- and 5-year OS were 0.807, 0.711 and 0.667, respectively, which were higher than those of the other signatures, which validates the previous results (Fig. 12C).

To forecast the survivability of patients with liver cancer, a nomogram including factors such as age, sex, stage and risk score was created to predict probability of OS at 1, 3, and 5 years in the TCGA cohort. In addition, calibration plots were constructed to evaluate the predictive power of the nomogram (Fig. 12D). Similarly, two nomograms were also constructed for the ICGC and GSE14520 cohorts (Fig. S3). These all indicate the good predictive power of the model.

Among the five signature genes, CDCA8, MYCN, HAVCR1 and TXNRD1 have previously been demonstrated to have a biological regulatory function in liver cancer (vide infra), but TP53I3 has been poorly studied in liver cancer. Therefore, the role of TP53I3 in liver cancer cells was evaluated using cellular experiments. TP53I3 was knocked down in HepG2 and MHCC97H cells using siRNA, and the transfection efficiency was verified by western blot analysis and RT-qPCR (Fig. 13A and B). To explore the impact of TP53I3 on the proliferation of liver cancer cells in vitro, CCK-8, EdU and colony formation analyses were performed. The results showed that the proliferation ability and colony formation of the liver cancer cells was significantly suppressed after TP53I3 depletion (P<0.05; Fig. 13C-E), which indicates that TP53I3 promotes the proliferation of liver cancer cells.