The RNA sequencing data and corresponding clinical information for LIHC were obtained from the TCGA database (https://portal.gdc.cancer.gov/) (356 samples). Data for 212 samples from the International Cancer Genome Consortium Liver Cancer (ICGC-LIHC) cohort were sourced from the LIHCDB V2.0 database (http://lifeome.net:809/#/home), and data for another 206 samples from an ICGC-LIHC cohort were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Specifically, the dataset GSE124751, comprising 206 samples, was retrieved from GEO. The “limma” package was used to correct for batch effects between data sets. After excluding patients with incomplete survival data, a total of 781 tumor samples were retained for subsequent analysis (9).

To identify genes significantly associated with histone acetylation, we utilized the GeneCards database (https://www.genecards.org/), a comprehensive and integrative database of human genes. We performed a search using the keyword “Histone Acetylation.” GeneCards calculates a Relevance Score for each gene, representing the strength of the association between the gene and the search term based on weighted evidence from diverse sources, including experimental data, genetic pathways, and literature text mining. To ensure the reliability and specificity of the retrieved gene set, we applied a filtering threshold to exclude genes with weak associations. Only genes with a Relevance Score greater than 30 were retained. This rigorous screening process resulted in a final set of 588 histone acetylation-related genes for subsequent analysis.

The full single-cell RNA sequencing (scRNA-seq) datasets GSE149614 (6) and GSE166635 (7) were acquired in the GEO database. The Seurat package (8) was used for quality control and downstream analysis of scRNA-seq data. We excluded cells expressing fewer than 200 genes and genes detected in fewer than 5 cells and retained cells with a gene count between 300 and 5000. To ensure adequate sequencing depth while removing technical outliers, the total RNA count per cell (nCount_RNA) was kept above 1000 but below the 97th percentile of the dataset’s count distribution. The PercentageFeatureSet function was used to calculate the percentage of mitochondrial gene expression in each cell, excluding cells with >20% mitochondrial gene expression or >5% hemoglobin gene expression. The final curated dataset contains 86,409 cells in 23 samples.

The “RunHarmony” function was used to correct batch effects in the scRNA-seq data. The optimal number of principal components (PCs) for downstream analysis was determined by the “ElbowPlot” function. Cell clustering was performed using the “FindClusters” function, followed by dimensionality reduction visualization using “RunUMAP” (dims=1:10). Cell populations were manually annotated based on classical cell type markers reported in the literature.

Frozen tissues and corresponding non-tumor tissues from five patients diagnosed with liver malignancies were collected from the Second Affiliated Hospital of Xi’an Jiaotong University. All patients had a pathological diagnosis of hepatocellular carcinoma (LIHC) after surgical resection. Samples were used for WB and qRT-PCR experiments. All participants provided informed consent before the commencement of the study, which received approval from the Ethics Committee of the Second Affiliated Hospital of Xi’an Jiaotong University.

MHCC-97H, Hep3B, and LO2 cells were purchased from the Cell Bank of the Chinese Academy of Sciences and cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 μg/ml streptomycin in a cell culture incubator at 37°C, 5% CO2.

The Kaplan-Meier method was used to plot the overall survival (OS) and disease-free survival (DFS) curves, and the log-rank test was performed to compare the survival differences between the high-risk and low-risk groups (P < 0.05 was considered statistically significant). In the ICGC-LIRI-JP cohort, subgroup analyses were stratified by age (≤60 years/>60 years), gender, viral etiology (HBV/HCV), TNM stage (III/III-IV), and fibrosis stage (F0-F2/F3-F4). In the GSE124751 cohort, subgroups were stratified by age (≤60 years/>60 years) and tumor mutational burden (TMB, divided into high/low TMB groups using the median as the cutoff). The prognostic value of the risk model was validated in each subgroup. All survival analyses were performed using the “survival” package in R software (version 4.2.1), and Kaplan-Meier curves were plotted using the “survminer” package.

The “FindMarkers” function of the Seurat package (8) was applied to identify differentially expressed genes (DEGs) in cancer and normal tissues. Statistical significance was determined using the Wilcoxon test (threshold |logFC|>0.5 and corrected p-value <0.05, default for remaining parameters). The DEGs of cancer and normal tissue cells were intersected with histone acetylation-related genes to obtain 66 genes and further intersected with all genes of TCGA-LIHC, ICGC-LIHC, and GSE124751 tumor samples to obtain 65 genes, defined as histone acetylation-related differentially expressed genes (HAc-related DEGs).

To build a robust prognostic signature, we conducted a systematic and unbiased exploration to identify the optimal modeling strategy. The entire cohort (n=781) was first randomly partitioned into a training set (70% of patients) and an independent test set (30% of patients). Using the Mime R package, we then constructed 101 distinct prognostic models based on the 11 core prognostic genes. These models were generated by creating unique combinations of 10 machine learning algorithms: Lasso, Ridge, Stepwise Cox (stepCox), CoxBoost, Random Survival Forest (RSF), Elastic Net (Enet), Partial Least Squares Regression for Cox (plsRcox), Supervised Principal Components (SuperPC), Generalized Boosted Models (GBM), and Survival Support Vector Machine (survival-SVM).

Each of the 101 combinations represented a distinct pipeline for feature selection and/or model training. All models were developed and tuned exclusively on the training set. For algorithms requiring hyperparameter tuning (e.g., RSF), a 10-fold cross validation was performed only within the training set to identify the optimal parameters. The final selection of the best-performing model was based purely on its predictive accuracy on the unseen independent test set, as measured by the Concordance Index (C-index). The C-index is a key metric for evaluating survival models, quantifying the concordance between the predicted risk order of patients and their actual survival time order. It is the accepted standard for clinical prognostic studies as it effectively handles the censored data common in survival analysis. This data-driven approach ensures that the chosen model has the best-demonstrated generalizability, free from researcher bias.

Patients in the training and test sets were categorized into high-risk and low-risk groups based on their risk scores. The “subplot” function of the “Mime1” package (17) was used to analyze the differences in overall survival (OS) between risk groups (log-rank test, p < 0.05). The diagnostic potential of the histone acetylation-related DEGs was assessed using the “timeROC” (18) and “pROC” packages (19) for analyzing receiver operating characteristic (ROC) curves. Risk score distributions were analyzed and visualized using boxplots.

To further explore the prognostic value of the model, the prognostic ability of the model was assessed using OS in different subgroups. All analyses were done in the “survival” package, and the results were visualized using the “ggplot2” package.

Differentially expressed genes (DEGs) with corrected p-value < 0.05 and |logFC| > 0.5 were screened based on RNA-seq data. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the “clusterProfiler” package (20). Based on the logFC-ranked DEGs, the gene set of the MSigDB database (v7.5.1) was used for Gene Set Enrichment Analysis (GSEA). The pathways ranked in the top 10 by normalized enrichment score (NES) were visualized, and correlations between the risk model and core genes within these pathways were analyzed. All visualizations were generated using the “enrichplot” and “ggplot2” packages.

The “ESTIMATE” package (21) was used to calculate the immunity score, stroma score, and tumor purity for all tumor samples (n=781) and to compare the differences between groups. The relative proportion of immune cells was estimated using the “CIBERSORT” package (22). Immune cell activity was assessed by the “xCell” package (23), and immune cell expression levels were assessed. The correlation of prognostic genes with the immune score, stroma score, and tumor purity was visualized using the “corrplot” package. In addition, risk models were analyzed for correlation with immune checkpoint-related genes, chemokines, chemokine receptors, and immune checkpoint blockade therapy (ICB).

Drug sensitivity prediction was performed using the “OncoPredict” package (24). The model was trained based on the GDSC2 dataset, and IC50 values were calculated for the samples. The Wilcoxon test was used to compare the drug sensitivity differences between groups.

Based on the somatic mutation data from the TCGA dataset, we used the “maftools” package (25) for mutation analysis. We also analyzed the tumor mutation load (TMB) between the two groups. The frequency of mutations and their distribution in specific genes were analyzed by “maftools” package.

Molecular data from three hepatocellular carcinoma (HCC) cohorts (TCGA-LIHC, MSK 2024, and CLCA 2024) were obtained from the cBioPortal database (https://www.cbioportal.org/). Mutation types (including point mutations, insertions/deletions) and copy number alterations (CNAs, such as amplifications, deep deletions, shallow deletions) of key genes in the risk model (e.g., HLA-B, NEU1, CDK4) were extracted. The “maftools” package in R software was used to generate an Oncoprint plot visualizing the gene alteration profiles of individual patients across cohorts. A stacked bar plot was created using the “ggplot2” package to compare the overall frequencies of major alteration types (amplification, deep deletion, mutation) among the three cohorts. For the NEU1 gene in the MSK-IMPACT cohort, the Kruskal-Wallis test was performed to compare the distribution of mutation counts across different CNA types (amplification, gain, diploid, shallow deletion; p < 0.05 was considered statistically significant).

Ferroptosis refers to impaired intracellular lipid peroxide metabolism and toxic lipid production, which induces cell death. m6A is a methylation of RNA, i.e., on the 6th nitrogen atom of adenine (A) in RNA, which affects mRNA stability, translational efficiency, alternative splicing, and localization. We analyzed the prognostic risk model to correlate with ferroptosis and m6A-related genes in LIHC. Ferroptosis-related genes were derived from Ze-Xian Liu et al. Systematic analysis of the abnormalities and functions of ferroptosis in cancer (26). m6A-related genes were derived from Juan Xu et al. on the molecular characterization and clinical significance of m6A regulators across 33 cancer types (27). Statistical significance (p < 0.05) was estimated using an independent t-test for normally distributed variables.

Intercellular communication networks are systematically analyzed and resolved at the cellular level using the “CellChat” (28) package. The method is based on a database of ligand-receptor interactions that quantifies and predicts major cellular communication pathways. A probabilistic statistical framework is used to infer cell type-specific ligand-receptor pair communication probabilities and quantify the strength of communication and the number of receptor-ligand pairs between different cell types. Identify key cellular communication patterns in the tumor microenvironment, with special attention to signaling pathways associated with NEU1 and differences between groups with high and low NEU1 expression.

The “CytoTRACE” (29) algorithm was used to assess the developmental potential of single cells, and the CytoTRACE score was calculated based on the number of genes detected in each cell, with higher scores indicating cells with higher differentiation potential. The distribution of CytoTRACE scores was compared between tumor epithelial cells and normal epithelial cells to assess the effect of cancer on the differentiation status of cells. In cancer epithelial cell populations, the CytoTRACE scores of high and low NEU1 expression groups were compared to explore the correlation between NEU1 expression levels and cell differentiation potential. The endothelial cell differentiation trajectories were constructed by applying the “Monocle3” package (30), and the cell distribution was visualized by the UMAP downscaling method and the master map was calculated to infer the cell differentiation trajectories. Along the inferred differentiation trajectory, the dynamic changes in the expression of NEU1 and its related genes were analyzed. Combining the proposed temporal information with differential expression analysis, the potential regulatory role of NEU1 in endothelial cell differentiation was revealed.

The NEU1 gene was mapped by the Coremine Medical Ontology Information Retrieval Platform (www.coremine.com/medical/) to screen its top 10 potential regulatory molecules. Download NEU1 structure files from the PDB database (https://www.rcsb.org/) to obtain target protein result files. Obtain the structure files of active compounds from PubChem database (https://pubchem.ncbi.nlm.nih.gov/), CB-Dock2 platform (31–33) to optimize the structure (removing water molecules, ligands), and molecular docking simulation to evaluate the binding role of receptor proteins to small-molecule ligands after hydrogenation and balancing the charge on receptor proteins.

Total RNA was extracted from tissues using TRIZOL reagent (Invitrogen), and the process was performed according to the instructions. RNA was reversed to cDNA using the PrimeScript RT kit (Takara). qRT-PCR experiments were completed using the Takara SYBR Premix Ex Taq II kit, and the process was performed according to the instructions. The final results were corrected for the expression of GAPDH.

The primer sequences used in this study were as follows:

Frozen tissues were lysed using RIPA buffer (Beyotime, China) containing a protease inhibitor mixture. Protein concentration was detected using the BCA Protein Quantification Kit (Beyotime, China). Protein samples were separated on 8% SDS-PAGE gels and transferred to PVDF membranes blocked with 5% non-fat powdered milk at low temperature for 2 h. The PVDF membranes were incubated with diluted primary antibodies and incubated overnight at 4°C (NEU1: 1:2000, proteintech, USA; GAPDH: 1:1000, CST, USA).

On the following day, after three washes with TBST, the PVDF membranes were incubated with HRP-conjugated secondary antibody (1:2000) for 2 hours. After three washes with TBST, protein bands were visualized using an ECL detection kit (Thermo Fisher, USA; 32106) and imaged with a Bio-Rad ChemiDoc Imaging System. Band intensity was quantified using ImageJ software.

NEU1-specific siRNA (Santa Cruz Biotechnology, sc-106297) was transfected into Hep3B and MHCC-97H cell lines using Lipofectamine™ 3000 (Thermo Fisher Scientific, L3000008) to knock down NEU1 expression. Knockdown efficiency was assessed by RT-qPCR and Western blotting.

For RT-qPCR analysis, cells were seeded in 12-well plates and transfected with 40 pmol siRNA per well. Total RNA was isolated 24 h post-transfection using TRIzol reagent (Invitrogen) and reverse-transcribed into cDNA with a PrimeScript RT kit (TaKaRa, RR036A). qPCR amplification was performed in 10 μL reactions using SYBR Green PCR Master Mix (TaKaRa, RR420B). Melting curve analysis confirmed primer specificity. The relative mRNA expression levels were normalized to GAPDH and calculated using the 2−ΔΔCt method. All primer sequences were synthesized by Sangon Biotech (Shanghai, China).

For Western blotting, cells were harvested 48 h after transfection. Total protein was extracted using RIPA lysis buffer (Beyotime, China; P0013B) supplemented with a protease inhibitor cocktail (Roche, Switzerland; 4693116001; 1 mg/mL). Protein concentrations were quantified using a BCA assay kit (Beyotime; P0010), and 20-30 μg of protein per sample was separated on 8% SDS-PAGE gels. Proteins were transferred to PVDF membranes (Millipore, USA), followed by blocking with 5% non-fat powdered milk for 2 h at room temperature. Membranes were incubated overnight at 4°C with primary antibodies against NEU1 (Proteintech, 67032-1-Ig; 1:1000) and GAPDH (Proteintech; 60004-1-Ig; 1:3000), followed by incubation with HRP-conjugated goat anti-rabbit secondary antibodies for 2 h at room temperature. Protein bands were visualized using an ECL detection kit (Thermo Fisher, USA; 32106) and imaged with a Bio-Rad ChemiDoc Imaging System. Band intensity was quantified using ImageJ software.

Transfected Hep3B and MHCC-97H cells were plated in 96-well culture plates at a density of 2×10³ cells per well. Following incubation under standard culture conditions (37°C, 5% CO₂) for 12, 24, and 36 hours, the culture medium was replaced with 100 μL of DMEM containing CCK-8 reagent (Beyotime Biotechnology, China; C0038) at a 1:10 (v/v) dilution. After 1 hour of additional incubation, optical density measurements were performed using a microplate reader at a wavelength of 450 nm.

Transfected Hep3B and MHCC-97H cells (1000–2000 cells/well) were plated in 6-well plates and cultured in a complete medium for 14 days. Colonies were fixed with 4% paraformaldehyde for 20 min, stained with 1% (w/v) crystal violet (Beyotime, China; C0121) for 20 min, and washed three times with deionized water. Colonies containing more than 50 cells were counted using an inverted microscope.

The ibidi Culture-Insert 2 Well (ibidi GmbH, Germany, 80209) was placed in a 60-mm Petri dish. Transfected Hep3B and MHCC-97H cells were seeded into the left and right chambers of the insert at a density of 2.5×10⁵ cells/well, respectively. After 12 hours of incubation in a serum-free medium, the culture insert was carefully removed. Cell migration patterns were documented at 0 h and 24 h using an inverted microscope.

Transwell inserts (Corning, USA; 3422) were coated with Matrigel matrix (Beyotime, China; Cat# C0376-5ml) diluted in DMEM at a 1:8 ratio (v/v). After 12 hours of serum starvation, cells were trypsinized and resuspended in a serum-free medium. Subsequently, 5×10⁴ cells per well were seeded in the upper chamber, while the lower chamber was filled with 600 μL of complete medium. Following 36 hours of incubation, cells were fixed with 4% paraformaldehyde for 20 minutes and stained with 1% crystal violet (Beyotime, China; C0121) for 20 minutes. After being washed three times with deionized water, non-invading cells on the upper membrane surface were gently removed using cotton swabs. Migrated cells were imaged with an inverted microscope, and five random fields per well were analyzed.

All statistical analyses and visualizations were done by RStudio (version 4.3.3) and GraphPad Prism 8.0. Comparisons of normally distributed continuous variables between two groups were performed using the independent samples t-test, and non-normally distributed variables were tested using the Wilcoxon rank sum test. The Welch one-way ANOVA was used to evaluate comparisons between several groups. Each experiment was repeated three times, and all in vivo experiments were performed with at least five animals per group. Data are expressed as mean ± standard error of the mean (SEM), with p < 0.05 considered statistically significant.

Figure 1 summarizes the overall workflow of this study. We first performed quality control and downstream analysis on two datasets, GSE166635 and GSE149614, from the Gene Expression Omnibus (GEO) database, and filtered the data to obtain a total of 86,409 cells, which were categorized into three groups of G1, G2M, and S phases, respectively (Figure 2A). Subsequently, the “RunHarmony” function was applied for batch removal, and principal component analysis (PCA) was used for dimensionality reduction, which successfully removed the batch effect while preserving the differences in the data (Figures 2B, C). We further analyzed the highly variable genes in the samples, and the results showed that 15 genes, including IGHG1 and IGKC, were significantly enriched in the high-expression-high-variation region (Figure 2D). Using the “ElbowPlot” function, we found that the standard deviation curve flattened out after the harmony dimension was about 10, indicating that the contribution of additional dimensions to the variation of the data gradually decreased. Therefore, we chose the first 10 harmony dimensions for downstream clustering and dimensionality reduction visualization analysis, balancing the trade-off between dimensionality reduction and computational complexity (Figure 2E). Following unsupervised clustering of the scRNA-seq data, we employed UMAP dimensionality reduction for visualization. This process partitioned the cells into seven distinct clusters, which were subsequently annotated based on canonical marker genes as: B cells, Endothelial cells, Fibroblast/Mesenchymal cells, Hepatocytes, Macrophages/Monocytes, Plasma cells, and T/NK cells. To demonstrate the effectiveness of our cross-platform normalization for the bulk transcriptomic data, we performed UMAP on the integrated patient samples. The resulting plot shows a homogenous mixture of samples from all cohorts, confirming that batch effects were successfully removed. Finally, analysis of cellular composition revealed that Hepatocytes, Macrophages/Monocytes, and T/NK cells were the most abundant populations in LIHC samples, whereas T/NK cells and Macrophages/Monocytes were the dominant populations in normal liver tissue (Figures 2F-H, Supplementary Figure 1).

In addition, to leverage the high resolution of single-cell RNA sequencing, we performed a cell-type-specific differential expression analysis. Specifically, within each annotated cell type (e.g., T cells, hepatocytes), we compared cells originating from LIHC tumor tissue with those from adjacent normal tissue. The lists of differentially expressed genes from each cell-type comparison were then merged and deduplicated, yielding a final set of 1,625 unique DEGs for subsequent analysis (Figure 2I, Supplementary Table 1).

To identify a set of candidate prognostic genes for subsequent analysis, we implemented a multi-step filtering process. We began with a list of 588 genes highly associated with histone acetylation, obtained from the GeneCards database (Figure 3A, Supplementary Table 2). Next, to ensure these genes were not only biologically relevant and actively dysregulated in the liver cancer microenvironment, but also measurably present in our study cohorts, we performed a crucial intersection. We intersected these 588 genes with both (1) a list of 1,625 differentially expressed genes (DEGs) identified from single-cell analysis, and (2) the set of all genes available in our combined transcriptomic dataset (TCGA, ICGC, and GSE). This rigorous process yielded a final set of 65 genes, which were defined as histone acetylation-related differentially expressed genes (Hac-DEGs) and were carried forward for all subsequent prognostic model development (Figure 3B, Supplementary Table 3).

Next, we went through 101 machine-learning combinations and obtained the C-index of each combination in all training and validation sets. The results showed that the StepCox[forward]+RSF combination had the highest C-index, which was 0.95 and 0.62 in the training and validation sets, respectively, suggesting that it had the best performance in predicting the prognosis of LIHC (Figures 3C, D). As a result, a histone acetylation-associated risk model containing 11 genes (HLA-B, HEXB, CDK4, ACAT1, NAA10, B2M, HSPD1, NPM1, PON1, NEU1, and CFB) was constructed (Figure 3E). We evaluated the expression of these 11 genes in the TCGA-LIHC dataset (Supplementary Figure 2) and their respective roles in assessing the overall survival (OS), disease-specific survival (DSS), progression-free interval (Progression Free Interval (PFI), and Disease-free-survival (DFS) (Supplementary Figure 3).

A risk score was assigned to each patient using the model, and all patients were included in either the high-risk or low-risk group based on the score. The OS of patients in different subgroups was analyzed in the training set using KM curves, and the results showed that patients in the high-risk group had a significantly worse prognosis (p<0.001, 95% CI:17.36-34.27); this result was also confirmed in the validation set (p=0.014, 95% CI:1.27-3.71), which suggests that the model demonstrates robust performance in evaluating the overall survival of LIHC patients. (Figures 3F, G). In addition, time-dependent ROC analysis showed that the model also had high value in predicting the overall survival of LIHC patients at 1, 3, and 5 years in both the training and validation sets (AUC: 0.948, 0.953, and 0.930 for the training set; AUC: 0.644, 0.628, and 0.605 for the validation set) (Figures 3H, I). We depicted the risk factor plots of the model in the training and validation sets, respectively, to demonstrate the relationship between the model score and the survival status of the LIHC patients, where green dots indicate the survival status and red dots indicate the death status. In both the training and validation sets, the number of red points increased as the risk score increased, indicating that patients with high scores on the model had a higher risk of death (Figures 3J, K).

Moreover, we incorporated an additional independent validation cohort, OEP000321, to verify the generalization ability of the models. In this rigorous external validation, the StepCox[forward] + RSF model achieved an AUC of 0.74, demonstrating the highest predictive accuracy among candidate models (Supplementary Figures 4A, B).

Instead of relying on a single metric, we evaluated the stability of the models over time. We calculated the average time-dependent AUC across all three datasets (training set, testing set, and the new OEP000321 validation set) for the 101 machine learning combinations. At the 1-year horizon, the StepCox[forward] + RSF model ranked 2nd among all combinations (Supplementary Figure 4C). At the 3-year horizon, the model ranked 1st, demonstrating superior mid-term predictive power (Supplementary Figure 4D). At the 5-year horizon, the model also ranked 1st, further confirming its long-term robustness (Supplementary Figure 4E).

We evaluated the prognostic value of the risk model in different subgroups of TCGA-LIHC using R software. The results showed that the model had strong prognostic significance across ten distinct subgroups stratified by age (>60 or ≤60 years), sex (female or male), pathologic stage (I, II, or III/IV), and T stage (T1, T2, or T3) (all p < 0.001; Figure 4). We then performed survival analysis across clinical subgroups within the entire ICGC-LIRI-JP cohort (Supplementary Figure 5) and validated the prognostic performance in the clinical subgroups of the entire GSE124751 cohort (Supplementary Figure 6). These findings suggest that this model has strong prognostic value in LIHC and is expected to be another new model for predicting the prognosis of LIHC patients.

We divided 781 LIHC patients into high-risk and low-risk groups based on the scoring of the risk model, compared the differential genes between the two groups in R software, and screened a total of 768 differential genes, including 367 up-regulated genes and 401 down-regulated genes (Figure 5A). To understand the potential enrichment pathways of these genes, we performed GO and KEGG enrichment analyses. The results of GO analysis showed that these differential genes were mainly enriched in the small molecule catabolic process, response to xenobiotic stimulus, and collagen-containing extracellular matrix (Figure 5B). KEGG analysis showed that these differential genes were mainly enriched in Complement and coagulation cascades, Cell cycle, and Retinol metabolism pathways (Figure 5C). GSEA enrichment analysis showed that the model was significantly associated with some pathways in LIHC, e.g., ADIPOGENESIS (NES = 1.836, p < 0.001), CHOLESTEROL HOMEOSTASIS (NES = 1.83, p < 0.001), COAGULATION (NES = 1.979, p < 0.001), FATTY ACID METABOLISM (NES = 1.918, p < 0.001), MTORC1_SIGNALING (NES = 1.791, p < 0.001), MYC_TARGETS_V1 (NES = 2.208, p < 0.001), OXIDATIVE PHOSPHORYLATION (NES = 2.225, p < 0.001), PEROXISOME (NES = 1.769, p < 0.001), PROTEIN SECRETION (NES = 1.75, p < 0.001), XENOBIOTIC METABOLISM (NES = 1.936, p < 0.001), and the model scores were significantly correlated with the core genes of the above pathways (ALDH2: rho = -0.376, p < 0.001; ADH4: rho = -0.397, p < 0.001; SERPING1: rho = -0.383, p < 0.001; SMS: rho=0.399, p < 0.001; EIF2S2: rho=0.401, p < 0.001; PA2G4: rho=0.433, p < 0.001; ACAT1: rho=-0.398, p < 0.001; TOP2A: rho=0.344, p < 0.001; CLTA: rho=0.41, p < 0.001; F11: rho=-0.392, p < 0.001) (Figures 5D–M). These results suggest that the modulation of downstream pathways by this model may be by targeting these molecules to produce positive or negative effects, providing new targets for subsequent anti-LIHC drug development.

We first analyzed the composition of immune cells in all samples using the CIBERSORT algorithm with R software. The results showed that the immune cells in the samples were predominantly T cells and macrophages, with the highest percentage of M0-type macrophages and the smallest percentage of neutrophils (Figure 6A). Analysis of the correlation between immune cells in the microenvironment revealed a significant positive correlation between CD4 memory activated T cells and CD8 T cells (Cor=0.36), as well as a significant negative correlation between CD4 follicular helper T cells and CD4 memory resting T cells (Cor=-0.41) (Figure 6B). Examination of the correlation between each gene in the model and tumor stroma score, immune score, ESTIMATEScore, and tumor purity demonstrated a strong correlation between HLA-B and B2M with the tumor stroma score, immune score, ESTIMATEScore, and tumor purity (Figure 6C). Further analysis indicated that a higher model score corresponded to greater tumor purity and lower tumor stroma score, immune score, and ESTIMATEScore (p < 0.01) (Figure 6D).

We analyzed the correlation between the risk model score and immune cell infiltration, immune checkpoint-related gene expression, chemokine and chemokine receptor infiltration with R software, and the results showed that the infiltration level of many kinds of immune cells in the microenvironment changed significantly with the higher risk score of the model, the expression level of immune checkpoint-related genes, such as MYD1, TNFRSF4, and other immune checkpoint-related genes, were significantly up-regulated. The infiltration levels of chemokines such as CXCL5, CXCL9, CXCL10, CXCL13, CXCL14, and chemokine receptors such as CCR2, CCR4, CCR5, CCR7, CCR9, CCR10, CXCR3, CXCR6, XCR1, etc., also changed significantly (p < 0.05) (Figures 6E–H).

To assess the clinical potential of risk models in immunotherapy, we analyzed the ICI responses of the high-risk and low-risk groups using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm with R software. The results showed that the ICB scores for CD274 and PDCD1LG2 were lower in the high-risk group, indicating that the high-risk group had a better response to targeted CD274 and PDCD1LG2 therapy (p < 0.001) (Figure 6I).

We analyzed the correlation between risk model scores and sensitivity to 25 clinically used chemotherapeutic drugs based on the GDSC2 database. The results demonstrated that the higher the risk model score, the lower the IC50 of Axitinib, Paclitaxel, Telomerase Inhibitor IX, and Vinblastine, indicating greater sensitivity to these drugs. Conversely, for Camptothecin, Carmustine, Cisplatin, Cyclophosphamide, Cytarabine, Dactinomycin, Entinostat, Erlotinib, Fludarabine, Foretinib, Gefitinib, Gemcitabine, Irinotecan, Niraparib, Oxaliplatin, Palbociclib, Pyridostatin, Rapamycin, Savolitinib, Talazoparib, and Venetoclax, higher IC50 values correlated with an increased susceptibility to develop resistance to these drugs (p < 0.05) (Figure 7). These findings suggest that risk model scores are instructive for personalized therapeutic drug selection in patients with LIHC.

Based on the somatic mutation data from the TCGA-LIHC dataset, we analyzed the mutational profiles associated with the risk model using R software. The results showed that missense mutation was the most common type of alteration (Figure 8A). The top 30 most frequently mutated genes in the entire cohort, the low-risk group (n = 149), and the high-risk group (n = 206) are presented in Figures 8B–D, respectively.

Tumor mutation burden (TMB) was defined as the total number of somatic gene coding errors, base substitutions, gene insertion or deletion errors detected per million bases. We analyzed the TMB scores between the high-risk and low-risk groups of the risk model and showed that the high-risk group had higher TMB scores (p < 0.05) (Figure 8E).

To systematically investigate the mutation frequency and spectrum of key genes in our model across diverse molecular contexts, we performed a direct comparison using cBioPortal between the model’s development cohort (TCGA-LIHC) and two other large-scale HCC studies (MSK 2024 and CLCA 2024). This comparative analysis revealed distinct molecular architectures among the three cohorts. The TCGA cohort was characterized by a high frequency of Copy Number Alterations (CNAs), predominantly gene amplifications. In contrast, the MSK 2024 cohort exhibited a mutation-driven profile, featuring specific alterations such as truncating mutations in HLA-B. The CLCA 2024 cohort displayed a relatively genomically inactive background with the lowest overall frequency of CNAs (Supplementary Figure 7).

We analyzed the correlation between the risk model score and ferroptosis-related genes, and the results showed that the higher the risk model score, the higher the expression levels of GPX4, AIFM2, SLC7A11, GSS, ACSL4, TMEM164, TFRC, FTH1, DHCR7, SLC3A2, SLC1A5, SLC39A7, HSPB1, PCBP2 were upregulated, while the expression levels of NFE2L2 and SLC40A1 were downregulated (p < 0 05), suggesting that the risk model may promote or inhibit ferroptosis by influencing the expression of these genes (p < 0 05). were also correspondingly up-regulated, while the expression levels of NFE2L2 and SLC40A1 were correspondingly down-regulated (p < 0.05), suggesting that the risk model may promote or inhibit ferroptosis by affecting the expression of these genes (Figure 9A).

We performed a correlation analysis between risk model scores and m6A methylation-related genes, and the results showed that the higher the risk model scores, the expression levels of METTL3, HNRNPC, RBM15, COPS5, DSN1, HNRNPA2B1, YTHDF2, IGF2BP2 and FABP5 were correspondingly up-regulated, whereas the expression levels of GMP6A, and ZC3H13 were correspondingly down-regulated (p < 0.05), suggesting that the risk model may promote or inhibit m6A methylation by affecting the expression of these genes (Figure 9B).

After analyzing the expression of 11 genes in the risk model in TCGA-LIHC, eight genes, HLA-B, HEXB, CDK4, NAA10, B2M, HSPD1, NPM1, and NEU1, were found to be highly expressed in the tumor tissues of TCGA-LIHC (Supplementary Figure 2). Analysis of the differential genes between the two groups, the high-risk group and the low-risk group, revealed that the upregulation of NEU1 was the most significant (Figure 5A), and prognostic analyses of these 11 genes in TCGA-LIHC revealed that NEU1 had a good prognostic function (Supplementary Figure 3). Therefore, we concluded that NEU1 may be a key risk factor in this model.

To clarify the role of NEU1 in cellular communication, we first analyzed its expression at the single-cell level using R. NEU1 was expressed in B cells, endothelial cells, fibroblasts, mesenchymal stromal cells, hepatocytes, monocytes/macrophages, plasma cells, and T/NK cells. Compared with non-tumor tissues, NEU1 expression levels in these cell types were significantly higher in LIHC tissues (p < 0.05) (Figures 10A–C). In LIHC tissues, the NEU1 high-expression group had a higher frequency and intensity of interactions compared to the NEU1 low-expression group (Figure 10D). Changes in the expression level of NEU1 could lead to alterations in the intensity of communication between a variety of cells, with the highest frequency of communication between endothelial cells (Figure 10E). Figure 10F demonstrates the potential mechanisms of communication, including ligand-receptor interactions and direction of communication, between different cells in the NEU1 high-expression and NEU1 low-expression groups. For example, in the NEU1 high-expression group, monocytes/macrophages may regulate changes in fibroblast/mesenchymal cell function or phenotype through SPP1 binding to CD44 or integrins (ITGAV+ITGB1).

NOTCH signaling pathway, TNF signaling pathway, VEGF signaling pathway, and PDGF signaling pathway are common histone acetylation-related pathways, and our study showed that these four pathways were most significant in endothelial cells in both NEU1 high- and low-expression groups (Figure 10G). Analysis of the interaction strength of histone acetylation-related pathways in different cells revealed that the interaction strength was highest in endothelial cells and low in fibroblasts/mesenchymal stromal cells when NEU1 was highly expressed, and it remained highest in endothelial cells and increased in fibroblasts/mesenchymal stromal cells when NEU1 was low-expression (Figure 10H).

We used the “CytoTRACE” package to quantify cell differentiation potential and compared the CytoTRACE scores of LIHC endothelial cells with those of normal endothelial cells. The results showed that the differentiation degree of LIHC endothelial cells was lower than that of normal endothelial cells (Figures 11A, B). Further analysis revealed that the expression level of NEU1 was significantly up-regulated in endothelial cells of LIHC tissues (Figure 11C) and that the endothelial cells of LIHC had a higher differentiation potential (Figure 11D). Among the endothelial cell subpopulations of LIHC, the expression level of NEU1 was significantly upregulated in low-differentiated endothelial cells of LIHC, which often had higher differentiation potential, which may be an important reason why the high expression of NEU1 suggests a poor prognosis for patients with LIHC (Figures 11E–H).

To reveal the trajectory of endothelial cell differentiation and the potential role of NEU1 in this process, we used Monocle3 to perform the proposed time-series analysis. The analysis showed that normal endothelial cells formed highly aggregated cell clusters and exhibited high stability, whereas endothelial cells of LIHC exhibited heterogeneity that was significantly different from that of normal endothelial cells (Figure 11I). Further analysis revealed that the expression level of NEU1 was significantly up-regulated mainly in endothelial cells of LIHC and some normal endothelial cells, whereas it was low-expressed in most normal endothelial cells (Figure 11J).

Pseudo-temporal trajectory analysis revealed a stable differentiation pathway of endothelial cells along the left-to-right direction, demonstrating a sequential process of transformation of normal endothelial cells to LIHC endothelial cells (Figure 11K). Notably, NEU1 expression levels showed progressive up-regulation on this differentiation trajectory, suggesting that NEU1 may play a key role in the transformation of normal endothelial cells to a malignant phenotype (Figure 11L).

Ten drugs highly related to NEU1 were screened through the Coremine database, from which the active compounds with structural information available in PubChem were selected and subsequently analyzed by molecular docking in the CB-Dock2 platform. The results showed that a variety of compounds exhibited good binding activities (Figure 12A), among which 2’-(4-methylumbelliferyl)-α-D-N-acetylneuraminic acid, N-Acetylneuraminic acid, N acetylneuraminolactose, 2-deoxy-2,3-dehydro-N-acetylneuraminic acid, and Oseltamivir all showed binding energy values below -5 kcal/mol, suggesting that these compounds have superior binding ability to NEU1 (Figures 12B–F).

We analyzed the expression of the aforementioned 11 genes in LIHC in clinical specimens using qRT-PCR. Compared with adjacent non-tumor tissues, mRNA levels of HLA-B, HEXB, NAA10, B2M, NPM1, and NEU1 were significantly upregulated in tumor tissues, while PON1 and ACAT1 were downregulated. No differential expression was observed for CFB, consistent with database analyses (Figures 13A–I). However, there was no statistically significant difference between CDK4 and HSPD1 in our sample (Supplementary Figure 8). Based on machine learning models and qPCR validation, NEU1 was selected for further investigation. Western blot analysis of LIHC in clinical specimens showed significantly elevated NEU1 protein expression in tumor tissues compared with paired adjacent non-tumor tissues (Figures 13J, K).

Similarly, NEU1 mRNA levels were significantly higher in LIHC cell lines (Hep3B and MHCC-97H) than in the normal hepatocyte line LO2 (Figure 14A). To investigate the functional role of NEU1, we transfected Hep3B and MHCC-97H cells with specific siRNA targeting NEU1. We found that the si-NEU1–3 fragment exhibited robust knockdown efficacy at both the mRNA and protein levels (Figures 14B–E) and was therefore selected for subsequent experiments. CCK-8 assays (Figure 14F), clone formation assays (Figures 14G, H), wound healing assays(Figures 14I, J), and Transwell invasion assays(Figures 14K, L) demonstrated that NEU1 knockdown significantly suppressed proliferation, migration, and invasion in Hep3B and MHCC-97H cells.