Establishment and evaluation cuproptosis-related gene signature for predicting the prognosis and immunotherapy response of hepatocellular carcinoma

Abstract

Background

This study aims to develop a novel cuproptosis-related model through bioinformatics analysis, providing new insights into HCC classification. It also explores the correlation between the cuproptosis-related risk score and factors such as prognosis, tumor mutation burden (TMB), biological function, tumor microenvironment (TME), and immune efficacy.

Methods

We performed unsupervised clustering of cuproptosis-related gene expression profiles from TCGA and GEO to identify molecular subtypes and differentially expressed genes. Prognostic models were constructed using univariate, Lasso, and multivariate Cox regression analyses. HCC patients were classified into high-risk and low-risk subgroups, and the model’s prognostic value was assessed through survival analysis, ROC curves, and nomograms. Immune checkpoint, drug sensitivity, and IPS were used to evaluate immunotherapy response. The model’s predictive ability was further validated with the ICGC database and IMvigor210 cohort. Finally, key gene expression and biological functions were validated in human tissues and HCC cell lines.

Results

The cuproptosis-related gene risk score model (CRGRM), based on GMPS, DNAJC6, BAMBI, MPZL2, ASPHD1, IL7R, EPO, BBOX1, and CXCL9, independently predicted HCC prognosis and immune response. Clinical correlation and ROC curve analysis demonstrated its accuracy in predicting 0.5-, 1-, 3-, and 5-year survival. The risk score also strongly correlates with immunotherapy response and serves as a reliable treatment predictor. Drug sensitivity analysis revealed that the low-risk group was more sensitive to dasatinib, imatinib, and gefitinib. In vitro, BAMBI knockdown significantly inhibited HCC cell proliferation and metastasis.

Conclusions

This model demonstrates potential in predicting prognosis and immunotherapy response, providing insights into personalized treatment strategies for HCC. Additionally, our study suggests that BAMBI may serve as a novel biomarker and potential therapeutic target for HCC.

Introduction

Primary liver cancer, including hepatocellular carcinoma (HCC), is a leading malignancy in the digestive system [1], ranking sixth in new cases and second in deaths worldwide [2]. HCC accounts for 85–90% of primary liver cancers [3] and is strongly linked to HBV infection, chronic alcohol use, and family history [4]. Most patients are diagnosed at an advanced stage, making immunotherapy crucial. Our study identifies novel immunotherapy targets for HCC and evaluates their potential efficacy.

In addition to traditional forms of cell death such as ferroptosis, autophagy, and pyroptosis, cuproptosis has emerged as a crucial mechanism in tumor progression [5]. Copper, a trace element, is vital for numerous biological processes. Recent studies have demonstrated that copper induces cell death by binding to fatty acylation proteins in the tricarboxylic acid cycle, thereby linking mitochondria to the process of cell death [6]. By disrupting this cycle, copper induces proteotoxic stress, which accelerates programmed cell death. Moreover, mitochondrial dysfunction is closely associated with drug resistance in cancer, contributing to a poor prognosis [7]. This further underscores the profound impact of cuproptosis on tumor progression, highlighting the need for continued exploration.

Copper (Cu) levels are elevated in hepatocellular carcinoma (HCC) patients, and free serum copper is linked to the incidence and progression of various cancers, including HCC [8]. Studies show that COMMD10 disrupts Cu-Fe homeostasis in HCC, enhancing radiosensitivity through the HIF1α/CP loop [9]. Immunotherapy has become a key treatment for liver cancer, improving overall survival, yet the role of cuproptosis in HCC progression remains unclear, particularly its potential in guiding immunotherapy. Copper ions, essential for immune function, influence immunity and immunotherapy response. This study aims to explore the relationship between cuproptosis-related genes, HCC prognosis, and immune response using multi-omics analysis from public databases.

We developed a cuproptosis-related gene risk model (CRGRM) comprising nine genes. ROC curve and prognostic nomogram analyses demonstrated its strong predictive value for HCC prognosis. Additionally, we explored the model’s correlation with immunotherapy by assessing its association with TMB, immune checkpoints, and drug sensitivity.

Materials and methods

Data collection and gene set acquisition

We obtained hepatocellular carcinoma (HCC) transcriptome data (424 samples), clinical data, and somatic mutations from the TCGA database (https://portal.gdc.cancer.gov/repository), including 370 samples with complete clinical information. Expression and clinical data were also collected from the GSE76427 cohort in the GEO database (https://www.ncbi.nlm.nih.gov/gds/). We then merged the data from the two platforms using the “limma” and “sva” packages. Based on previous literature, we combined the cuproptosis-related gene set with the MsigDB v7.0 database (http://www.gsea-msigdb.org/gsea/msigdb/), identifying 49 relevant genes after removing duplicates [6, 10]. External validation was performed using gene expression and clinical data from the GSE14520 cohort and the International Cancer Genome Consortium (ICGC) LIRI-JP project (https://dcc.icgc.org/projects/LIRI-JP). Additionally, transcriptome and clinical data from the IMvigor210 dataset were used to validate the risk model's predictive efficacy for immunotherapy.

Identification of differentially expressed and prognostic genes

Differential gene expression between hepatocellular carcinoma and normal tissues was analyzed using the Wilcoxon rank sum test, based on data from the TCGA dataset and processed with the “limma” package. Univariate Cox regression and Kaplan–Meier (KM) survival analyses were employed to identify genes associated with prognosis.

Unsupervised cluster analysis

Based on the expression levels of cuproptosis-related genes, consensus clustering analysis, cumulative distribution function (CDF), and consensus matrix were used to determine the optimal number of clusters.

Analysis of immune traits

Single-sample gene set enrichment analysis (ssGSEA) was performed to assess differences in immune function pathways and immune cell infiltration between groups [4]. The “ESTIMATE” R package was used to calculate stromal, immune, and estimated scores, as well as tumor purity, within the tumor microenvironment (TME), and compare these scores across subtypes. Immune correlations of model genes and total risk scores were analyzed using the CIBERSORT algorithm. Additionally, the expression of common immune checkpoints (ICI), m6A methylation, and human leukocyte antigen (HLA) genes were compared among different subgroups.

Functional analysis of differentially expressed genes

GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analyses evaluate the biological functions and pathways of differentially expressed genes (DEGs). GSVA (Gene Set Variation Analysis) assesses gene set enrichment across samples, offering a comprehensive view of biological differences. GSEA (Gene Set Enrichment Analysis) determines whether predefined gene sets show significant expression differences. Both GSVA and GSEA were conducted using the “clusterProfiler” and “org.Hs.eg.db” R packages.

Evaluation of immunologic efficacy

The immunophenoscore (IPS) data for HCC samples were obtained from the TCIA-LIHC project (https://tcia.at/patients) to predict patient response to biological agents such as CTLA4 and PD-1 inhibitors. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was used to predict TIDE scores and analyze score differences across risk groups, aiding in the assessment of immunotherapy sensitivity. Additionally, we explored the relationship between risk scores and the six immune subtypes described in previous literature.

Establishment and validation of prognostic models

We analyzed differentially expressed genes across clusters based on the optimal grouping of cuproptosis-related genes. Univariate regression, LASSO Cox, and multivariate Cox regression analyses were employed to identify prognostically significant genes and construct a risk model by assigning risk coefficients. HCC patients were stratified into high- and low-risk groups based on the median risk score from the training cohort. Risk curves were visualized using the “pheatmap” R package. Survival analysis revealed significant prognosis differences between the groups, and a time-dependent ROC curve assessed the model's sensitivity and specificity. Model accuracy was further validated using C-index analysis.

Analysis of the prognostic value of the model and construction of prognostic nomogram

Univariate and multivariate regression analyses were performed to identify independent prognostic predictors for HCC. Risk score, clinical stage, gender, and age were incorporated to construct a prognostic nomogram, which estimated the probability of patient outcomes at 0.5, 1, 3, and 5 years. Prediction accuracy was validated using a calibration curve.

Clinical relevance analysis

To assess the clinical relevance of the CRGRM, we conducted an association analysis of key clinical features using the chi-square test and visualized the results with boxplots and a heatmap.

Waterfall plot of mutated genes

Waterfall plots of high-risk and low-risk groups were established by “maftools” R package to analyze the differences of potential tumor mutations.

Prediction of drug susceptibility

Drug sensitivity analysis was performed by using pRRophetic and ggplot2 software to estimate the IC50 value of each HCC sample and compare the differences in drug sensitivity between high and low risk groups.

Tissue specimens and immunohistochemical staining

Tissue samples were collected from 36 HCC patients at the Affiliated Hospital of Xuzhou Medical University, with approval from the hospital’s Ethics Committee (Ethics approval: XYFY2022-KL486-01). Paraffin-embedded tissue sections were baked at 65 °C for 2 h, deparaffinized with xylene, and hydrated with gradient ethanol. Antigen retrieval was performed using citrate solution, followed by blocking with 10% goat serum. Sections were incubated overnight with anti-BAMBI (1:200) at 4 °C, washed with PBS, and incubated with secondary antibodies for 20 min at room temperature. DAB staining was applied for 3 min, followed by hematoxylin re-staining for 1 min. Images were captured after resin sealing and drying.

Cell culture

HepG2 cells (Pricella, China) and MHCC-97H cells (Beyotime, China) were verified by short tandem repeat (STR). HepG2 cells were cultured in MEM with glutamine and 10% FBS, while MHCC-97H cells were cultured in DMEM with glutamine and 10% FBS.

Real-time fluorescent quantitative PCR

Total RNA was extracted from tissues and cells using the FastPure R Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, China), and reverse transcription quantitative PCR (RT-PCR) was performed with the SYBR Premix Ex Taq™ II kit (Vazyme, China). The reaction conditions included 40 cycles at 95 °C for 1 min, with 10 s at 95 °C and 40 s at 57 °C. Relative expression levels were calculated using the 2-ΔΔCt method, with GAPDH as the endogenous control. Primer sequences are provided in Supplementary Table 1.

Western blot

Proteins were extracted from tissues and cells using RIPA and PMSF, with GAPDH as the loading control. Protein signals were detected using ECL reagents. Primary antibodies used included anti-BAMBI, NRP1, VTCN1 (B7H4), and GAPDH (all from Abcam). BAMBI interference was performed using a viral vector from Zebrafish Biotech, with shRNA sequences: sh-BAMBI#1: 5ʹ-GGG CAC GGG AAG CTG GAA T-3ʹ and sh-BAMBI#2: 5ʹ-GCG CCT GCT TCT CTA GAC T-3ʹ.

Cell counting kit-8 assay

HepG2 and MHCC-97H cells from both transfection and control groups were seeded in 96-well plates and incubated. Cell viability was measured at specified time points using the CCK8 kit (GLPBIO, USA).

Clone formation test

Cells stably transfected with control and shRNA were seeded into small culture dishes. After 2 weeks, colonies were stained with crystal violet.

Edu assay

Transfected control cells and stable shRNA-transfected cells were seeded in 96-well plates for 24 h. EdU solution from the EdU kit (Ribobio, China) was added, and cells were incubated for 2 h. After fixation with 4% paraformaldehyde in PBS for 30 min, staining was performed according to the manufacturer’s protocol. Images were captured using a Nikon microscope, and positive cells were quantified using ImageJ software.

Wound healing assay

Cell migration was evaluated using a wound healing assay. Transfected control cells, sh-BAMBI#1, and sh-BAMBI#2 stable cells were seeded in 6-well plates. Once cells reached 100% confluence, a wound was made in the center using a 200 μL pipette. After washing with PBS, serum-free medium was added, and cells were cultured in the incubator. Images were captured at 0, 12, and 48 h after scratching, and wound areas were quantified using ImageJ software.

Cell invasion assay

Knockdown and control cells were resuspended in serum-free medium, seeded into 4 × 6 Transwell chambers, and cultured on Matrigel-coated upper chambers. The lower chamber contained 12% FBS medium. After 24 h of incubation, cells adhering to the lower membrane were fixed with 4% paraformaldehyde, stained with crystal violet, and the upper membrane was wiped to remove non-adherent cells. Images were captured using a light microscope, and cell counts were performed using ImageJ.

Statistical analysis

All analyses were conducted using R software (version 4.1.3, https://www.r-project.org/). Student’s t-test was used for comparing two groups of continuous, normally distributed variables, while the Wilcoxon rank-sum test was applied for multiple ordered groups. Statistical significance was defined as p < 0.05, with * indicating p < 0.05, ** p < 0.01, and *** p < 0.001.

Results

Data merging and differential expression of cuproptosis-related genes (CRG)

Gene expression data from the TCGA-LIHC and GSE76427 datasets were integrated using the “limma” and “sva” R packages to eliminate batch effects and standardize profiles. Forty-six cuproptosis-related genes were selected for further analysis. Comparison of their expression levels between tumor and normal tissues revealed significant differences (|log2 FC|> 0.585, all P < 0.05) (Fig. 1A). Univariate Cox regression identified 20 genes, including LIAS, FDX1, SLC31A1, ATOX1, and AQP1, as favorable prognostic factors for HCC, while the other 29 were linked to poor prognosis (Fig. 1B). Survival analysis confirmed that 34 cuproptosis-related genes were significantly associated with HCC prognosis (Supplementary Fig. 1–2). Additionally, protein–protein interaction (PPI) networks were constructed to explore gene interactions involved in cuproptosis -related processes, providing insights into the molecular mechanisms of cuproptosis in HCC (Supplementary Fig. 3).

Fig. 1

Data merging and differential expression of cuproptosis-related genes (CRG). A Differential expression of CRGs in hepatocellular carcinoma and normal tissues. B The prognostic network map

Classification of HCC subtypes and exploration of biological functions

Unsupervised consensus clustering of HCC samples based on 46 cuproptosis-related genes identified two subtypes, Cluster A and Cluster B, with the most consistent clustering method selected (Fig. 2A). Differences in biological functions and immune characteristics between the subtypes were explored. Principal component analysis (PCA) clearly distinguished the two subtypes (Fig. 2B), and survival analysis showed a significant difference in overall survival (Fig. 2C). A heatmap of clinicopathological features revealed notable differences in tumor stage and histologic grade between the clusters (Supplementary Fig. 4). Gene set variation analysis (GSVA) evaluated the pathway enrichment differences between the two subtypes (Supplementary Fig. 5). Differential gene expression analysis identified 2474 DEGs (Fig. 2D), and subsequent GO and KEGG pathway analyses revealed functional enrichment (Fig. 2E, F).

Fig. 2

Classification of HCC subtypes and exploration of biological functions. A Unsupervised consensus clustering analysis. B The PCA plots showed the cluster could clearly distinguish HCC patients into two groups based on their expression profiles. C The Kaplan–Meier curve between the two clusters. D The Venn diagram was utilized to display DEGs between cluster A and cluster B. E GO pathway enrichment of differentially expressed genes between cluster A and cluster B. F KEGG pathway enrichment of differentially expressed genes between cluster A and cluster B.*P < 0.05; **P < 0.01

Immune profiling analysis based on consensus cluster analysis

In this study, the ESTIMATE algorithm was used to calculate immune scores, stromal scores, estimate scores, and tumor purity within the tumor microenvironment (TME). Box plots visualized differences in immune microenvironment characteristics across subtypes, revealing that cluster A had significantly higher immune, stromal, and estimate scores compared to cluster B (Fig. 3A), but lower tumor purity (Fig. 3B). Single-sample Gene Set Enrichment Analysis (ssGSEA) further confirmed higher immune cell enrichment in cluster A (Fig. 3C). Additionally, significant differences in immune functions were observed between the clusters, with Type I and II interferon responses enriched in cluster B, and other immune functions more prominent in cluster A (Supplementary Fig. 6). To assess potential immunotherapy responses, we compared the expression of immune checkpoint inhibitors (ICI), HLA, and m6A-related genes. The results showed higher expression of these genes in cluster A, suggesting that patients in this group may have a better response to immunotherapy (Fig. 3D–F).

Fig. 3

Immune function analysis among different clusters. A The box plots to show the changes of immune microenvironment in different clusters. B Violin plots shows differences in tumor purity between clusters. C The boxplot displaying immune cell enrichment status between different cluster groups. The expression levels of immune checkpoint genes (D), HLA genes (E) and M6A genes (F) in different clusters were compared

Consensus clustering analysis of prognostic genes (uniSigGenes) and establishment of CRGRM

To explore the potential of cuproptosis-related genes (CRGs) in hepatocellular carcinoma (HCC), we conducted univariate Cox regression analysis on differentially expressed genes (DEGs) from previous clusters, identifying a prognosis-related gene expression profile (uniSigExp) with a significance threshold of P < 0.05. A total of 485 HCC patients were randomly divided into a training set (n = 243) and a validation set (n = 242), with clinical information summarized in Table 1. Cluster analysis of uniSigGenes revealed two subtypes (Cluster A and Cluster B) with the highest consistency (Fig. 4A, B). Survival analysis indicated that Cluster A had significantly worse overall survival (OS) than Cluster B (Fig. 4C). Box plots showed differential expression of cuproptosis-related genes between the two clusters (Supplementary Fig. 7). A combined heatmap displayed the expression of DEGs alongside clinical features (age, gender, T stage, clinical stage, histologic grade) in both clusters (Supplementary Fig. 8). LASSO Cox regression analysis identified 19 prognostic genes in the training set while minimizing overfitting (Fig. 4D, E). Multivariate Cox regression refined the model to 9 DEGs most associated with prognosis, which were used to construct a cuproptosis-related risk model (CRGRM). Based on these genes, we calculated a cuproptosis-related risk score (CRS) for each patient, stratifying them into high- and low-risk groups based on the median score. The construction of the CRGRM and survival distribution across subtypes was visualized using a Sankey plot (Fig. 4F).

Table 1 Clinical information of the patients with HCC and the comparability of the groups

Fig. 4

Consensus clustering analysis of prognostic genes (uniSigGenes) and establishment of CRGRM. A Concordance matrix of subtypes. B The cumulative distribution function based on the sign indicated that the optimal number of subtypes was 2. C The Kaplan–Meier curve for cluster A and cluster B. D, E LASSO Cox regression analysis based on uniSigGenes for screening the prognostic model genes. F The sankey diagram illustrates the construction process of CRGRM and the distribution of survival status among different subtypes

Validation of cuproptosis-related genes risk model (CRGRM)

Box plots were used to display the differences in risk score distributions across CRG clusters (Fig. 5A, B). Survival analysis revealed significantly worse overall survival in the high-risk group compared to the low-risk group (Fig. 5C–E). The risk curve further indicated that higher risk scores were associated with increased mortality and shorter survival times (Fig. 5F–H). Differential expression of CRGs between high-risk and low-risk groups was shown in a box plot (Fig. 6A). The risk score demonstrated high specificity and sensitivity for predicting HCC prognosis, with AUC values for 0.5, 1, 3, and 5-year survival in the entire cohort of 0.792, 0.801, 0.749, and 0.712, respectively (Fig. 6B). To validate the model, we assessed its predictive accuracy in external datasets. In the GSE14520 cohort, the AUC for overall survival at 0.5, 1, 3, and 5 years were 0.645, 0.662, 0.639, and 0.604, respectively (Fig. 6C). In the ICGC cohort, the AUC values for the same time points were 0.717, 0.775, 0.670, and 0.798, respectively (Fig. 6D). The classification ROC curve, which included other clinical characteristics and the risk score (RS), showed that the RS AUC was significantly higher than for age, gender, and tumor stage, indicating superior predictive power of the CRGRM (Fig. 6E). Finally, the c-index was calculated to evaluate model discrimination, and the results showed that the CRGRM outperformed other clinical factors in prediction accuracy (Fig. 6F).

Fig. 5

Validation of cuproptosis-related genes risk model. A, B The boxplot presents the distribution of risk scores for HCC patients based on different GRG clusters. The survival curves based on entire group (C), GSE14520 cohort (D), and ICGC cohort (E). The risk curve of all group (F, I, L), Training group (G, J, M), and validation group (H, K, N)

Fig. 6

Evaluation of the prognostic value of CRGRM. A The boxplot to reveal the differential expression of CRG in the high-risk and low-risk groups. B–D The ROC curve of 0.5, 1, 3, and 5 years in the HCC patients. C The ROC curve of 0.5, 1, 3, and 5 years in the GSE14520 cohort. D The ROC curve of 0.5, 1, 3, and 5 years in the ICGC cohort. E The classification ROC curve were used to evaluate the predictive probability of OS for HCC patients with factors including risk score, sex, age, and tumor stage. F C-index analysis showed that the risk score had high prediction accuracy

Clinical relevance, independent prognostic value and nomogram construction of CRGRM

The risk score was significantly correlated with the clinicopathological characteristics of HCC patients, including gender, age, T stage, histological grade, survival status, and clinical stage (Supplementary Fig. 9). In external validation, the risk score showed positive correlations with AFP level, CLIP grade, tumor size, and stage in the GSE14520 cohort (Supplementary Fig. 10), and with survival status and clinical stage in the ICGC cohort (Supplementary Fig. 11). To assess the clinical applicability of the risk model, survival analysis was performed within different clinical subgroups of HCC patients (Fig. 7), revealing that the low-risk group consistently had better overall survival (OS) compared to the high-risk group. Univariate and multivariate Cox regression analyses confirmed that clinical stage and risk score were independent predictors of HCC prognosis (Fig. 8A). In both GSE14520 and ICGC datasets, the risk score remained a significant independent prognostic factor (Fig. 8B, C). Based on these analyses, we constructed a nomogram incorporating gender, age, clinical stage, and risk score to predict 1-, 3-, and 5-year survival probabilities (Fig. 8D). The calibration curve indicated high consistency between the predicted and actual survival values, demonstrating the nomogram’s reliability for clinical prognosis prediction (Fig. 8E).

Fig. 7

Clinical relevance of cuproptosis-related prognostic models. A Boxplot reflecting the risk scores of HCC patients by different gender. B Boxplot of risk score distribution based on age of HCC patients. C Boxplot of risk scores for HCC patients at different stages. D Boxplots of risk scores for HCC patients with different tissue grades. E Boxplot of risk scores for HCC patients based on different survival status. F Boxplot of risk scores for HCC patients at different clinical stages. G–I The K–M survival curve based on different clinical subgroups

Fig. 8

Independent prognostic value and nomogram construction of CRGRM. A, B Univariate and multivariate Cox regression analyses demonstrated that clinical stage and risk score were independent predictors of prognosis in HCC patients. C, D The forest plot revealed that risk score was also an independent risk factor for HCC prognosis in the GSE14520. E, F The forest plot revealed that risk score was also an independent risk factor for HCC prognosis in the ICGC cohorts. G A nomogram was used to predict the 1 -, 3 -, and 5-year survival probabilities of HCC patients. H The calibration curve proved that the nomogram had high accuracy in predicting the prognosis of patients

Analysis of immunologic properties of CRGRM

To assess immunotherapy efficacy, we examined immune traits across risk groups, focusing on immune scores, immune cell infiltration, and immune microenvironment changes. First, the CIBERSORT algorithm was applied to analyze the infiltration of 23 immune cell types in each HCC sample. We then explored the correlation between model genes, overall risk score, and immune cell infiltration, visualizing the results through a correlation heatmap (Fig. 9A) and scatter plot (Supplementary Fig. 12). Next, ESTIMATE algorithm results showed that the immune score, estimated score, and stromal score were significantly higher in the low-risk group, suggesting better immunotherapy response and higher sensitivity in this group (Fig. 9B). Tumor purity differences between risk groups were also analyzed (Fig. 9C). Finally, ssGSEA analysis revealed that the low-risk group exhibited greater immune cell enrichment and stronger tumor immune responses, indicating higher immunotherapy potential (Fig. 9D).

Fig. 9

Analysis of immunologic properties of CRGRM. A The correlation heatmap reflecting the correlation between model genes and risk scores and immune cells. B Violin plot shows the differences in stromal score, immune score,and estimate score of HCC patients in different risk groups. C Violin plot shows the differences in tumor purity of HCC patients in different risk groups. D Comparison of immune cell enrichment in low-risk and high-risk HCC patients

The potential role of CRGRM in immunotherapy

We used correlation heatmaps to explore associations between risk scores and immune checkpoint gene expression (Supplementary Fig. 13). Boxplots revealed that BTLA, LAG3, CD27, ICOS, LILRB2, and IDO1 were significantly upregulated in low-risk HCC patients, while NRP1 and VTCN1 were elevated in high-risk patients (Fig. 10A). Notably, NRP1 and VTCN1 expression was positively correlated with BAMBI, as confirmed by TIMER2.0 (Supplementary Fig. 14). ssGSEA analysis demonstrated significant differences in 11 of 13 immune functions between the two risk groups, all enriched in the low-risk group (Fig. 10B). The low-risk group also exhibited a higher immunophenoscore (IPS) and better response to immunotherapy, indicating the model’s ability to predict immune checkpoint inhibitor (ICI) response (Fig. 10C). Furthermore, HLA gene expression was significantly higher, while M6A genes were lower in the low-risk group (Fig. 10D, E). Our analysis also revealed a positive correlation between risk score and cancer stem cell content (Fig. 10F). In the IMvigor210 cohort, metastatic urothelial carcinoma patients treated with anti-PD-L1 agents were divided into high- and low-risk groups based on the model. The low-risk group had a higher proportion of PD-L1-expressing tumor-infiltrating immune cells and tumor cells (Supplementary Fig. 15). Although overall survival (OS) differences between the groups were not significant (P = 0.362), the low-risk group showed a trend towards better survival (Supplementary Fig. 16). Collectively, these results suggest that low-risk patients may have a better response to immunotherapy and may benefit from various treatments.

Fig. 10

The potential role of CRGRM in immunotherapy. A Comparison of immune checkpoint gene expression between high-risk and low-risk HCC patients. B Comparison of 13 immune functions between high-risk and low-risk HCC patients. C Violin plots revealing differences in IPS between risk groups. Boxplots were used to elucidate the differential expression of HLA genes (D) and M6A genes (E) in high-risk group and low-risk group, respectively. F Scatter plot of the association between cancer stem cells and risk scores

Tumor immune dysfunction and exclusion (TIDE) and immune subtypes analysis

We calculated the TIDE score for each patient using the TIDE algorithm and visualized the differences between groups with a violin plot (Fig. 11A). The results indicated that low-risk patients were more likely to benefit from immunotherapy, with a higher response rate compared to high-risk patients (Fig. 11B, C). Furthermore, immune subtype distribution differed significantly between the groups, with a higher proportion of inflammatory responses in the low-risk group (Fig. 11D). Boxplots further confirmed the significant association between risk scores and immune subtypes (Fig. 11E).

Fig. 11

TIDE and immunophenotype analysis. A Violin plot shows the difference in TIDE score between the different groups. B Bars presents the proportional distribution of immunotherapy responses by different risk groups. C Comparison of risk scores for different immunotherapy response states. D Bars indicates the proportional distribution of immune subtypes by different risk groups. E Comparison of risk scores of different immune subtypes in HCC patients

Association between risk scores and gene mutations

The high-risk group had a higher mutation probability than the low-risk group, and TP53 was the gene with the highest mutation frequency (Fig. 12A, B). In addition, survival analysis showed that patients with high TMB tended to have a worse prognosis, and combined survival analysis also showed that TMB combined with risk score could more accurately predict the prognosis of patients (Fig. 12C).

Fig. 12

To assess the tumor mutation burden (TMB) in different groups. Waterfall plots reflect the frequency of gene mutations in HCC patients in the high-risk (A) and low-risk (B) groups. (C) Stratified survival analysis of TMB combined with risk groups

Functional enrichment analysis

We identified differentially expressed genes (DEGs) between model subgroups using the “limma” R package. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses revealed that DEGs were enriched in “condensed chromosome,” “chromosomal region,” and “microtubule motor activity” (Fig. 13A), and were significantly associated with the “p53 signaling pathway” and “ECM-receptor interaction” (Fig. 13B). Additionally, GSVA analysis identified highly enriched pathways in the model subgroups (Fig. 13C, D). To further validate these findings, we performed Gene Set Enrichment Analysis (GSEA), which complemented the GO and KEGG results (Supplementary Fig. 17).

Fig. 13

Functional enrichment analysis of DEGs. A Bar plot of GO enrichment analysis. B Bar plot of KEGG pathways enrichment analysis. C GSVA analysis based on KEEG gene function sets. D GSVA analysis based on GO gene function sets

Evaluation of drug susceptibility based on risk groups

We used the “pRRophetic” R package to calculate IC50 values for 138 antitumor drugs and visualized the distribution differences between risk groups using box plots. The results revealed lower IC50 values for Imatinib, Erlotinib, Dasatinib, Vorinostat, and Gefitinib in the low-risk group, indicating higher drug sensitivity. Conversely, the high-risk group showed greater sensitivity to Gemcitabine, Sorafenib, Etoposide, and Doxorubicin (Fig. 14). These findings provide valuable insights for clinical chemotherapy drug selection and offer new perspectives on personalized treatment.

Fig. 14

To explore the correlation between risk scores and drug sensitivity

BAMBI was identified as a prognostic marker for HCC

Previous studies suggest that BAMBI promotes tumor growth and protects cancer cells, with its expression elevated in cancer tissues [11]. The role of BAMBI in different cancers remains unclear, and data on its involvement in hepatocellular carcinoma (HCC) are limited. Our analysis indicates that high BAMBI expression in HCC is associated with poor prognosis.

Validation of expression levels of model genes

Tissue samples from 36 HCC patients at the Affiliated Hospital of Xuzhou Medical University were used to validate model gene expression. Real-time PCR revealed higher expression of GMPS, DNAJC6, BAMBI, MP2L2, ASPHD1, and CXCL9 in tumor tissues, while IL7R, EPO, and BBOX1 were more highly expressed in non-tumor tissues (Supplementary Fig. 18). IHC confirmed increased BAMBI expression in HCC tissues (Fig. 15A), and Western blotting corroborated these findings (Fig. 15B), supporting our bioinformatics results. Further real-time PCR testing showed efficient BAMBI silencing using two shRNAs (Fig. 15C), and knockdown of BAMBI also reduced NRP1 and VTCN1 expression (Fig. 15D, E).

Fig. 15

Verify the mRNA expression levels of the model genes in the tissues. A BAMBI representative IHC stained images in HCC tissue and adjacent tissues (magnification: left, 20×; right, 40×). B Comparison of the BAMBI protein expressions in adjacent tissues and HCC tissues. C The bars showed that BAMBI could be efficiently silenced by two independent shRNAs. D Bars reflect the effect of inhibition or non-inhibition of BAMBI expression on VTCN1 protein expression levels. E Bars reflect the effect of inhibition or non-inhibition of BAMBI expression on NRP1 protein expression levels

Knockdown of BAMBI reduced the activity, proliferation and migration of HCC cells

This study used sh-BAMBI#1 and sh-BAMBI#2 to knock down BAMBI expression and construct stable HCC cell lines for functional validation. Western blotting confirmed efficient BAMBI silencing by both shRNAs (Fig. 16A). Cell viability and Edu assays revealed that BAMBI knockdown suppressed HepG2 and MHCC-97H cell proliferation (p < 0.05, Fig. 16B–E). Wound healing assays showed delayed closure in BAMBI-depleted cells (Fig. 16F, G). Transwell assays demonstrated reduced cell invasion in both cell lines (Fig. 16H, I), and clonogenic assays confirmed reduced colony formation (Fig. 16J, K). Moreover, BAMBI knockdown led to decreased expression of immune checkpoints VTCN1 and NRP1 at the protein level. In conclusion, BAMBI promotes HCC cell proliferation and migration, suggesting it as a potential therapeutic target for HCC.

Fig. 16

BAMBI promotes the proliferation, invasion and migration of HCC cells. A Western blot analysis showed that in HepG2 and MHCC-97H cells BAMBI can be stably knocked down by both shRNAs, and the expression of NRP1 and VTCN1 was positively correlated with BAMBI expression. B, C CCK8 assay showed that the viability of HepG2 and MHCC-97H cells was significantly inhibited after BAMBI knockdown compared with the control group. D, E Edu assay showed that the proliferation rate of HepG2 and MHCC-97H cells was significantly inhibited after BAMBI knockdown compared with the control group. F, G Wound healing assay showed that the rate of wound healing was significantly inhibited in HepG2 and MHCC-97H cells after BAMBI silencing. H, I Transwell assays indicated that BAMBI knockdown was inhibited Cell invasion ability of HepG2 and MHCC-97H cells. J, K The colony formation ability of the control group and BAMBI knockdown group was compared

Discussion

The incidence of liver cancer has continued to increase over the past few decades, and liver-cancer related mortality has increased by more than 2% per year since 2007 [12]. HCC is associated with poor survival outcomes and has traditionally had few available treatment options [13]. Biomarkers are crucial for diagnosing, prognosticating, and predicting treatment response, enhancing patient stratification and clinical outcomes. This study aims to identify effective biomarkers based on cuproptosis to advance HCC treatment and immunotherapy.

As a newly discovered mode of cell death, the induction of cuproptosis is different from traditional programmed cell death [14], and its mechanism has not been fully confirmed. Studies have demonstrated that cuproptosis is significantly linked to prognosis and immune landscape changes in bladder cancer [15], breast cancer [16] and head and neck cancer [17]. Fang et al. found that higher serum copper levels and an elevated copper/zinc ratio may correlate with poorer HCC survival [18]. However, no study has explored the role of cuproptosis in targeted drug selection and tumor stratification in HCC.

Studies have shown that BAMBI plays a key role in osteosarcoma pathogenesis and progression [19]. Our research also suggests that BAMBI affects HCC proliferation and metastasis, with its knockdown delaying HCC progression. Additionally, BAMBI correlates with immunotherapy targets such as NRP1 and VTCN1, making it a potential therapeutic target to improve HCC prognosis. Yang et al. found that DNAJC6 promotes HCC progression by inducing epithelial-mesenchymal transition [20]. GMPS expression has been linked to Gleason score in prostate cancer [21], while BBOX1 depletion inhibits TNBC cell growth without affecting normal breast cells, highlighting its potential as a therapeutic target for triple-negative breast cancer [22].

Tumor mutation burden (TMB) serves as a biomarker for predicting clinical response to immune checkpoint inhibitors (ICI) and reflects their underlying mechanism [23]. High TMB is associated with poorer overall survival (OS). We also found that combining TMB with survival analysis better differentiates patient prognosis. The waterfall plot revealed a higher incidence of TP53 mutations in the high-risk group and CTNNB1 mutations in the low-risk group. TP53 mutations, common in human cancers, contribute to tumorigenesis through various mechanisms [24].

This study found that immune checkpoint gene expression was higher in the low-risk group, which also exhibited greater immune cell infiltration. Additionally, the low-risk group had higher IPS, suggesting a better response to immune checkpoint inhibitors (ICI). TIDE analysis and drug sensitivity evaluation confirmed that our model could predict immunotherapy effectiveness and guide clinical drug selection. Gene function analysis revealed that the high-risk group was enriched in metabolism-related pathways, while the low-risk group showed greater enrichment in immune and drug metabolism pathways. Further analysis of immune subtypes [25] revealed that the high-risk group had higher proportions of C1 (wound healing) and C2 (IFN-γ-dominant) subtypes, while the low-risk group had more C3 (inflammation) and C4 (lymphocyte depletion) subtypes. C4 is enriched in specific hepatocellular carcinoma (LIHC) subtypes with a prominent macrophage signature, Th1 suppression, and higher M2 response [26]. Using the IMvigor210 cohort, we showed that the low-risk group had more tumor immune cell infiltration and PD-L1 expression, suggesting greater potential benefit from immunotherapy. This model can guide immunotherapy selection and predict HCC prognosis. In vitro, BAMBI expression was significantly higher in HCC tissues than in normal liver tissues. BAMBI knockdown significantly affected HCC cell proliferation and migration. We also confirmed the positive association of BAMBI with VTCN1 (B7H4) and NRP1, offering new insights for immune-targeted therapy.

Compared to the study by Peng et al. [27], we screened more potential cuproptosis-related genes using gene functional databases, significantly enhancing the predictive value of our model, which was validated across multiple datasets. In addition, we identified a novel prognostic marker for hepatocellular carcinoma (HCC), BAMBI, and thoroughly investigated its upstream and downstream mechanisms in tumor progression. The research by Wang et al. provided valuable insights into multi-sample data standardization, and inspired us to analyze the immunotherapeutic efficacy of our model from multiple perspectives, thereby increasing the clinical application potential of biomarkers [28, 29]. Overall, the strength of this study lies in its integration of multi-omics data, which, for the first time, reveals the association between cuproptosis and the HCC immune microenvironment, offering new strategies for precision treatment of liver cancer. Our cuproptosis-based risk score model, validated using TCGA, GEO, and other datasets, demonstrates high predictive accuracy. Furthermore, BAMBI, identified as an emerging immune therapy target, was experimentally shown to significantly inhibit HCC cell proliferation and migration, further supporting its potential in immunotherapy.

The cuproptosis-based risk scoring model developed in this study demonstrated strong performance across multi-omics datasets. However, its clinical utility and correlation with immunotherapy efficacy require further validation with larger clinical cohorts. Additionally, the specific mechanisms and signaling pathways of BAMBI in hepatocellular carcinoma remain to be elucidated through in vivo experiments and molecular studies. Furthermore, the limited sample sizes in the IHC and WB experiments may affect the statistical reliability of the results.

In summary, we developed and validated a cuproptosis-related risk model using multiple datasets, demonstrating strong performance in predicting prognosis and immunotherapy response. The study also confirmed the model's gene expression in HCC and identified BAMBI as a potential therapeutic target.