Integrated Machine Learning and Multi-Omics Identifies a Novel Molecular Signature for Improving the Prognosis of Hepatocellular Carcinoma

Introduction

The incidence of hepatocellular carcinoma (HCC) is on the rise, making it the third most prevalent cause of cancer-related mortality globally.^1,2 HCC typically presents with an insidious onset and rapid progression, contributing to a global five-year survival rate of less than 30%.³ Over the past decades, despite the establishment of a multidisciplinary treatment approach—including surgery, radiotherapy, radiofrequency ablation, molecular targeted therapy, and transcatheter arterial chemoembolization—morbidity and mortality have remained largely unchanged, and the overall prognosis remains poor.⁴ Given the significant heterogeneity of HCC, predictive, preventive, and personalized medicine (PPPM) represents a critical strategy for improving patient outcomes.⁵ Although several biomarkers related to HCC have been identified in recent years,^6–8 their clinical application remains limited. First, bulk transcriptomic analyses often fail to capture the pronounced heterogeneity of the tumor microenvironment (TME). Second, the robustness and generalizability of biomarker signatures developed across independent research groups show considerable variability in validation cohorts,^9,10 thereby diminishing their reliability for guiding individualized therapeutic decisions.

Currently, the intricate interactions between immune and cancer cells have emerged as a crucial focus of research efforts.^11,12 The TME consists of stromal cells, fibroblasts, immune cells, endothelial cells, and cancer cells.¹³ The significant association of the immune cells and various tumor cells in the TME have been demonstrated.^14–16 Cancer cells can maintain an immunosuppressive microenvironment and evade antitumor responses through multiple mechanisms, such as secreting immunosuppressive factors or recruiting regulatory T cells (Tregs).^17,18 Meanwhile, chemotherapy, radiotherapy, and targeted therapies have been shown to enhance the immunogenicity of tumor cells and stimulate immune responses, thereby improving the efficacy of anticancer treatments.^19–21 These findings suggest that immune-related genes (IRGs) serve as pivotal molecular mediators bridging tumor heterogeneity and personalized therapeutic strategies, underscoring their substantial translational and research significance.

With the advancement of sequencing technology, the integrated analysis of multi-omics has provided a new approach for exploring the biological characteristics of tumors. Herein, we employed a multi-omics approach combined with machine learning to investigate the potential clinical value of IRGs in HCC and identify an immune-derived molecular signature (IDMS). Our research indicates that the IDMS could be used as an innovative biomarker for predicting the prognosis and immunotherapy response in HCC, thereby contributing to improving patient outcomes through targeted prevention and personalized medical care.

Materials and Methods

Acquisition and Processing of Data

The single-cell dataset GSE162616²² was sourced from the Gene Expression Omnibus (GEO), converted into a Seurat object, and screened to retain high-quality cells (mitochondrial gene content < 20%, number of features > 300, log₁₀^GenesPerUMI > 0.80). The functions “Harmony” was employed to eliminate the impact of batch effects, followed by data normalization. The “FindMarkers” function was then employed to identify single-cell differentially expressed genes (scDEGs) with |logFC| > 0.25 and an adjusted p-value < 0.05. Bulk RNA-seq data were downloaded from The Cancer Genome Atlas (TCGA), and the associated clinical information was retrieved from the UCSC Xena platform. After excluding samples lacking clinical information, a final cohort of 424 patients from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) was obtained. To robustly identify key genes exhibiting consistent differential expression while minimizing the risk of overlooking biologically relevant molecules due to excessively stringent statistical thresholds, this study applied a DEG screening criterion of |log₂FC| > 0.3 and an adjusted p-value < 0.05. Additionally, three independent datasets (GSE14520, GSE76427, GSE121248) were obtained from GEO,^23–25 and a comprehensive dataset composed of 2,483 IRGs was retrieved from the ImmPort database. The full gene list is provided in Table 1.

Table 1 Immune-Related Genes (IRGs)

Weighted Gene Co-Expression Network Analysis (WGCNA)

The “WGCNA”²⁶ was utilized to compute correlation coefficients between gene pairs and apply these weighted coefficients to construct a gene co-expression network consistent with a scale-free topology. Subsequently, the TCGA-LIHC dataset was analyzed to calculate gene variance, and the top 5,000 genes with the highest median absolute deviation were selected. We then assessed correlations between HCC and control groups across various modules. Genes within modules that exhibited an absolute correlation coefficient (|r|) greater than 0.50 were selected for further analysis.

Identification of Immune-Derived Molecule Signature and Prognosis Model

We obtained hub genes by intersecting of IRGs, scDEGs, DEGs and module genes from WGCNA. Hub genes were selected by five machine learning algorithms: the least absolute shrinkage and selection operator (LASSO), Boruta, random forest (RF), learning vector quantization (LVQ), and Bagged Trees. The overlapping genes detected by these five algorithms were utilized to determine the IDMS for HCC. Subsequently, using GSE14520 as the validation set, six machine learning algorithms—LogitBoost, support vector machine (SVM), Naive Bayes (NB), RF, k-nearest neighbors (KNN), and AdaBoost—were applied to construct HCC classification models using the TCGA-LIHC dataset. Five-fold cross-validation was performed on each model to optimize para-meters. To guarantee the reliability of the results, the optimization process was repeated ten times, each time using different random seeds for the resampling. The model that yielded the highest area under the curve (AUC) on the GSE14520 validation set was selected as the optimal prognosis model. In addition, for assessing the predictive performance of the prognosis model, we compared the predictive capabilities of IDMS with those of previously reported HCC markers across the TCGA-LIHC, GSE76427, GSE121248, and GSE14520.

Survival Analysis and Establishing a Nomogram

To evaluate the prognostic efficacy of the model, the TCGA-LIHC dataset was stratified into high- and low-risk groups based on the median value of predicted values. Kaplan-Meier survival curves were identify differences in OS between the two groups.Time-dependent receiver operating characteristic (ROC) curves were generated to calculate AUC using “survivalROC”. Furthermore, clinical characteristics and the predicted value of the diagnostic model were incorporated into multivariate Cox regression analysis. Based on these factors, a nomogram was developed to assess the risk of HCC, and the predictive performance of the nomogram was evaluated using calibration curves.

Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

The “limma” and “clusterProfiler” packages were used to perform GSEA and GSVA on all genes in the TCGA-LIHC dataset. Functional enrichment differences between the high- and low-risk groups were then calculated. Genes with an adjusted p-value < 0.05 were considered statistically significant. All p-value adjustments were performed using the Benjamini-Hochberg method.

Analysis of Mutation and Immune Characteristics

SM data for HCC patients in the TCGA-LIHC dataset were processed using “VarScan”. CNV analysis was performed using GISTIC2.0.²⁷ The Mann–Whitney U-test was used to evaluate differences in microsatellite instability (MSI), tumor mutation burden (TMB), and tumor immune dysfunction and exclusion (TIDE) scores among the subgroups.

Analysis of Immune Infiltration and Drug Sensitivity

The abundance of immune cell types was evaluated using CIBERSORT. The Drug–Gene Interaction Database (DGIdb) was used to explore interactions between IDMS and small-molecule compounds for identifying potential therapeutic candidates. Predicted binding sites from CB-Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php) were used to visualize molecular docking between candidate drugs and their target proteins.

Quantitative Real‑time Polymerase Chain Reaction (qRT‑PCR)

Total RNA was isolated using the RNAiso plus RNA extraction kit (Takara,9019, Japan), following the manufacturer’s guidelines.The PrimeScript RT reagent kit (Perfect Real Time, Takara, RR037A, Japan) was employed to conduct reverse transcription. qRT‑PCR was conducted using the LineGene 9600Plus fluorescence qPCR system (Bioer, Hangzhou, China). Cycle threshold (Ct) values for each gene were normalized using GAPDH as an internal reference, and relative gene expression was calculated using the comparative Ct method. Data were derived from three biological replicates and analyzed employing Student’s t-test. The primer sequences and detailed cell lines are available in Tables 2 and 3.

Table 2 The Source of Cell Lines

Table 3 Primer Sequences for IDMS

Statistical Analysis

All data analyses were carried out with R software version 4.2.1 and the corresponding R packages. Student’s t-test was employed to compare data that followed a normal distribution,while the Mann–Whitney U-test was used for data with non-normal distribution.

Results

Research Process and Immune-Derived Genes in Single-Cell Transcriptome

The overall analysis workflow is shown in Figure 1A. Following quality control of scRNA-seq dataset, a total of 61,763 cells were retained. After being annotated with the “SingleR” package, we identified six distinct cell types and 590 scDEGs (Figure 1B and C).

Figure 1 Flowchart and clustering analysis of scRNA-seq (A) Overall analytical workflow. (B) UMAP plot illustrating six identified cell types. (C) Heatmap of marker genes for each cell type.

Identification of the Hub Genes in Bulk RNA-Seq

Using the “limma” R package, we identified 10,533 DEGs in the TCGA-LIHC dataset, including 9,175 upregulated and 1,358 downregulated genes (Figure 2A). Co-expression modules were constructed using WGCNA, and the optimal soft-thresholding power was determined to be 8 (Figure 2B). Correlations between module eigengenes and the HCC and control groups were evaluated based on the expression profiles of the top 5,000 genes with the highest median absolute deviation (MAD) (Figure 2C–E). And three modules containing genes were filtrated, namely: MEblue, MEblack, and MEbrown. Subsequently, seven hub genes—NR4A2, TMSB10, S100A10, NR4A1,MAP3K8, HSPA1A, and HSP90AA1—were identified for further investigation (Figure 2F).

Figure 2 Identification of the hub genes (A) Volcano plot of differential gene expression analysis in the TCGA-LIHC. (B) Selection of the optimal soft-thresholding power and assessment of network connectivity. (C) Module clustering results for the top 5,000 genes based on MAD. (D) Correlation between module eigengenes and HCC/control groups. (E) Gene dendrogram and module clustering results. (F) Venn diagram of IRGs, scDEGs, DEGs, and module genes.

Screening Immune-Derived Molecular Signature by Multi-Machine Learning

Five machine learning algorithms were applied to select the biological significance biomarkers in hub genes (Table 4). Seven genes (HSP90AA1, HSPA1A, MAP3K8, NR4A1, NR4A2, S100A10, TMSB10) were identified using LASSO algorithm (Figure 3A). Both the Boruta and LVQ algorithms selected seven genes, respectively (Figure 3B and C). Bagged Trees identified four genes (S100A10, NR4A1, HSP90AA1, TMSB10) (Figure 3D), while the random forest algorithm identified six genes (HSP90AA1, S100A10, NR4A1, NR4A2, TMSB10, MAP3K8) (Figure 3E). Finally, we selected four overlapping genes as IDMS: HSP90AA1, NR4A1, S100A10, and TMSB10 (Figure 3F).

Figure 3 Screening immune-derived genes (A) Gene frequency distribution derived from the LASSO algorithm. (B) Genes selected using LVQ algorithm.The red dotted line indicating the threshold of importance for gene screening. (C) Gene selection using the Boruta algorithm. (D) Variable importance and model performance from the Bagged Trees. (E) Variable contribution and model accuracy from the random forest model. (F) The Upset diagram shows the intersection of the screening results of five machine learning algorithms and determines the final key genes.

Table 4 The Parameters of the Five Machine Learning Algorithms for Feature Selection

Identifying Prognosis Model by Integrated Machine Learning Analysis

Based on 4 characteristic genes identified above, we evaluated the predict potential of IDMS using 6 machine learning algorithms. The ROC curve analysis showed that SVM model could effectively distinguish HCC from the control groups, with an AUC value of 0.957 (Figure 4A and B). Furthermore, compared with four previously published HCC signatures developed by Liu et al, Shen et al, Qian et al, and Zhou et al,^22,28–30 the IDMS exhibited consistently high predictive accuracy across multiple independent datasets, including TCGA-LIHC, GSE76427, and GSE14520. In the GSE121248 dataset, the predictive performance of IDMS was comparable to that of the biomarkers developed by Qian et al, with a higher meanAUC observed for IDMS (Figure 4C and D). Based on these findings, we selected the SVM as the optimal risk model for HCC. Subsequently, we performed annotation of IDMS and found that these genes were located on human chromosomes 1, 2, 12, and 14 (Figure 4E).

Figure 4 Identification of risk model and comparison with other HCC signatures (A) AUC values of risk models constructed using six machine learning algorithms. (B) ROC curve of the SVM model. (C) Circos plot comparing the performance of various published HCC signatures. (D) Heatmap showing predictive performance of the risk model across datasets. (E) Chromosomal mapping of IDMS.

Predictive Performance of SVM for HCC Prognosis

Using the median predicted value generated by the SVM model, a risk factor plot was created to visualize stratification between high- and low-risk groups (Figure 5A). The analysis showed that the expression levels of S100A10, TMSB10 and HSP90AA1 are significantly elevated in the high-risk group, and the mortality is also higher (Figure 5B). Kaplan-Meier survival analysis revealed that increased expression of these genes, along with the predicted value, was significantly associated with poor prognosis (Figure 5C–F). Moreover, time-dependent ROC analysis demonstrated that these markers had strong prognostic value at 1-, 3-, and 5-year intervals (Figure 5G–J).

Figure 5 Prognostic performance of the SVM for HCC (A) Risk factor plot within risk subgroups. (B) Compare gene expression levels among risk subgroups. (C-F) Kaplan–Meier survival curves for predicted value (C), HSP90AA1 (D), S100A10 (E), and TMSB10 (F). (G–J) Time-dependent ROC curves for predicted value (G), HSP90AA1 (H), S100A10 (I), and TMSB10 (J). ***p < 0.001.

Establishment of a Nomogram Integrating with Clinical Characteristics

Clinical characteristics from the TCGA-LIHC dataset were integrated with the predicted values of the SVM model for multivariate Cox regression analysis. Based on this analysis, a nomogram was constructed to estimate the survival probability of patients with HCC (Figure 6A). Among all the included variables, the predicted value of SVM demonstrated the strongest prognostic value, while the effect of gender was relatively weak. The calibration curve showed that the nomogram exhibited optimal clinical predictive performance at the 5-year (Figure 6B). Additionally, DCA demonstrated that the net clinical benefit of the prognostic model followed the pattern: 5-year > 3-year > 1-year (Figure 6C–E).

Figure 6 Establishment and evaluation of the nomogram (A) Nomogram integrating predicted value with clinical characteristics. (B) Calibration curves evaluating 1-, 3-, and 5-year survival predictions. (C–E) DCA of the prognostic model at 1-year (C), 3-year (D), and 5-year (E).

Underlying Molecular Mechanisms of IDMS

To investigate the biological processes underlying the prognostic relevance of IDMS, we performed GSEA and GSVA using transcriptomic data. The results revealed that oncogenic signaling pathways—including PI3K_AKT_MTOR, NOTCH, and WNT_BETA—were significantly activated in the high-risk group (Figure 7A–E and Table S1). In contrast, the low-risk group exhibited marked alterations in bile acid and fatty acid metabolism (Figure 7F and G and Table S2). These findings indicate a strong association between HCC risk stratification and dysregulated biological processes and metabolic pathways relevant to cancer progression.

Figure 7 Transcriptome characteristics of the subgroups (A) Ridge plot of biological functions identified by GSEA. (B–E) GSEA enrichment plots showing upregulation of hypoxia metagene (B), glycolysis/gluconeogenesis (C), fatty acid metabolism (D), and CDKN1A-mediated apoptosis via TP53 (E).(F and G) Heatmap (F) and box plot (G) showing GSVA results comparing risk subgroups.

Analysis of Mutation and Immune Characteristics in Different Subgroups

SMs in IDMS were analyzed and visualized using “maftools” (Figure 8A). The results showed two types of SMs, with missense mutations being the most prevalent. In addition, single nucleotide polymorphisms represented the predominant mutation type, and C >T transitions were the most frequently observed single nucleotide variants in HCC samples. CNVs were present in all IDMS-related genes in the TCGA-LIHC samples (Figure 8B and C). Among them, S100A10 showed the highest frequency of CNV amplification and the lowest frequency of CNV deletion, particularly in the variant group compared to the non-variant group. We further evaluated immunotherapy sensitivity in HCC samples by comparing risk groups classified by SVM model. No statistically significant were observed in TIDE scores or TMB among risk subgroups. In contrast, the MSI score was significantly higher in the high-risk group (Figure 8D and F). MSI is an established biomarker for predicting tumor response to immunological therapy.³¹ These findings indicate that patients in the high-risk group might have a greater likelihood of benefiting from immunotherapeutic strategies.

Figure 8 Mutation and immune characteristics in IDMS and its subgroups. (A) SM landscape of IDMS in the TCGA-LIHC. (B and C) CNV profiles of IDMS in TCGA-LIHC. (D–F) Comparison of immunotherapy response indicators between high- and low-risk groups: (D) TIDE scores, (E) MSI scores, and (F) TMB scores.

The Correlation Between the IDMS and Single-Cell Characteristics

A large amount of evidence shows an inseparable connection between the immune microen- vironment and HCC. In this study, the CIBERSORT algorithm was employed to quantify the abundance of 22 immune cell types and to compare their distribution across risk subgroups (Figure 9A). Additionally, a more in-depth analysis of immune cell infiltration uncovered intricate correlations among various immune cell populations. For instance, Plasma cells are positively correlated with naïve B cells, whereas monocytes are negatively correlated with follicular helper T cells (Figure 9B). Interestingly, Tregs exhibited a significant negative correlation with NR4A1, while displaying significant positive correlations with S100A10 and TMSB10. Moreover, HSP90AA1 demonstrated a significant negative association with M1-type macrophages, suggesting a potential role of IDMS in immune tolerance mechanisms (Figure 9C).

Figure 9 The correlation between the IDMS and single-cell characteristics (A) Stacked bar chart showing the distribution of immune cell types across risk subgroups. (B) Heatmap illustrating correlations among immune cell infiltration levels. (C) Heatmap showing correlations between immune cell abundance and IDMS expression levels.*p <0.05.

Analysis of Drug Sensitivity and Molecular Docking

To predict potential therapeutic agents for high-risk patients, we carried out a comprehensive drug sensitivity analysis using DGIdb and the DrugBank database to screen for drugs or small-molecule compounds interacting with IDMS (Figure 10A and Table 5). In addition, candidate therapeutic compounds and their target genes were evaluated using the Connectivity Map (CMap) (Figure 10B and Table 6). Based on compound susceptibility scores, the top two drugs were selected for molecular docking. The results indicated that HSP90AA1 exhibited moderate binding affinity with geldanamycin (Figure 10C). This interaction involved amino acid residues TYR493, GLU527, VAL530, GLN531, LEU533, LYS534, THR540, LEU541, VAL542, SER543, LYS546, THR603, TYR604, GLY605, TRP606, THR607, and MET610, which engaged in hydrogen bonding, weak hydrogen bonding, and hydrophobic interactions. Furthermore, HSP90AA1 also demonstrated moderate binding affinity with alvespimycin (Figure 10D). The interacting residues included TYR493, VAL530, GLN531, LYS534, THR540, LEU541, VAL542, SER543, LYS546, THR603, TYR604, and GLY605, forming similar hydrogen bonds and hydrophobic contacts. Finally, we evaluated the expression of IDMS in five cell lines and found HSP90AA1, NR4A1, and S100A10 significantly upregulated in tumor cell lines versus the liver cell line (Figure 10E–H).

Figure 10 Drug sensitivity and molecular docking analysis of IDMS (A) Interaction network between IDMS and candidate drugs or small-molecule compounds based on DGIdb and DrugBank analysis. (B) Predicted interactions between NR4A1, HSP90AA1, and candidate drugs based on CMap analysis. (C–D) Molecular docking results of HSP90AA1 with geldanamycin and alvespimycin, showing overall docking structure and interaction strength diagrams from left to right. The protein is color-coded from green (hydrophilic) to red (hydrophobic) to reflect amino acid properties. Hydrogen bonds are indicated by dashed blue lines; weak hydrogen bonds by light blue dashed lines; and hydrophobic interactions by gray dashed lines. (E–H) Validation of mRNA expression levels of HSP90AA1 (E), NR4A1 (F), TMSB10 (G), and S100A10 (H) in one liver cell line (THLE-2) and four HCC cell lines (Hep-3B, BEL-7405, MHCC97H, and HCCLM3) via qRT-PCR.

Table 5 List of Drugs Related IDMS From DGIdb Database

Table 6 List of Drug-Related IDMS and Drug Susceptibility Scores in CMap Database

Discussion

HCC remains a major malignancy that significantly threatens human health. Previous studies investigating the molecular basis of HCC have identified numerous potential biomarkers for diagnosis and prevention.^32–34 For instance, alpha-fetoprotein (AFP) is a commonly used biomarker for HCC diagnosis, however, its limited sensitivity and specificity restrict its clinical utility.³⁵ In addition, vascular endothelial growth factor A has demonstrated potential for predicting the efficacy of sorafenib treatment.³⁶ Promising molecular markers have also been identified in circulating tumor DNA (ctDNA) and circulating tumor cells.³⁷ Nevertheless, due to a lack of large-scale prospective studies, robust molecular biomarkers for clinical decision-making in HCC remain limited. At present, pathological classification remains a primary tool for determining treatment strategies. However, substantial differences in patient prognosis are frequently observed under this model. Therefore, in the context of HCC—characterized by limited diagnostic and therapeutic options and marked heterogeneity—identifying reliable molecular biomarkers to support the implementation of PPPM is a critical and urgent need.

Advancements in multi-omics technologies have facilitated the elucidation of tumor molecular mechanisms.³⁸ In this work, we integrated bulk and single-cell transcriptome datasets to explore potential biomarkers and characteristics of immune cell infiltration in HCC. To minimize bias stemming from a single algorithm, this study applied five machine learning algorithms with cross-validation to identify four immune-derived genes that exhibit greater stability and biological significance. Subsequently, we comprehensively analyzed the predictive performance of six machine learning approaches and selected the optimal risk model. Compared with the previous immune-related prognostic models constructed based on single-cell transcriptomes and bulk RNA-seq,^39,40 through cross-validation with multiple algorithms and independent validation across multiple cohorts, immune-derived molecular features with higher generalization ability were obtained.More importantly, survival analysis and nomogram indicated that IDMS exhibits a robust capacity for risk stratification and prognostic value. This result indicates that IDMS not only has better predictive performance in the prognosis of HCC, but can also more effectively identify high-risk patients and guiding the development of targeted surveillance and intervention strategies. For high-risk individuals, clinicians may consider increasing monitoring frequency and recommending lifestyle interventions such as routine liver ultrasounds, smoking cessation, and reduced alcohol intake. Furthermore, detecting ctDNA can enhance the early diagnosis of HCC and prevents its occurrence.⁴¹

In the past few years, immunotherapy has become a major focus of basic and clinical oncology research, with promising advances in its application to HCC.^42,43 However, challenges persist, including low objective response rates and adverse treatment effects.⁴⁴ Therefore, identifying patients most likely to benefit from immunotherapy is a key strategy for improving outcomes and personalizing treatment plans. Multi-omics and drug sensitivity analysis in this study indicated that patients in the high-risk group were more likely to benefit, potentially reducing overall treatment costs. Moreover, geldanamycin and alvespimycin may exhibit superior therapeutic efficacy. Geldanamycin inhibits the function of HSP90AA1, leading to the degradation of multiple tumor-related proteins, thereby preventing the proliferation of tumor cells and inducing apoptosis.^45,46 Alvespimycin inhibits tumor proliferation, metastasis and angiogenesis by suppressing pathways such as PI3K/AKT and MAPK, and enhances the sensitivity to chemotherapy.^47–49Taken together, these findings suggest that IDMS are superior to previous markers in the prediction and personalized treatment of HCC, and are more conducive to the clinical management of PPPM for HCC patients.

Through the analysis of heterogeneity of IDMS at the Inter-cellular level, we found that all genes included in the IDMS are implicated in immune evasion or tumorigenesis. Tregs are a subset of T cells with immunosuppressive functions. Within the tumor microenvironment (TME) of HCC, persistent infiltration of Tregs promotes an immunosuppressive niche, thereby facilitating tumor progression and immune escape.^18,50,51 NR4A1, a transcription factor, has been shown to promote the differentiation of naïve T cells into Tregs and enhance immunosuppression by modulating signaling pathways downstream of the T cell receptor (TCR).^52,53 S100A10 has been reported by Wang et al to activate the cPLA2 and 5-LOX axis, leading to CD8⁺ T cell exhaustion in HCC and promoting sustained tumor progression and metastasis.⁵⁴ TMSB10, a member of the β-thymosin family, is overexpressed in various cancers.^55,56 Knockdown of TMSB10 has been shown to facilitate the polarization of macrophages towards the M1 phenotype, thereby improving antitumor immunity.⁵⁷ HSP90AA1, a molecular chaperone, interacts with numerous client proteins and is involved in multiple biological processes related to cancer. Notably, plasma levels of HSP90AA1 are significantly elevated in HCC patients and exhibit superior diagnostic performance compared to AFP, with strong associations with tumor progression and metastasis.^58,59 Importantly, both S100A10 and TMSB10 are positively correlated with Treg abundance and are highly expressed in HCC tissues, suggesting that they may play key roles in mediating Treg-associated immune escape mechanisms. However, the molecular mechanisms involved are not yet fully elucidated. Moreover, the potential biases that may arise from the publicly available data also require further investigation.

We investigated the underlying mechanisms of the IDMS using multi-omics analysis. DEGs between the risk subgroups were primarily correlated with cell signaling pathways—such as the PI3K, WNT, and MYC pathways—as well as metabolic processes including glycolysis, xenobiotic metabolism, and fatty acid metabolism. These results may partially explain the differences in prognosis observed between IDMS-defined subgroups. Furthermore, metabolism-associated biological pathways may represent promising directions for future research, and the development of targeted inhibitors against these pathways could offer therapeutic benefit in preventing HCC progression.

Nevertheless, there are certain limitations in this investigation. First of all, all data were sourced from publicly databases. Although the immune-derived molecular signature evaluated in this study demonstrated consistent prognostic predictive performance across multiple independent datasets, their clinical utility requires further validation in large-scale, prospective cohort studies. Additionally, our findings are derived from computational bioinformatics analyses, and functional experimental validation is lacking to confirm the biological roles of the genes comprising the IDMS. We hope that the publication of this study will attract broader research interest in the IDMS and contribute to a deeper understanding of its function in the pathogenesis and progression of HCC.

Conclusion

We identified an IDMS that shows potential as a valuable tool for predictive, preventive, and personalized medicine in HCC. Through multi-omics analysis, this study may offer novel perspectives on the molecular mechanisms that underlie the occurrence and development of HCC.