Introduction

Hepatocellular carcinoma (HCC) is the predominant form of primary liver cancer, with its occurrence rate on the rise worldwide.1 Targeted therapy and immunotherapies have advanced rapidly, and the FDA have approved immune checkpoint inhibitors (ICIs) and vascular endothelial growth factor inhibitors as the primary treatment for advanced HCC. Yet overall survival rates remain low, with a five-year survival rate of below 20%.1–3 The poor prognosis of HCC is attributed to factors such as delayed diagnosis, high recurrence rates, and strong resistance to conventional chemotherapy and targeted therapies.4,5 Consequently, there is an immediate necessity to establish more accurate prognostic biomarkers and individualized therapeutic approaches.

In recent years, post-translational modifications (PTMs) have garnered significant attention for their regulatory roles in tumorigenesis and progression.6,7 Among these, neddylation, as an important ubiquitin-like modification, has gradually revealed its unique functions.8 Neddylation refers to the covalent linkage between the neural precursor cell expressed developmentally downregulated protein 8 (NEDD8) molecule and target proteins, which regulates protein stability, activity, subcellular localization, and their involvement in signal transduction processes.8–10 Previous studies have shown that abnormal activation of neddylation is present in various tumors, particularly in HCC, where it is significantly associated with cellular metabolic reprogramming, proliferation, and survival, and is closely linked to the activation of key oncogenic signaling pathways.1,3 Additionally, the emergence of NEDD8 activating enzyme inhibitor MLN4924 has provided a new targeted direction for HCC treatment. Studies have shown that this drug is associated with reduced tumor proliferation by inhibiting the function of the Cullin-RING E3 ligase (CRLs), which may be related to the induction of autophagy and apoptosis in HCC cells.1,8,11

Although the pro-tumorigenic role of the neddylation pathway in HCC has been reported, systematic studies on its associated genes in HCC prognosis remain limited. Current studies predominantly rely on bulk RNA sequencing data.12,13 Although they can identify differentially expressed genes (DEGs) on a global scale, they are unable to analyze cellular heterogeneity within tumor tissues, especially the signal regulation of non-tumor cells such as immune cells, tumor-associated fibroblasts, and vascular endothelial cells in the tumor microenvironment (TME).13 The swift advancement of single-cell RNA sequencing (scRNA-seq) technology has offered a novel perspective for identifying the expression patterns of neddylation-related genes (NRGs) in different cell types within the TME.10,14 Recent studies have demonstrated that dysregulation of neddylation is closely implicated in HCC initiation, progression, and therapeutic resistance, and targeting this pathway holds great potential for prognostic evaluation and individualized treatment.15–18 At the same time, the progression of HCC is highly dependent on the interactions between neoplastic cells and the immunological microenvironment. Factors such as immune cell infiltration composition, immune checkpoint expression, and matrix remodeling jointly influence disease progression and treatment response.2,3 However, most current prognostic models fail to integrate immune microenvironmental features with intrinsic tumor cell signaling pathways, making it difficult to reflect biological differences among patients and limiting their clinical application value.

To address these gaps, this study extensively examined the expression characteristics and survival associations of key NRGs in HCC, combined with data from the The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, and constructed and validated a robust prognostic model. The expression of prognostic genes in HCC samples was validated via real-time quantitative PCR (RT-qPCR). Subsequently, comprehensive immune-infiltration and drug-response analyses were performed to characterize immune landscape differences between risk subgroups and uncover potential therapeutic targets. At the single-cell level, the expression signatures of prognostic genes across distinct cellular subpopulations and the intercellular communication networks in HCC were further delineated. These efforts provide theoretical basis for precision treatment and risk stratification in HCC.

Materials and MethodsData Collection

The transcriptional data and clinical information for TCGA-liver hepatocellular carcinoma (LIHC) were retrieved from the TCGA database (https://portal.gdc.cancer.gov/),19 serving as the training set. This dataset comprised 421 samples, including 371 primary tumor tissues and 50 normal adjacent control tissues, with 365 cancer samples having complete survival information. Two additional datasets were delivered from the GEO database (https://www.ncbi.nlm.nih.gov/geo):19 GSE14520 and GSE149614. GSE14520, based on the GPL3921 sequencing platform, included 225 HCC liver tumor tissues and 220 normal tissue samples and was selected as the validation cohort for subsequent analysis. GSE149614, which was scRNA-seq data sequenced on the GPL24676 platform, included 10 HCC tumor tissues and 8 control samples. The GeneCards database (https://www.genecards.org/)20 was queried using the keyword “neddylation”, and 2,227 NRGs were downloaded with the screening criterion of score > 1.

Differential Expression Analysis

In the training cohort, the DESeq2 package (v 1.38.3)21 was executed to screen the DEGs between the HCC and normal groups, with |log2-fold change| > 1 and false discovery rate (FDR) < 0.05 considered as significantly different DEGs. Subsequently, the ggplot2 package (v 3.5.1)22 and the ComplexHeatmap package (v 2.14.0)23 were executed to draw volcano plot and heatmap, respectively, to intuitively display the differences in gene expression.

Identification of Intersection Genes and Enrichment Analysis

Based on the training cohort, the survival package (v 3.4–0)24 was utilized to conduct univariate Cox regression analysis on the DEGs in the HCC samples to identify survival-predictive and prognosis-associated genes, with hazard ratio (HR) ≠ 1 and P < 0.05, while also satisfying the proportional hazards (PH) assumption test. The ggvenn package (v 0.1.10)25 was further executed to visualize the convergence of DEGs associated with survival and NRGs, obtaining the intersection genes. Subsequently, the clusterProfiler package (v 4.6.2)26 was executed to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on the intersection genes, with the threshold set at adj.P < 0.05, and the enrichment results were depicted utilizing the ggplot2 package and the cnetplot function.

Construction of Risk Model

The glmnet package (v 4.1–8)27 was executed to implement least absolute selection and shrinkage operator (LASSO)-Cox regression analysis on the intersection genes, with 10-fold cross-validation conducted. On this basis, multivariate Cox regression analysis was conducted on the screening results (HR ≠ 1, P < 0.2), and stepwise regression was used to identify prognosis-predictive genes. After calculating the individual risk scores employing the predict function, the survminer package (v 0.5.0)28 and its surv_cutpoint function were employed to divide the HCC patients into groups of high and moderate risk (HRG and LRG) based on the risk scores (minprop = 0.25), and Kaplan-Meier (KM) survival curves were drawn between the groups. At the same time, risk curves were plotted, and heatmaps illustrating expression levels of prognostic genes across several risk categories were displayed. The timeROC program (v 0.4)29 was implemented to formulate receiver operating characteristic (ROC) curves to figure out the area under the curve (AUC) values for 1, 2, and 3 years as survival time points in order to evaluate the predictive performance of the prognostic model. Finally, the robustness of the risk model was further verified in the validation cohort.

Construction of Nomogram

The correlation between the risk score and 5 clinical features (age, stage, N, T, and gender) was assessed. Subsequently, the risk score and these 5 clinical features were included in univariate Cox regression analysis to verify whether they were independent predictive factors for overall survival (OS). Considering the independent prognostic variables that were found, the rms package (v 6.8–1)30 was executed to construct nomogram, and the results were visualized using the regplot package (v 1.1).31 To assess the predictive performance and accuracy of the model, calibration and ROC curves were prepared. The package (v 1.2)32 was executed to perform decision curve analysis (DCA) to assess the clinical benefit of the diagnostic model.

Immune Microenvironment Analysis

The estimate package (v 1.0.13)33 was executed to assess the differences in ImmuneScore, StromalScore, and EstimateScore among different risk categories. In the training cohort, the CIBERSORT package and the LM22 gene set34 were employed to calculate the proportions of 22 immune cell types. Bar graphs of the immune cell infiltration abundance percentages in various groups were created with the ggplot2 software. All immune cells and prognostic genes were examined utilizing Spearman correlation, and heatmap was generated. In addition, based on the 63 common immune checkpoints reported in the literature,35 differential expression analysis was performed to gauge the alterations in immune checkpoints between HRG and LRG.

Drug Sensitivity Prediction and Tumor Immune Dysfunction and Exclusion (TIDE)

The transcriptional data of the training cohort were standardized based on the mean expression values of normal samples and then applied to the TIDE website (http://tide.dfci.harvard.edu) to explore the differences in immune therapy sensitivity between HRG and LRG through TIDE scores. Based on the reference cell line expression matrix and drug treatment information provided by oncoPredict (v 1.2)36 for GDSC v2 database, drug sensitivity prediction was performed. Subsequently, Wilcox tests were conducted on the 50% inhibitory concentration (IC50) values of each drug between HRG and LRG, and the top 6 drugs with significant differences were displayed.

Human Protein Atlas (HPA)

The expression atlas of prognostic genes in normal and HCC tissues was deployed from the HPA (https://www.proteinatlas.org/).37 The image format was converted using ImageJ software (version 1.53e), and the mean optical density (MOD) values were calculated.

Mendelian Randomization (MR) Analysis

MR was conducted with prognostic genes as exposure variables and HCC as the outcome. The analysis strictly adhered to the MR reporting guidelines for observational studies (STROBE-MR) (Table S1).38 The genome-wide association study (GWAS) dataset for HCC (ieu-b-4953) was acquired from the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/),39 which included 168 HCC and 372,016 control European samples, involving 6,304,034 single nucleotide polymorphisms (SNPs). Additionally, cis-expression quantitative trait loci (cis-eQTL) data for prognostic genes (Table S2) were obtained from the eQTLGen Consortium (https://www.eqtlgen.org/) to further support the analysis. The instrumental variables (IVs) must meet three basic assumptions: (1) IVs must exhibit a robust association with the exposure factor being studied; (2) IVs must be independent of any known or unknown confounding variables; (3) IVs must solely affect the outcome via the exposure factor, without any direct causal paths. The SNP selection threshold was set at P < 5*10−6. The ieugwasr package (v 1.0.0)40 was implemented to exclude SNPs exhibiting linkage disequilibrium, with parameters configured at r2 = 0.1; kb = 100. Subsequently, the f-statistic of each genetic variant was calculated, and only those exhibiting an f-statistic over 10 were retained. Using the GWAS Catalog (https://www.ebi.ac.uk/gwas/), SNPs potentially related to the outcome GWAS traits were excluded at a threshold of P < 1*10−5. The harmonise_data function in the TwoSampleMR package (v 0.6.3)41 was employed, and IVs significantly correlated with the outcome variable were excluded. Subsequently, the mr function was used to assess causal relationships using various MR methods, with the inverse variance weighting (IVW)42 as the principal methodology, augmented by MR-Egger, weighted median, simple mode, and weighted mode. The forestploter package (v 1.1.2)43 was executed to draw forest plots for visualization. Sensitivity test was executed to assess the dependability of the MR analysis results. When the P-values of the heterogeneity test and horizontal pleiotropy test were greater than 0.05, it indicated the absence of horizontal pleiotropic effects and significant heterogeneity. In addition, leave-one-out analysis was conducted by recalculating the IVW causal estimates after removing one SNP at a time to verify whether the estimates were biased or driven by outliers.

ScRNA-Seq Analysis

The Seurat package (v 5.1.0)44 was executed to conduct strict quality control and preprocessing in the GSE149614 dataset. Genes covered by fewer than three cells and cells containing fewer than 200 genes were eliminated, as well as cells with more than 15% mitochondrial genes and cells with ≤ 200 or ≥ 7,000 genes, and counts ≤ 200 or ≥ 80,000. After standardization, 2,000 genes with high intercellular variation coefficients were extracted using the FindVariableFeatures function’s variance stabilizing transformation (vst) technique. Batch correction was performed using integrated.rpca. The ScaleData function was employed to normalize the data. The top 50 primary components were picked for further examination after the principal components with statistical significance emerged via the JackStrawPlot and JackStraw functions. The FindNeighbors and FindClusters functions of the Seurat package were executed to identify small cell clusters (resolution set at 0.2), and uniform manifold approximation and projection (UMAP) were employed for clustering analysis of cellular aggregates. The FindAllMarkers function was executed to identify characteristic DEGs in each cluster, and cell types were annotated according to the top 100–200 characteristic genes in each category using the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html). The spatial distribution and expression levels of prognostic genes across various cell types and tissues were then analyzed. To study the potential interactions between HCC and normal samples, the CellChat package (v 1.6.1)45 was implemented to analyze cell-cell communication.

Clinical Sample Collection

Between June and August 2025, eight pairs of hepatocellular carcinoma (HCC) tissues and corresponding adjacent non-tumor tissues (control) (at least 5 cm from the tumor margin) were procured from the First Affiliated Hospital of Dalian Medical University. HCC cases were pathologically confirmed, with exclusions for other types of liver cancer, other concurrent cancers, severe complications, resectable tumors, prior systemic therapy, negative composite biomarkers, Child–Pugh class C or specific class B liver dysfunction, substandard organ function, or high HBV-DNA without antiviral treatment. Control tissues were confirmed pathologically normal and excluded for tumor infiltration, inflammation, fibrosis, or other malignancies. All participants read and fully understood the pre-provided written informed consent forms. The Ethics Committee of the First Affiliated Hospital of Dalian Medical University approved the study (No. PJ-KS-KY-2025-689) and conducted in accordance with the Declaration of Helsinki.

RT-qPCR

Total RNA was extracted from both HCC tissues and matched adjacent non-tumor tissues utilizing the SteadyPure Rapid RNA Extraction Kit (AG21023, ACCURATE Biotechnology, Hunan, China), following the manufacturer’s instructions. Subsequently, cDNA was synthesized from the extracted total RNA utilizing the Evo M-MLV RT Kit with gDNA Clean for qPCR (AG11728, Accurate Biotechnology (Hunan) Co., Ltd., China). RT-qPCR was performed using the SYBR Green Premix Pro Taq HS qPCR Kit (AG11701, Accurate Biotechnology (Hunan) Co., Ltd., China) on a QuantStudio™ Real-Time PCR System (Applied Biosystems). The amplification process included an initial denaturation at 95°C for 30 seconds, succeeded by 40 cycles of 95°C for 5 seconds and 60°C for 30 seconds. The table of primer sequences was listed in Table 1. The expression levels of target genes were adjusted to the endogenous control β-actin and calculated using the 2−ΔΔCt technique. All reactions were conducted in triplicate to ensure technical reproducibility.

Table 1 Primer Sequences Used in This Study

Statistical Analysis

All analyses were performed in R version 4.4.1, with inter-group differences assessed using Wilcoxon test or Log rank test. Data from at least three independent experiments are expressed as the mean ± standard error of mean (SEM). All statistical analyses were conducted with GraphPad Prism (version 10.1; GraphPad Software, San Diego, CA). Normality and homogeneity of variance were examined prior to parametric tests. For comparisons between two groups, normally distributed data were analyzed employing independent samples t-tests, while non-normally distributed data were assessed with non-parametric tests (eg., Mann–Whitney U-test). Univariate Cox regression analyses were subjected to the PH assumption test. A P-value of less than 0.05 was defined as statistically meaningful.

ResultsBiological Function Analysis of Intersection Genes

The training cohort’s HCC and control samples’ gene expression levels were scrutinized using differential analysis, identifying a total of 1,921 DEGs, of which 1,284 were upregulated and 637 were downregulated (Figure 1A and B). Subsequently, these 1,921 DEGs were included in univariate Cox regression analysis, ultimately identifying 528 genes significantly pertaining to survival. Among these genes, 495 had consistent differential expression directions in HCC and normal tissues with the risk direction from the univariate Cox regression analysis, and thus were retained for further analysis (Table S3). To further explore the potential associations between HCC and NRGs, an intersection analysis was conducted between these 495 genes and the 2,227 NRGs. Through this analysis, 62 intersection genes were successfully identified (Figure 1C). To further comprehend how these intersection genes act biologically, enrichment analysis was deployed. The GO enrichment analysis results unveiled that these genes were augmented in 35 significant entries, with 27 significant entries in biological processes and 8 in cellular components. These enriched entries were closely related to DNA-related processes, such as “DNA unwinding involved in DNA replication”, “regulation of DNA-templated DNA replication”, “DNA duplex unwinding”, and “positive regulation of DNA-directed DNA polymerase activity”, “regulation of DNA primase activity”, which played important roles in cell proliferation and maintenance of genomic stability. Additionally, they were involved in the regulation of enzyme activities, such as “regulation of phosphatidylinositol 3-kinase activity”, “regulation of cyclin-dependent protein kinase activity”, and “regulation of cyclin-dependent protein serine/threonin kinase activity”, which were keys in cell cycle progression and signal transduction (Figure 1D). KEGG enrichment analysis additionally elucidated the probable functions roles of these genes in metabolic and signaling pathways. A total of 9 significant entries were enriched, including “tryptophan metabolism,” “p53 signaling pathway”, “valine, leucine and isoleucine degradation”, “fatty acid degradation”, and “pathways in cancer” (Figure 1E).

Figure 1 Identification and enrichment analysis of intersection genes. (A and B): volcano plot and heatmap of differentially expressed genes (DEGs) between hepatocellular carcinoma (HCC) and control groups in TCGA-LIHC. (C): intersection genes between neddylation-related genes (NRGs) and survival-associated DEGs. (D): Gene Ontology (GO) enrichment analysis of intersection genes. Biological process, BP; cellular component, CC; molecular function, MF. (E): Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of intersection genes.

Identification of 6 Prognostic Genes and Development of Risk Model

To further explore prognostic biomarkers among the 62 intersection genes, LASSO-Cox regression analysis was conducted. By setting Lambda.min = 0.064, 11 characteristic genes were identified (Figure 2A). Subsequently, multivariate Cox analysis further narrowed down to 6 prognostic genes: SOCS2, DIRAS2, LPL, KRT17, BFSP1, and POF1B (Figure 2B). To further validate the expression profiles of these prognostic genes, relevant data were obtained from the HPA. The results unveiled that the expression of SOCS2 in HCC tissues (MOD = 0.578) was lower than that in normal tissues (MOD = 0.339), consistent with the previous differential analysis results (Figure 2C). However, MR analysis failed to establish a causal link between these prognostic genes and HCC (Table S4).

Figure 2 Identification of prognostic genes. (A): least absolute selection and shrinkage operator (LASSO)-Cox regression analysis to screen for feature genes. (B): multivariate Cox regression analysis of 11 feature genes [hazard ratio (HR) ≠ 1, P < 0.2] with stepwise regression to identify 6 prognostic genes. (C): Human Protein Atlas (HPA) analysis of SOCS2. Left, normal tissue; right, HCC tissue.

Based on these 6 prognostic genes, a risk prediction model was developed in the training cohort. Risk score = h(t)× exp[SOCS2×(−0.382)+DIRAS2×0.189+LPL×0.172+KRT17×0.136+BFSP1×0.475+POF1B×0.114]. The risk curve intuitively demonstrated that the number of deceased samples in the HRG exceeded that in the LRG, signifying that the risk score could effectively distinguish the prognosis of patients (Figure 3A). Heatmap analysis further revealed the expression characteristics of different genes in the high-risk state: DIRAS2 and BFSP1 had higher expression levels in the HRG, while the expression trend of SOCS2 was lower in the HRG (Figure 3B). This differential gene expression pattern provides important clues for understanding the prognostic-relevant molecular characteristics of HCC. KM survival curve analysis was conducted to evaluate the predictive capability of the risk model (Figure 3C). According to the data, there was a substantial difference in survival among the HRG and LRG (P < 0.0001), with patients in the HRG having a much lower survival status. Additionally, ROC curve analysis indicated that the AUC values of the prognostic risk model at 1, 2, and 3 years were all greater than 0.700, demonstrating the model’s high predictive accuracy (Figure 3D). The robustness and diagnostic performance of the risk model was then validated in the validation cohort, with results consistent with those of the training cohort (Figure 3E–H).

Figure 3 Construction and validation of the risk model. (A): risk curve in the training set. (B): heatmap of prognostic gene expression in the training set. (C): Kaplan-Meier (KM) curve of the risk model in the training set. (D): evaluation of the risk model using receiver operating characteristic (ROC) curve in the training set. (E): risk curve in the validation set. (F): heatmap of prognostic gene expression in the validation set. (G): KM curve of the risk model in the validation set. (H): evaluation of the risk model using ROC curve in the validation set.

Figure 3 continued.

Further Improvement of Prognostic Model by Combining Risk Score with Clinical Pathological Features

To construct a more precise prognostic model, the relationship between 5 clinical features and the risk score was thoroughly analyzed. The results unveiled that clinical stage and T were correlated with the risk score (Figure S1A), suggesting that these clinical features are significantly associated with patient prognosis. Subsequently, these factors were incorporated into univariate Cox analysis, with stage, T stage, and risk score being significantly linked to survival (Figure S1B). Based on this, stage, T stage, and risk score were incorporated as independent prognostic variables into the nomogram to forecast the OS probability of patients (Figure 4A). Analysis of the calibration curve revealed the nomogram’s anticipated survival rates and the actual observed values were highly consistent, indicating that the model could precisely forecast patient survival (Figure 4B). Further ROC curve analysis demonstrated that the AUC values were all greater than 0.7, representing an improvement compared to the risk score alone, representing an improvement compared to the risk score alone and clinical stage alone (Figure 4C). A comparison of AUC values between our prognostic model and published HCC prognostic models indicated that our risk score achieved higher AUC values than clinical stage at all time points and performed comparably to previously reported signatures with AUC values in the range of 0.70–0.75, which is a common and clinically meaningful level for HCC prognostic models.46–48 The concordance index (C-index) of the clinical model (T stage + stage) was 0.659, which was elevated to 0.707 after integrating the risk score. Moreover, significant positive net reclassification improvement (NRI) and integrated discrimination improvement (IDI) values were observed at 1, 3, and 5 year survival endpoints (all P < 0.05), confirming that our nomogram provides significant incremental prognostic value beyond conventional clinical predictors (Tables 2 and 3). Additionally, DCA indicated that the model had good net benefits, further supporting its value in clinical prognostic assessment (Figure 4D).

Table 2 Concordance Index Comparisons

Table 3 Net Weight Classification Improvement Index (NRI) and Comprehensive Discrimination Improvement Index (IDI)

Figure 4 Construction of the nomogram. (A): nomogram constructed with stage, T, and risk score as variables. (B): calibration curve of the nomogram. (C): ROC curve of the nomogram. (D): decision curve analysis (DCA) curve of the nomogram.

Immune Features Between Different Groups

To investigate the disparities in the immunological microenvironment across various categories, immune scoring analysis was first conducted. The results unveiled a significant decrease in StromalScore in the HRG, suggesting that the infiltration of stromal cells might be lower in this group, indicating potential associations with TME remodeling (Figure 5A). Substantial discrepancies in the infiltration levels of eight immune cell types across risk categories were detected by further analysis of immune infiltration (P < 0.05) (Figure 5B and C). Specifically, the infiltration levels of “M0 macrophages” and “follicular helper T cells” were higher in the HRG, while “resting memory CD4+ T cells” and “resting mast cells” were higher in the LRG. These findings indicated variations in immune cell infiltration patterns among different categories. Additionally, correlation analysis showed that M0 macrophage had the most pronounced negative connection with SOCS2 (cor = −0.303, P < 0.0001) and the strongest positive correlation with DIRAS2 (cor = 0.291, P < 0.0001), revealing the complex predictive and correlative relationships between immune cells and gene expression (Figure 5D). Furthermore, the above immune infiltration patterns were successfully validated in the independent GSE14520 cohort, where increased infiltration of “M0 macrophages” and decreased infiltration of “resting mast cells” in the high-risk group were consistently replicated (Figure S2A–C). Notably, “M0 macrophages” showed the most significant negative correlation with SOCS2 in both the training and validation cohorts (Figure S2D).

Figure 5 Immune microenvironment analysis. (A): differences in ImmuneScore, StromalScore, and EstimateScore among different risk groups. (B): stacked bar chart of immune cell scores in high- and low-risk groups. (C): boxplot showing differences in immune cells between high- and low-risk groups. (D): correlation plot between prognostic genes and immune cells.

Note: *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Figure 5 continued.

Differences in Immune Therapy

To comprehensively assess the expression patterns of immunological checkpoints in HCC patients across various risk groups, differential analysis was conducted on 63 common immune checkpoint molecules. The results showed significant differences in 13 immune-activating molecules, 14 immune-inhibitory molecules, and 16 bi-directional immune molecules between the two categories (P < 0.05) (Figure 6A). Specifically, immune-activating molecules such as TNFRSF4, CD276, and HLA-DRA were significantly upregulated in the HRG, while immune-inhibitory molecules such as TDO2 and BTNL9 were significantly downregulated. These differentially expressed immune checkpoint molecules are associated with the disparities in immunological microenvironments across several risk categories, indicating that the model may predict differential efficacy of immune therapy. To further investigate the differences in immune therapy efficacy between HRG and LRG, the TIDE scores of tumor samples were assessed (Figure 6B). The results showed significant differences in TIDE scores between different risk groups, indicating that individuals with different disease risks might have different responses to immune therapy. Additionally, based on GDSC drug sensitivity analysis, significant differences in the IC50 values of 122 drugs were found among different risk groups (P < 0.05). To more intuitively display these results, boxplots were used to visualize the differences in the top 6 drugs (Figure 6C). Specifically, in the HRG, the IC50 values of Vinblastine_1004, Paclitaxel_1080, MG−132_1862, and Sepantronium bromide_1941 were significantly lower, indicating higher sensitivity to these drugs in HRG patients. Conversely, the IC50 values of Doramapimod_1042 and JAK1_8709_1718 were significantly higher in the HRG. These findings provide preliminary predictive references for drug selection and immunotherapeutic evaluation based on risk stratification, which deserve further prospective clinical validation.

Figure 6 Immune checkpoints and drug sensitivity. (A): differences in 63 immune checkpoint molecules among different risk groups. (B): tumor immune dysfunction and exclusion (TIDE) score between high- and low-risk groups. (C): boxplot of the top 6 commonly used chemotherapeutic drugs’ 50% inhibitory concentration (IC50) values between high- and low-risk groups.

Note: *, P < 0.05; **, P < 0.01; ***, P < 0. 001; ****, P < 0.0001.

Characteristics of Prognostic Genes from a Single-Cell Perspective

After identifying the prognostic gene signature, we further explored its expression distribution and predictive characteristics at the single-cell level. Strict filtering of the data was performed, successfully capturing the maximum variability of the top 2,000 highly variable genes between cells (Figure S3A and B). Based on this, the top 50 principal components were picked for examination analysis (Figure S3C). Through clustering analysis, 16 unique cell clusters with distinct characteristics were identified (Figure S3D). By comparing these cell clusters with known cell type marker genes, 16 different cell types were successfully annotated, including cholangiocytes, macrophages, endothelial cells, mesenchymal cells, B cells, mast cells, and others (Figure 7A and B). Bar plot was drawn to analyze the ratio of each cell type in the samples (Figure 7C). Additionally, it was determined that the overall expression levels of all genes were relatively low across different cell types (Figure 7D). However, it is worth noting that the SOCS2 gene had higher expression levels in endothelial cells and significant differential expression in most cell types (Figure 7E). Furthermore, cell-cell communication is an important mechanism for maintaining tissue homeostasis and regulating cell functions. The analysis of the communication network revealed the quantity of cell-cell interactions in HCC samples was significantly higher compared to normal samples, while the interaction strength was lower (Figure 8A and B). In HCC, the input and output signals of macrophages and exhausted T cells were significantly weakened (Figure 8C), which may be associated with altered functional states in the TME and dysregulated cell signaling.

Figure 7 Single-cell RNA sequencing (scRNA-seq) analysis. (A and B): annotated to 16 cell types by matching with marker genes. (C): proportion distribution of various cells in different samples. (D): expression of prognostic genes in different cells. (E): differences in prognostic genes in various cell types between HCC and normal samples.

Note: *, P < 0.05; **, P < 0.01; ****, P < 0.0001.

Figure 8 Cell communication analysis. (A and B): interaction and strength between HCC and normal groups. (C): changes in interaction strength of cell communication among different cells.

Figure 8 continued.

The Expression of Prognostic Genes Was Verified by RT-qPCR

To further validate the expression of prognostic genes (SOCS2, DIRAS2, LPL, KRT17, BFSP1, and POF1B) in a real-world clinical setting, we performed RT-qPCR analysis on a cohort comprising surgically resected tumor tissues and corresponding neighboring non-tumor tissues from individuals with HCC. The significant downregulation of SOCS2 and upregulation of DIRAS2, LPL, KRT17, BFSP1, and POF1B in tumor tissues validated the results of the bioinformatics differential analysis (Figure S4), providing further evidence for the predictive and prognostic value of these genes.

Discussion

Neddylation, a recently characterized post-translational modification, is significantly associated with immune regulation, metabolic reprogramming, and the pathogenesis of HCC.49 Emerging evidence suggests that alterations in neddylation are closely linked to HCC progression.50 Based on this, by integrating multiple analytical strategies, a stable and highly predictive prognostic model associated with neddylation was constructed, providing a tool for HCC risk stratification and personalized treatment. The model includes six prognostic genes: SOCS2, DIRAS2, LPL, KRT17, BFSP1, and POF1B. The gene expression of such prognostic in HCC was validated by RT-qPCR. Further analysis showed that the risk score model demonstrated good discriminatory ability and predictive performance within both the training set and outside verification set. Moreover, integrating the risk score with clinical parameters to construct a nomogram further improved the operational feasibility of model and clinical applicability.

In the HCC prognostic model constructed in this study, six prognostic genes demonstrated multidimensional biological significance in the emergence and progression of hepatic carcinoma. Prior research have shown that SOCS2, as an inhibitor of the JAK/STAT pathway, exerts antitumor effects in various solid tumors, and its elevated expression is frequently correlated with better prognosis.51–53 The study found that the expression of SOCS2 in HCC tissues was reduced,54 which was consistent with our RT-qPCR results. Mechanistically, SOCS2 overexpression has been reported to suppress proliferation, migration, and invasion of HCC cells.55 In line with previous evidence, our study also found that SOCS2 is significantly associated with favorable prognosis in HCC, suggesting that it may inhibit the release of pro-inflammatory factors in the TME by negatively regulating inflammatory signaling pathways, thereby slowing tumor progression.56 In contrast, KRT17, BFSP1, and LPL have been reported to be intimately linked with tumor invasiveness, growth rate, and immune evasion capacity in multiple cancer types.57–59 KRT17 is not only a component of the cytoskeleton but also acts as a signaling regulator involved in proliferation and immune regulation. HCC tumor tissue samples exhibited high levels of KRT17 expression,60 which was consistent with our results. The overexpression of KRT17 is commonly observed in epithelial tumors, where it activates the AKT/mTOR pathway, enhancing tumor cell survival and migration capacity.57,59 BFSP1 and LPL are associated with lipid metabolism and lipid droplet accumulation, respectively, and such metabolic reprogramming is an important feature of HCC progression.5,7,59 BFSP1 and LPL are upregulated in HCC samples,58,61 and our results also show their upregulated expression. Additionally, although DIRAS2 and POF1B have been reported less frequently in liver cancer research, the literature has indicated that members of the DIRAS family may participate in cell cycle regulation and Ras signaling pathway inhibition, while POF1B plays a role in cell junctions and polarity maintenance, potentially influencing tumor metastasis potential by regulating cell adhesion DIRAS2.62,63

The TME plays an essential role in regulating the response to anti-tumor therapies in HCC. The examination of immune infiltration demonstrated substantial disparities in M0/M1 macrophages between risk groups, and scRNA-seq also annotated macrophages. Macrophages constitute a substantial ratio of immunological cells and exert key functions in the TME of HCC.64 M0 macrophages are characterized as immature macrophages capable of polarizing