The present study incorporated one scRNA-seq dataset and two bulk transcriptomic profiles from human HCC specimens. From the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), the GSE242889 dataset was retrieved, which included single-cell transcriptomes from five patients with HCC with corresponding peritumoral controls (15). Among the five patients, three were men and two were women, with a median age of 55 years (range, 42–60 years). Bulk transcriptomic profiles with corresponding clinical data were acquired from two independent sources: The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset (https://portal.gdc.cancer.gov/) and the GEO dataset, GSE76427 (16). For analytical purposes, the TCGA-LIHC cohort served as the training cohort, whereas the GSE76427 cohort was used for subsequent validation. Subsequently, data cleanup was performed. In the TCGA-HCC dataset, only the samples in which the histological type was HCC and the new event type was locoregional recurrence, extrahepatic recurrence or intrahepatic recurrence were included. In addition, for these two bulk RNA-seq datasets, cases missing either RFS data (follow-up time/status) or RNA-seq profiles data were excluded, even when clinical annotations were available. Finally, a total of 342 patients with HCC and 50 normal controls from the TCGA-LIHC cohort were included. The cohort comprised 230 men and 112 women, with a median age of 61 years (range, 16–90 years). Additionally, 108 patients with HCC from the GSE76427 cohort were analyzed, consisting of 86 men and 22 women, with a median age of 64 years (range, 14–93 years).

The single-cell transcriptome dataset, GSE242889, was merged using the Harmony algorithm (17). The samples were then analyzed as described in our previous study (18), with modifications. First, quality control measures excluded cells demonstrating either insufficient sequencing depth (<200 detected genes) or potential multiplets (>4,000 genes), along with those showing elevated mitochondrial transcript proportions (≥50%). Second, cellular clustering was performed using the FindClusters algorithm in Seurat (resolution parameter, 0.8) (19), which incorporates the first 17 principal components (PCs) as inputs.

After the cells in the GSE242889 dataset were classified and annotated, hepatocyte data were extracted for subsequent analysis. The FindMarkers algorithm in Seurat was employed to identify DEGs in hepatocytes derived from different origin samples (19). The statistical thresholds for DEGs were as follows: |log₂fold-change (FC)|>0.5 and adjusted P<0.05. In the TCGA-LIHC dataset, differential expression analysis comparing tumor and normal samples was performed using the ‘limma’ R package (version 3.44.3) (20), and genes with |log₂ FC|>1 and adjusted P<0.05 were regarded as DEGs. Subsequently, the associations between the gene expression profiles and HCC recurrence outcomes were assessed using the ‘survfit’ function in the ‘survival’ package (https://github.com/therneau/survival; version 3.2.7). The log-rank test result was interpreted with P<0.05 indicating statistical significance. The overlapping genes of these results were defined as the DEGs associated with RFS.

Survival-related DEGs were utilized as candidate biomarkers to construct a prognostic signature. Least absolute shrinkage and selection operator (LASSO) Cox penalized regression analysis was performed using the R package ‘glmnet’ (version 4.1) (21), and genes demonstrating non-zero coefficients were incorporated into the final risk score calculation. The risk score for each patient with HCC was calculated as follows: Risk score=Σ (βi × Expi), where βi represents the non-zero coefficients from LASSO Cox regression and Expi denotes the expression levels of selected genes. Subsequently, patients were dichotomized into high- and low-risk groups based on median risk score stratification, followed by Kaplan-Meier survival analysis.

To identify the underlying mechanisms that affect HCC, the DEGs between the high- and low-risk subgroups were obtained. Functional enrichment analysis was then performed using the ‘clusterProfiler’ package (version 4.7.1) to identify significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms among these DEGs (22). Gene Set Enrichment Analysis (GSEA) was performed for the genes in the prognostic signature to further elucidate their biological functions. Additionally, the relative abundances of 22 immune cell subtypes within HCC samples were quantified using the CIBERSORT algorithm (23), and the immune-related score was calculated using the ‘estimate’ R package (version 1.0.13) (24).

To isolate independent recurrence predictors, univariate Cox regression was first applied, followed by multivariate analysis for 155 patients with HCC with complete clinical variables. These variables included age, sex, body mass index, tumor (T) stage, lymph node (N) stage, metastasis (M) stage, grade, pathological stage, α-fetoprotein level, Child-Pugh score and the risk score, which were classified using the appropriate guidelines (25–27).

Following prognostic signature development, external validation was performed using the independent GSE76427 cohort. The prognostic ability of the signature was first assessed using time-dependent receiver operating characteristic (ROC) curve evaluation at 6-, 12- and 24-month intervals. Patients were subsequently dichotomized into high- and low-risk subgroups using the median risk score as the cutoff threshold, followed by Kaplan-Meier survival curve analysis.

The present study was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University (Nanchang, China; approval no. 2024–078). Tissue samples were collected at the Second Affiliated Hospital of Nanchang University (Nanchang, China) from February 2024 to June 2024. Patients that were diagnosed with HCC through pathological examination were included in the present study, while any patients with coexisting tumors were excluded from the study. A total of five HCC tissue samples and five adjacent tissue samples were obtained following documented informed consent from all study participants. Among the patients, there were four men and one woman, with a median age of 51 years (range, 38–65 years).

RNA was extracted and reverse transcribed from the HCC and normal tissues. For qPCR, the primers used are presented in Table I. As the immunoglobulin λ constant 2 (IGLC2) cDNA sequence was not available on the National Center for Biotechnology Information website, it was retrieved and aligned from Ensembl (http://ensemblgenomes.org/). The RT process proceeded with an initial primer annealing step at 25°C for 5 min, followed by the main cDNA synthesis reaction at 42°C for 30 min, which was completed by a brief 5-sec heat inactivation step at 85°C. The qPCR thermocycling conditions were as follows: Initial denaturation at 95°C for 30 sec. This was followed by 40 cycles of denaturation at 95°C for 15 sec and annealing/extension at 60°C for 30 sec. Finally, a melting curve analysis was conducted by heating the samples from 65 to 95°C. To ensure reliability, each experimental procedure was performed in triplicate, and the study was independently performed three times. Gene expression levels were quantified using RT-qPCR with a PrimeScript^™ RT Master Mix (cat. no. RR036A; Takara Bio, Inc.) on an ABI 7500 Real-time PCR system detection instrument (Applied Biosystems; Thermo Fisher Scientific, Inc.), with relative expression determined using the 2^−ΔΔCq method (28).

R software (https://www.r-project.org/; version 4.2.3) was utilized for statistical analysis. The data are presented as mean ± SD from experiments conducted more than three times. Differences between two groups were analyzed using a two-sided unpaired Student's t-test. The Kaplan-Meier method was employed to plot survival curves, and differences were assessed using the log-rank test. The Cox proportional hazards model was used to determine independent factors. P<0.05 was considered to indicate a statistically significant difference.

From the original 60,284 cells in the GSE242889 dataset, the quality control pipeline retained 57,219 cells that passed all quality thresholds for downstream processing (Fig. 1A-D). The single-cell data were subsequently segregated into 30 distinct clusters, from which 24,630 characteristic marker genes were detected (Fig. 1E; Table SI). The 10 most significant markers per cluster were visualized using a heatmap (Fig. 1F). The results revealed that clusters 0, 9, 10, 11, 17, 18, 22, 26 and 28, containing 20,317 cells, were hepatocytes; clusters 3, 6, 12, 19, 20 and 25, containing 10,299 cells, were monocytes; clusters 1, 7 and 21, containing 6,824 cells, were T cells; clusters 2 and 24, containing 8,693 cells, were macrophages; clusters 5, 16, 23 and 29, containing 4,667 cells, were B cells; clusters 8, 13 and 27, containing 3,749 cells, were endothelial cells; cluster 14, containing 1,299 cells, was tissue stem cells; cluster 15, containing 1,091 cells, was neutrophils; and cluster 24, containing 280 cells, was common myeloid progenitors (Fig. 1G).

Among the 20,317 hepatocytes in the GSE242889 dataset, 17,560 cells were from HCC samples, and 2,757 cells were from adjacent nontumor samples. A total of 2,320 DEGs were identified in these cellular populations, comprising 752 significantly upregulated genes and 1,568 significantly downregulated genes (Fig. 2A). In the TCGA-LIHC dataset, 5,586 DEGs were detected, with 2,637 demonstrating increased expression and 2,949 showing decreased expression levels (Fig. 2B). In addition, the analysis revealed that 4,971 genes were significantly associated with HCC recurrence outcomes. By intersecting these results, 53 RFS-related DEGs were identified that were incorporated into the prognostic signature (Fig. 2C).

Through regularization with LASSO Cox regression, a risk score prognostic signature based on six genes [cyclin-dependent kinase inhibitor 2A (CDKN2A), complement factor H-related 3 (CFHR3), cytochrome P450 family 2 subfamily C member 9 (CYP2C9), high mobility group box 2 (HMGB2), IGLC2 and Jupiter microtubule-associated homolog 1 (JPT1)] was constructed (Fig. 2D and E). The risk score was calculated as follows: Risk score=(0.009 × expression of CDKN2A)-(0.019 × expression of CFHR3)-(0.010 × expression of CYP2C9) + (0.019 × expression of HMGB2)-(0.012 × expression of IGLC2) + (0.048 × expression of JPT1). Time-dependent ROC analysis demonstrated robust predictive accuracy for RFS, with 6-, 12- and 24-month area under the curve (AUC) values of 0.747, 0.742 and 0.721, respectively (Fig. 3A). Furthermore, patients in the high-risk subgroup exhibited significantly worse survival outcomes than those in the low-risk subgroup (Fig. 3B).

Differential expression analysis identified 525 genes distinguishing high- and low-risk patients. KEGG pathway enrichment revealed predominant involvement in the ‘cell cycle’, ‘metabolism of xenobiotics by cytochrome P450’, ‘retinol metabolism’, ‘drug metabolism-cytochrome P450’, ‘bile secretion’ and ‘chemical carcinogenesis-DNA adducts’ (Fig. 4A). Furthermore, GO analysis was used to characterize the functional attributes of the identified DEGs systematically (Fig. 4C). Biological processes were predominantly associated with ‘nuclear division’, ‘chromosome segregation’ and ‘carboxylic acid biosynthetic process’, cellular components were enriched in ‘spindle’, ‘chromosomal region’ and ‘chromosome, centromeric region’, and molecular functions were significantly represented by ‘ATPase activity’, ‘iron ion binding’ and ‘monooxygenase activity’. The results of the GSEA analysis of the six genes are presented in Fig. S1. In addition, CIBERSORT analysis revealed that the high-risk cohort had significantly fewer naive and memory B cells, memory resting CD4⁺ T cells and resting mast cells than the low-risk cohort, indicating immunosuppression (Fig. 4B). Moreover, the stromal score, immune score and estimate score of the high-risk group were significantly lower than those of the low-risk group (Fig. 4D).

Univariate Cox regression revealed that T stage, pathological stage and the risk score were significant predictors of HCC recurrence (Fig. 5A). However, subsequent adjustment demonstrated that only the prognostic signature retained independent predictive value, suggesting that HCC management could be augmented by incorporating the molecular signature assessment for more precise recurrence prediction (Fig. 5B).

As an externally validated dataset, the GSE76427 dataset was used to independently assess robustness and clinical applicability of the prognostic signature. Time-dependent ROC analysis demonstrated strong discriminative ability at 6 months (AUC, 0.840), with a maintained predictive value at 12 months (AUC, 0.596) and 24 months (AUC, 0.557) (Fig. 3C). Consistent with these findings, Kaplan-Meier analysis revealed significantly worse outcomes in the high-risk patients compared with that in the low-risk patients (Fig. 3D), indicating the robust predictive capacity of the integrated single-cell and bulk transcriptome-derived signature.

The RT-qPCR results revealed that CDKN2A, HMGB2 and JPT1 were significantly upregulated in HCC tumor tissues compared with in adjacent normal tissues (Fig. 6A-C), whereas CFHR3, CYP2C9 and IGLC2 were significantly downregulated in tumor tissues, compared with that in adjacent normal tissues (Fig. 6D-F). Notably, these experimental findings corroborated those of the prognostic signature, reinforcing its validity and implying potential functional roles for the signature genes.