Introduction

Primary liver cancer is the third leading cause of cancer-related deaths worldwide in 2022. Hepatocellular carcinoma (HCC) is the main histological subtype of primary liver cancer, accounting for 75-85% of cases1. The main treatment protocols of HCC include surgical resection, liver transplantation, chemotherapy, radiotherapy, targeted therapy and immunotherapy, among which surgical resection is the most effective treatment2,3,4. However, only 5–10% of patients can receive hepatic resection treatment5. Because of the insidious onset of HCC, approximately 80% of patients are diagnosed at intermediate or advanced stages, at which time the opportunity for optimal surgical intervention is missed6. For patients with advanced HCC, therapeutic options are further constrained by tumor progression and drug resistance7. The emergence of immunotherapy has brought new hope to HCC patients. Immune checkpoint inhibitors (ICIs) therapy is one of the most mature immunotherapy technique at present, and its application has significantly enhanced the survival rate of HCC patients8. However, ICIs therapy is applicable for a small subset of HCC patients due to complex pathogenesis, immunosuppressive tumor microenvironment (TME), and immune resistance9. Furthermore, HCC is characterized by high intratumoural heterogeneity, further complicating its management10. Consequently, the prognosis for HCC remains poor, with a five-year survival rate of only 18%11.

NAD(P)H dehydrogenase quinone 1 (NQO1), also known as DT diaphorase, is a flavoenzyme and located on chromosome 16q22 in humans12,13. NQO1 plays a pivotal role in reprogramming the TME by promoting the secretion of cytokines such as IL6, CXCL10 and GM-CSF, thereby contributing to the formation of an immunosuppressive microenvironment14. When NQO1 was knocked out in mice, CD8+ T cell infiltration increases significantly, and the efficacy of PD-1 monoclonal antibody therapy is markedly enhanced15. As a two-electron reductase, NQO1 catalyzes the reductions of quinones to the hydroquinones, using NADH or NADPH as the electron donor16. When the resulting hydroquinone is unstable, it can undergo auto-oxidation, generating reactive oxygen species (ROS) or producing compounds capable of DNA alkylation. This property lays the foundation for the development of NQO1-targeted antitumor drugs17,18. For instance, β-lapachone and deoxybenzoquinone directly induce DNA damage and kill cancer cells by generating excessive ROS through redox cycling. In contrast, mitomycin C promotes cancer cell apoptosis by alkylating DNA19. Moreover, NQO1 is overexpressed in a variety of cancer tissues (e.g., pancreatic, breast, and lung cancers) compared with normal tissues, and its overexpression is correlated with poor prognosis of patients20,21,22. This makes NQO1 a potential key molecular target in cancer development, treatment, and prognosis. Similarly, NQO1 is upregulated in HCC and has been closely linked to tumor size, progression, and poor prognosis23,24,25. However, the current understanding of the role of NQO1 in HCC remains limited, particularly regarding its involvement in the tumor microenvironment, prognosis, and immunotherapy.

In recent years, rapid advances in sequencing technology have provided strong support for identifying key molecular features driving the onset or progression of HCC. This study aimed to investigate the role of NQO1 in HCC from a multi-omics perspective. Firstly, we examined the expression of NQO1 at the tissue and cellular levels using multiplex immunofluorescence and single-cell transcriptome analysis, and investigated the role of NQO1 in TME by spatial quantification analysis. Then, core NQO1-related genes (NRGs) associated with HCC were identified, and an NRGs-related risk-scores in hepatocellular carcinoma (NRSHC) was constructed based on these genes to predict survival outcome and response to immunotherapy of HCC patients. To further enhance predictive accuracy, we integrated the risk model with other clinical factors to create a nomogram and developed an online prediction tool. We hope that our work can provide valuable insights for clinical decision-making in HCC.

Methods

Data acquisition and processing

Six tumor samples from GSE149614 dataset were downloaded from the Gene Expression Omnibus(https://www.ncbi.nlm.nih.gov/geo/) for single-cell transcriptome analysis. Cells with gene expression numbers fewer than 200 or more than 2500 or with a proportion of mitochondria exceeding 5% were excluded. Ultimately, a total of 14,912 cells were included in subsequent analysis. The gene expression, somatic mutation and clinical data from the GDC TCGA-LIHC datasets were obtained visa the UCSC Xena platform (https://xenabrowser.net/)26,27. A total of 337 tumor samples with complete data were included. To develop prognostic model, the TCGA-LIHC dataset was divided randomly into training (N = 238) and validation sets (N = 99) in a 7:3 ratio using the R package “caret“28. Besides, the gene expression data from GSE5423629 dataset was downloaded to validate the expression of NRGs in HCC, while the gene expression and clinical data from GSE1452030 dataset were used for external validation of the prognostic model. The 4573 NRGs were obtained from the GeneCards database (https://www.genecards.org/)31. All data were normalized and standardized to ensure reliable and comparable results. Detailed information about the cohorts used for modeling is summarized in Table 1.

Table 1 Descriptive characteristics of study cohorts for modeling.
Full size table

Dimensionality reduction, clustering, and annotation of cell types

Principal component analysis was performed to reduce the dimensionality of data. Cell clustering was conducted using the FindNeighbor and FindCluster functions. Subsequently, the identified cell clusters were visualized using uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) algorithm. Finally, the FindAllMarker function was used to identify differentially expressed genes (DEGs) between each cluster and others. The top 50 DEGs in each cluster were manually annotated using the CellMarker 2.0 database (http://117.50.127.228/CellMarker/)32. The above analyses primarily performed using the “Seurat” package33.

Cell trajectory and cell-cell communication analysis

After cell clusters were identified and annotated, the expression of NQO1 in each cell type was visualized using the FeaturePlot and VlnPlot functions in the “Seurat” package. Next, the cell type with the highest NQO1 expression was selected for cell trajectory analysis using the “monocle2” package34. DEGs were first identified using the differentialGeneTest function, and dimensionality reduction was performed with the DDRTree algorithm. And then cell trajectories were constructed based on the expression trends of the DEGs. Finally, the expression of NQO1 along the cell developmental trajectory was visualized using the plot_genes_violin and plot_genes_jitter functions.

The “Cellchat” package35 was used to perform intercellular communication analysis. The “Secreted Signaling” subset of CellChatDB.human was selected as the ligand-receptor interaction database. Overexpressed ligand-receptor pairs were identified, and the processed gene expression data were mapped onto the protein-protein interaction network. The computeCommunProb function was then used to calculate the communication probability. The communication information between cell types was visualized in circle graphs.

Immunohistochemistry and Immunofluorescence analysis

A human HCC tissue microarray (TMA) that consisted of 90 tumor and 90 paired non-tumor samples from 90 patients with HCC, with a total of 180 tissue spots, was used to detect proteins expression. Specifically, two tissue spots (1.5 mm diameter) per patient were used. All specimens were fixed in formalin and embedded in paraffin. The clinical information of the patients was shown in Table 2. The target genes were fluorescently labeled using the AlphaXTSA 7-color immunohistochemical staining kit (Alpha X Bio, AXT37100041). For immunostaining on tissues, slides were progressively dewaxed with dewaxing reagent (LEAGENE, DZ2011), rehydrated with different concentrations of ethanol, and then boiled for 25 min in EDTA solution (ZSGB-BIO, ZLI-9079) for antigen retrieval. After cooling, slides were blocked with a 5%BSA blocking solution and incubated for 15 min. The blocking solution was then removed, and the slides were incubated with primary antibody against AFP (HCC marker, Thermofisher, TA809081, 1:150) for one hour at 37 °C. After washing the sections three times with PBS, the slides were incubated with the secondary antibody (NQO1, Proteintech, 67240-1-Ig, 1:200) for one hour at 37 °C. To detect M1 and M2 macrophages, anti-CD68 antibody (General macrophage marker, ZSGB-BIO, ZM0060, 1:300), anti-CD80 antibody (M1 marker, Proteintech, 66406-1-Ig, 1:100), and anti-CD163 antibody (M2 marker, ZSGB-BIO, ZM0428, 1:200) were used, following the same procedure as previously described. Cell nuclei were counterstained with DAPI. Finally, Antifade Mounting Medium (Phygene, PH0429) was applied to the slides to reduce fluorescence quenching.

Scanning and imaging were performed using ZEISS AXIOSCAN 7 multispectral imaging system and ZEN 3.3 software. Cell quantification and spatial distance analysis were conducted using StrataQuest 7.1.1.145 software. After excluding samples with missing data, 86 samples were included in the cellular quantitative analysis. Among these, 74 pairs of samples were used for paired t-test analysis, and 86 tumor samples were used for differential analysis and survival analysis. The clinical information of the samples is provided in Table 2. For the survival analysis, the surv_cutpoint function from the “survival” package was utilized to determine the optimal cut-off point for the variable. Based on this cut-off, the variable was categorized into high-risk and low-risk groups to evaluate the impact of different variable levels on the survival probability of patients.

All experiments were approved by the Medical Ethics Committee of the First Affiliated Hospital of the University of South China (Approval No.2024LL0320001) and were conducted in accordance with Declaration of Helsinki. Informed consent was obtained from all subjects and/or their legal guardians.

Table 2 Clinical characteristics of the TMA cohort.
Full size table

Development, validation and comparison of the NRSHC

Using the R package “limma“36, differentially expressed NRGs were identified with a threshold of false discovery rate q < 0.05 and |fold change|>1 between HCC and normal tissue in the TCGA-LIHC dataset. The differential expression of NRGs was visualized using the “pheatmap” and “ggplot2” packages37. Weighted gene correlation network analysis (WGCNA) was performed with the “WGCNA” package38 to identify key modules and genes associated with HCC. Univariate Cox regression analysis was conducted to identify NRGs with significant prognostic values (P < 0.05). Subsequently, the common NRGs identified through the three feature selection methods mentioned above were subjected to functional enrichment analysis for the Kyoto encyclopedia of genes and genomes (KEGG) and Gene Ontology (GO) using the R package “clusterProfiler“39, to determine their functions and roles in HCC. Next, the least absolute shrinkage and selection operator (LASSO) Cox regression was performed to screen critical NRGs, which were then used for multivariate Cox regression analysis to construct a prognostic risk assessment model for HCC. The analyses were carried out using the R package “glmnet“40. The resulting model formula was as follows:

$$Risk~score = \mathop \sum \limits_{{i = 1}}^{n} \left( {Expression~of \times coefficient~of~gene} \right)$$

Expression and coefficient of gene represent respectively mRNA expression level and regression coefficient values of NRGs. The risk model was first constructed in TCGA-train cohort, and then tested in TCGA-test and GSE14520 cohort to confirm the accuracy. Finally, we compared the model constructed in this study with existing HCC models from three dimensions: predictive accuracy, clinical benefit, and survival risk, to comprehensively evaluate the model’s reliability. Time-dependent Receiver Operating Characteristic (ROC) curve analysis and the C-index were used to assess the model’s predictive accuracy, implemented through the R packages “timeROC“41 and “survcomp”, respectively. Clinical benefit was evaluated using decision curve analysis (DCA). Additionally, restricted mean survival time and median survival time were employed to compare the predictive performance for survival risk. To analyze the characteristics of the model, HCC samples from each cohort were classified into high-risk and low-risk groups based on the median risk-scores. The comparison analysis between the two risk groups used the data of TCGA-all cohort.

Construction of nomogram and its web application tool

In order to improve the predictive power of the risk model, a prediction nomogram was constructed based on the risk model and clinical risk factors of HCC using the R package “regplot“42. Calibration curves were generated with the the R package “rms” to evaluate the nomogram’s accuracy. Time -dependent ROC curves were used to estimate the predictive performance of the risk model and nomogram, and DCA was conducted using the R package “ggDCA“43 to evaluate their clinical benefit. To enhance the practical application of nomograms, a web-based prediction tool―NRSHC TOOL (http://wys.helyly.top/cox-t/cox.html) was developed. If the relevant information of an HCC patient is entered into the NRSHC TOOL, the predictive results of his drug sensitivity and survival probability will be displayed.

Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

The “c2.cp.kegg.v7.4.symbols” file was downloaded from Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb)44,45. The GSEA analysis was performed with the “limma” and “clusterProfiler” package to analyze the potential signaling pathways of two risk subgroups46. The five most significant pathways of KEGG in each subgroup were visualized using enrichment plots. Additionally, GSVA was conducted to assess the functional differences between the high- and low-risk groups using the R package “GSVA“47. A heatmap was generated to display the 50 pathways with the most significant differences between the two groups.

Immune cell infiltration analysis

First, the ESTIMATE algorithm48 was applied to calculate stromal scores, immune scores, ESTIMATE scores and tumor purity of patients. The differences in these metrics between the high-risk and low-risk groups were then compared. Next, the infiltration levels of 29 immunocytes in each sample were quantified with the single-sample gene set enrichment analysis. Finally, the proportions of 22 immune cell types were quantified using the CIBERSORT algorithm49, and their differences in the two groups were compared using bar plots. The correlation between risk-scores and immune cell types was assessed through spearman correlation analysis.

Prediction of treatment response

We identified 25 immune checkpoint genes from previous studies and analyzed their differences in mRNA expression levels between the two risk groups. Tumor mutational burden (TMB) analysis was performed using the R package “maftools“50. The twenty genes with the most mutations frequencies in HCC samples were visualized using oncoplots for each risk group. Immunophenoscore (IPS) data from 371 patients with HCC were downloaded from The Cancer Immunome Atlas database (https://www.tcia.at/)51 and compared between different risk groups. These analyses aimed to predict the response to immunotherapy in HCC patients. Besides, the IC50 values of eight chemotherapy drugs (cisplatin, docetaxel, doxorubicin, gemcitabine, paclitaxel, vinorelbine, mitomycin C, bleomycin) of HCC in different risk groups were evaluated using the R package “pRRophetic“52 to predict the efficacy of chemotherapy in different subgroups.

Characteristic analysis of critical NRGs

We performed a detailed feature analysis of the three NRGs ultimately selected. A forest plot illustrating the prognostic features and a lollipop chart displaying the regression coefficients were generated using the SRPLOT website (https://www.bioinformatics.com.cn/en)53. The mutation and variation status of each gene were analyzed via the cBioPortal platform (https://www.cbioportal.org/)54. Moreover, the promoter methylation and protein expression level of each gene were evaluated through UALCAN website (https://ualcan.path.uab.edu/)55.

Statistical analysis

Student’s t test and Wilcoxon rank-sum test were used to compare differences of variables. Univariate and multivariate Cox regression analyses were conducted to evaluate the impact of variables on the prognosis of patients. The log-rank test was used to assess survival differences among patients. With the exception of analyses performed using free online tools, all statistical analyses were conducted using R software (version 4.1.3 and 4.3.1). All statistical P values were two-sided and P<0.05 was considered statistically significant.

Results

Landscape of NQO1 expression in HCC

The workflow of this study is illustrated in Fig. 1. We first investigated the expression profile of NQO1 in HCC. Multiple immunofluorescence analysis revealed that the fluorescence expression intensity of NQO1 and AFP in HCC tissues was significantly higher than that in normal tissues (Fig. 2A). Notably, there was a statistically significant difference in the proportion of NQO1+AFP+ cells between tumor and normal tissues, with a greater distribution of NQO1+AFP+ cells observed in HCC tissues (Fig. 2B). To examine the expression of NQO1 at the cellular level, single-cell transcriptome analysis was implemented. After the preprocessing, we obtained 14,912 cells, which were classified into 25 clusters (Fig. 2C). Subsequently, DEGs were identified between clusters to annotate cell types (Fig. 2D), resulting in the identification of six major cell types (Fig. 2E). Among these, NQO1 was found to be highly expressed in HCC cells (Fig. 2F–G). Consequently, we further explored the differentiation trajectory of the HCC cells and the expression profile of NQO1 within these cells. The results showed that the HCC cells transitioned from rapid growth to gradual apoptosis (Fig. 2H). Furthermore, while the overall expression of NQO1 was low, its high expression period was mainly concentrated in the middle and late stages of HCC cell life (Fig. 2I). In summary, NQO1 was found to accumulate in HCC cells within HCC tissues.

Fig. 1
figure 1

The workflow of this study. (A) NQO1 expression in HCC tissues and HCC cells. (B) Role of NQO1 in HCC microenvironment. (C) Feature dimensionality-reduction of NRGs. (D) Construction and application of a prognostic risk model for HCC based on NRGs. (E) Analysis of characteristics of the high and low risk groups.

Full size image
Fig. 2
figure 2

Landscape of NQO1 in HCC. (A) Difference in fluorescence expression of AFP, NQO1, and NQO1+AFP+ between HCC samples and paired normal samples. (B) Statistical differences in expression of NQO1+AFP+ between HCC samples and paired normal samples. (C) The UMAP and t-SNE plots show 25 cell clusters formed by clustering. (D) The heatmap exhibits the top 5 marker genes in each cell cluster. (E) The UMAP plot displays major cell types from 6 HCC samples. (F–G) The expression level of NQO1 in six cell types. (H) The trajectory maps of HCC cell colored by state and pseudotime. (I) The expression changes of NQO1 during the growth of HCC cells.

Full size image

NQO1+HCC cells communicate frequently with M2 macrophages

We next investigated the role of NQO1 in TME. Among the communication networks between the six cell types, HCC cells demonstrated the highest number and strength of interactions with tumor-associated macrophages (TAMs) (Fig. 3A–B). Given on our previous findings that NQO1 was highly expressed in HCC cells, HCC cells were divided into NQO1+HCC cells and NQO1−HCC cells to examine the differences in their interactions with TAMs. The results revealed that NQO1+HCC cells had a significantly stronger interaction with TAMs than NQO1−HCC cells (Fig. 3C–D). Furthermore, by subdividing TAMs into M1 and M2 macrophages, we observed that NQO1+HCC cells communicate frequently with M2 macrophages (Fig. 3E–F).

Fig. 3
figure 3

Cell-cell communication analysis. (A–B) Circle plot shows the number and strength of interactions between cell types. (C–D) The number and strength of interactions between cancer cells stratified by NQO1 and other cell types. (E–F) The number and strength of interactions between cancer cells and macrophage.

Full size image

NQO1+HCC cells May promote M2 polarization of TAMs

The intercellular communication analysis revealed that NQO1+ HCC cells and M2 macrophages frequently interacted within the TME. To further validate this finding, multiplex immunohistochemistry and spatial quantification analyses were conducted. Fluorescence imaging of NQO1, AFP, CD68, CD80, and CD163 is presented in Fig. 4A. The distances between NQO1+AFP+ cells and CD68+CD80+ cells (M1 macrophages), as well as CD68+CD163+ cells (M2 macrophages) in tumor samples were calculated and visualized (Fig. 4B). The results demonstrated that the number of M2 macrophages was significantly higher than that of M1 macrophages within 0–50 μm and 50–100 μm ranges from NQO1+AFP+ cells (Fig. 4C–D). These findings suggest that NQO1+ HCC cells may play a crucial role in the M2 polarization of TAMs. We also analyzed the correlation between the prognosis of HCC patients and the number of NQO1+AFP+ cells and the number of M2 macrophages. The results showed that the number of both was inversely associated with the survival probability of HCC patients, indicating that higher counts of NQO1+AFP+ cells and M2 macrophages correspond to poorer patient prognosis. (Fig. 4E–F).

Fig. 4
figure 4

Multiple immunofluorescence and spatial quantitative analysis. (A) Fluorescence expression of all markers in representative tumor samples. (B) Schematic plot of the distribution of M1 and M2 near NQO1+AFP+ cells in representative region. Gray represents DAPI, purple represents NQO1+AFP+ cells, yellow represents CD68+CD80+ cells (M1), and green represents CD68+CD163+ cells (M2). In addition to this, cyan indicates the connection between NQO1+AFP+ cells and M1, and white indicates the connection between NQO1+AFP+ cells and M2. (C) Counts of M1 and M2 cells within 0-100 μm from NQO1+AFP+ cells in representative region. (D) Differences in the number of M1 and M2 cells within 0–50 μm and 50–100 μm from NQO1+AFP+ cells after logarithmic transformation. (E) Difference in overall survival time between HCC patients with high NQO1+AFP+ cell counts and low NQO1+AFP+ cell counts. (F) Difference in survival of HCC patients with different M2 cell numbers within 0–50 μm distance from NQO1+AFP+ cells.

Full size image

Identification and functional enrichment of the key NRGs

The previous results confirmed that high expression of NQO1 was associated with poor prognosis in HCC. Therefore, we further investigated the prognostic value and clinical application of NRGs in HCC. Differential expression analysis showed 541 NRGs were differentially expressed in HCC tissues and normal tissues in the TCGA-LIHC dataset, and among which 329 genes were upregulated in HCC tissues (Fig. 5A–B). A similar result was observed in GSE54236 dataset (Fig. 5C). Next, we performed WGCNA analysis to further reduce dimensionality. The dendrogram of NRGs and the corresponding color modules are shown in Fig. 5D. The heatmap indicated that the MEblue module had the highest correlation with HCC, showing a strong positive correlation with HCC tissues (R = 0.78, P = 4e − 17) and a negative correlation with normal tissues (Fig. 5E). This finding suggests that the 198 NRGs in the MEblue module were hub genes associated with HCC. Meanwhile, univariate cox regression analysis identified 984 NRGs related with HCC prognosis. By intersecting the 541 DEGs, 198 hub genes, and 984 prognostic genes, we identified 87 key NRGs, which were subsequently used for enrichment analysis (Fig. 5F). KEGG pathway enrichment analysis revealed that the key NRGs were predominantly involved in pathways related to cell growth and metabolism, viral infections, and liver diseases. These pathways included cell cycle, DNA replication, cellular senescence, human T-cell leukemia virus 1 infection, viral carcinogenesis, hepatocellular carcinoma, and hepatitis B (Fig. 5G). Similarly, GO enrichment analysis showed that these genes were enriched in biological processes such as DNA replication and regulation, mitotic cell cycle phase transition, and molecular functions such as catalytic and ATP-dependent activity acting on DNA. They were also localized to chromosomal regions, nuclear chromosomes, and protein-DNA complexes (Fig. 5H).

Fig. 5
figure 5

Identification of the core NRGs and functional enrichment analysis. (A–B) The differential expression of NRGs in TCGA-LIHC. (C) The differential expression of NRGs was tested in GSE54236. (D–E) The dendrogram and heatmap depicts respectively the clustering of genes and module-trait relationships from WGCNA analysis. (F) The overlapping NRGs of the three feature selection methods. (G) The enriched pathway of the core NRGs. (H) The role of the core NRGs in cellular components, molecular functions, and biological processes.

Full size image

Construction and validation of prognostic risk model in HCC—NRSHC

To develop a prognostic risk model for HCC, we first performed LASSO regression analysis, which identified 3 NRGs based on the optimal lambda value (Fig. 6A–B). Univariate Cox regression analysis further demonstrated that these 3 NRGs had significant prognostic value (Fig. 6C). Consequently, a risk model was constructed using multivariate Cox regression analysis in the TCGA-train cohort, and the formula was: The NRGs-related risk-scores in hepatocellular carcinoma (NRSHC) = TUBG1 × 0.423 + CAD × 0.415 + BSG × 0.344 (Fig. 6D). Subsequently, each cohort were stratified into high- and low-risk groups according to the median risk score. Time-independent ROC curve analysis revealed that the area under the curve (AUC) for NRSHC at 1, 3, and 5 years in the TCGA-train cohort was approximately 0.7. Survival analysis indicated that patients in the high-risk group had a significantly lower survival probability in the TCGA-train cohort (Fig. 6E).

To assess the reliability of the model, we conducted internal and external validation. In the TCGA-test cohort, the AUC values for the 1-, 3-, and 5-year survival endpoints were 0.791, 0.618, and 0.630, respectively, and patients in the high-risk group exhibited worse overall survival (Fig. 6F). Consistent results were observed in the TCGA-all cohort (Fig. 6G). External validation was performed using the GSE14520 cohort, which demonstrated AUC values of approximately 0.6 at 1, 3, and 5 years. Similarly, the high-risk group had a significantly worse prognosis (Fig. 6H). In a word, the NRSHC model was constructed successfully, exhibiting high sensitivity, specificity as well as robust predictive capability.

Fig. 6
figure 6

The construction of the NRSHC. (A–B) Dynamic process of screening modeling NRGs through Lasso regression. (C) The forest plot shows the impact of 3 NRGs on the prognosis of HCC. (D) Regression coefficients of 3 NRGs. (E–H) Time-independent ROC curve analysis and survival analysis of the NRSHC model in TGCA-train, TCGA-test, TCGA-all, and GSE14520 cohorts.

Full size image

The nomogram taking the NRSHC into account has good predictive performance

We compared the prognostic value of the NRSHC, gender, age, grade, pathological tumor-node-metastasis (pTNM) Stage using univariate and multivariate Cox regression analyses (Fig. 7A,B). Forest plots revealed that age, the NRSHC and pTNM stage were risk factors for HCC. Additionally, we investigated the association between the NRSHC and clinical characteristics. The results demonstrated that female patients with HCC had higher risk scores than male patients, and individuals with higher tumor grade and advanced pTNM stage also exhibited elevated risk scores (Suppl. 1A–D).

Based on the results from the prognostic analysis, a nomogram integrating age, pTNM stage, and NRSHC was developed to better predict the survival outcome of HCC patients. The nomogram indicated that the higher age, pTNM stage, and risk score corresponded to a higher total score, which in turn was associated with an increased probability of mortality at 1, 3, and 5 years (Fig. 7C). The calibration curve indicated that the OS predicted by the nomogram is close to the actual OS (Fig. 7D). We then compared the predictive accuracy and clinical benefit of the nomogram with other models. Regarding prediction accuracy, the effect of the nomogram was the best, followed by the NRSHC (Fig. 7E–G). And compared with other models, the nomogram provided a modest additional net benefit at 1-, 3-, and 5-years, underscoring its potential clinical utility (Fig. 7H–J).

Fig. 7
figure 7

Construction and evaluation of the nomogram. (A–B) Univariate and multivariate COX regression analyses are used to identify the effect of the NRSHC and clinical factors on the prognosis of HCC. (C) Schematic nomogram of survival outcome prediction for a representative sample. The sample has pTNM stage I, 55 years of age, and is classified as a high-risk group in the NRSHC model, with total score of 74.5, and the predicted 1 -, 3 -, and 5-year mortality rates are 7.47%, 17.2%, and 25.7%, respectively. (D) The calibration curve shows the relationship between the predicted values of the nomogram and the actual values. (E–G) The AUC value of the nomogram and other models at 1, 3, 5 years. (H–J) Clinical net benefit of the nomogram and other models at 1,3,5 years.

Full size image

Comparison of NRSHC and existing HCC prognostic models

We selected three HCC prognostic models based on different gene sets for comparison with the NRSHC model to validate its predictive performance. These three models were: the prognostic signature constructed by Wang et al.56 based on fatty acid metabolism-related genes, the prognostic model developed by Xing et al.57 using pyroptosis-related genes, and the prognostic signature established by Xu et al.58 based on DHRS1-related genes. First, we performed ROC curve analysis for each model. The results showed that the nomogram incorporating NRSHC had the highest AUC values for 1-years, 3-years, and 5-years predictions, followed by the NRSHC model alone (Suppl. 2A–B). The C-index calculation results were consistent with the ROC analysis (Suppl. 2C), indicating that the NRSHC model outperformed the other three models in terms of predictive accuracy. Additionally, DCA analysis demonstrated that both the nomogram and NRSHC had significantly greater clinical benefits than the models proposed by Wang, Xing, and Xu (Suppl. 2D). Finally, we conducted an in-depth comparison of the survival risk prediction capabilities of each model. Whether using restricted mean survival time or median survival time as the evaluation metric, all five signatures, including the nomogram, were identified as risk factors for poor prognosis in HCC patients. Notably, patients with higher model scores tended to have worse survival outcomes (Suppl. 2E–F). In conclusion, these findings demonstrate that NRSHC exhibits superior predictive performance and high reliability, making it a promising tool for clinical application.

Function and immune characteristic of the NRSHC in HCC

To investigate the function characteristics of distinct risk-scores, we performed GSEA and GSVA functional enrichment analysis. The results revealed that the primary signaling pathways enriched in the high-risk group were associated with biosynthesis of RNA, DNA damage repair, and cellular material transport. In contrast, the low-risk group was predominantly enriched in metabolism pathways, such as cytochrome P450 enzyme system, fatty acid and amino acid metabolism, and steroid biosynthesis (Fig. 8A, Suppl1E). Then, the immune scores of 363 HCC samples were calculated and compared. Violin plots exhibited the low-risk group had higher stromal and immune scores and lower tumor purity, which may reflect a more favorable prognosis (Fig. 8B). Furthermore, the low-risk group had a higher proportion of naïve B cell, resting memory CD4 T cells, M2 macrophages, and resting mast cells. Conversely, the high-risk group displayed a markedly higher fraction of memory B cells, activated memory CD4 T cells, regulatory T cells, and M0 macrophages (Fig. 8C). A comprehensive heatmap provided immune profiles of the NRSHC (Fig. 8D). In the end, the correlation between the NRSHC and 22 immune cell types was analyzed. Although the correlation coefficients were relatively low, there was a significant positive correlation between the NRSHC and most T-lymphocyte subtypes (Fig. 8E). These findings suggest that the immune system may be more active in the high-risk group.

Fig. 8
figure 8

The pathways and immune infiltration analysis of the NRSHC. (A) The top five KEEG pathways enriched in high-risk and low-risk groups. (B) Comparison of tumor purity, ESTIMATE score, immune score, stromal score between two risk groups. (C) The fraction of 22 immune cells in two groups. (D) The difference in levels of immune infiltration in different risk groups. (E) The correlation between the NRSHC and infiltration of ten immune cells.

Full size image

The NRSHC is a powerful predictor of immunotherapy response in HCC

The previous findings imply that HCC patients with high risk-scores may exhibit a “hot” tumor phenotype. Based on this, we hypothesized that the patients in high-risk group would benefit more from immunotherapy and investigated to validate this assumption. First, differential expression analysis showed that immune checkpoint genes, with the exception of TMIGD2 and IDO2, were highly expressed in the high-risk group (Fig. 9A). Secondly, TMB was higher in the high-risk group overall. Notably, TP53 was the most frequently mutated gene in the high-risk group, whereas TNN exhibited the highest mutation frequency in the low-risk group. (Fig. 9B–C). Then, IPS was used to assess the immunotherapy effect of different risk subgroups. The results displayed that IPS was higher in the high-risk group than low-risk group when CTLA-4 or PD-1 was inhibited (Fig. 9D). These findings demonstrate that the patients with high risk-scores were more likely to benefit from ICIs therapy. However, predictive analysis of drug sensitivity revealed that patients in the low-risk group exhibited greater sensitivity to chemotherapeutic agents, suggesting that these patients may respond more effectively to chemotherapy (Fig. 9E).

Fig. 9
figure 9

Prediction of treatment response by the NRSHC. (A) Differences in expression levels of immune checkpoint genes between different risk groups. (B–C) TMB analysis of HCC patients in high-risk and low-risk groups. (D) Difference in IPS between the two groups. (E) Differences in IC50 values of chemotherapy drugs between the two risk groups.

Full size image

The comprehensive feature of BSG, CAD, TUBG1

Next, we analyzed the characteristics of the critical NRGs (BSG, CAD, TUBG1) on mRNA expression, genetic variation, methylation level and protein expression. These NRGs were highly expressed in HCC, and exhibited varying degrees of genetic alterations. Among them, BSG had the highest altered frequency, followed by TUBG1 and BSG (Fig. 10A–B). The primary types of genetic variations observed included amplification, mutation, and deep deletion (Fig. 10C). At the level of promoter methylation, BSG exhibited higher expression in HCC samples, whereas CAD and TUBG1 showed higher expression in normal samples (Fig. 10D). In terms of protein expression, there was no difference in the expression level of BSG between tumor and normal samples. However, both CAD and TUBG1 were expressed at higher levels in HCC samples (Fig. 10E).

Fig. 10
figure 10

The comprehensive characteristics of BSG, CAD and TUBG1. (A) Differences in mRNA expression levels of BSG, CAD and TUBG1 between HCC and normal tissues. (B) An oncoplot shows the overview of genetic variants of 3 NRGs. (C) Type of variants of 3 NRGs. (D) Differences in promoter methylation levels of 3 NRGs between HCC and normal tissues. (E) Differences in protein expression levels of 3 NRGs between HCC and normal tissues. Meaning of statistical symbols: ***P < 0.001, ns non-statistical.

Full size image

The design of online predictive sites―NRSHC TOOL

To enhance the clinical utility of this study, we developed a web-based prediction tool named “NRSHC TOOL” to predict the immunotherapy response and survival outcomes of HCC patients (Fig. 11). NRSHC TOOL is a decision support system that integrates gene expression profiles, age, and pTNM stage, featuring three core functions. First, it performs risk stratification based on the expression levels of three genes (TUBG1, CAD, and BSG) in tumor samples. Second, it provides quantitative prognostic evaluation by combining risk level, age and pTNM stage. Finally, it predicts the patient’s potential sensitivity to ICIs therapy.

The tool features a simple and intuitive interface, ensuring ease of use for clinicians, patients, and their families. The operation process is as follows: Users begin by entering the expression levels of TUBG1, CAD, and BSG from the HCC patient into the “Genes” field to calculate the risk scores. The system will then classify patients based on the median risk scores, with scores greater than 0.904179 categorized as high-risk, while those below are classified as low-risk. High-risk patients are predicted to have a poorer prognosis but may exhibit greater sensitivity to immunotherapy. After confirming the risk group, the user input the patient’s age and pTNM stage to obtain their 1-year, 3-year, and 5-year survival probabilities. NRSHC TOOL provides a scientific basis for personalized treatment of HCC, assisting clinicians in optimizing therapeutic strategies and enhancing treatment outcomes, and offers valuable disease management guidance to patients and their families.

Fig. 11
figure 11

Webpage interface of NRSHC TOOL. If a user enters the expression levels of TUBG1, CAD, and BSG, age, and pTNM stage of a HCC patient, then he will get the results of whether the patient is sensitive to pd1, pd-L1, and CTLA-4 and his survival probability of 1, 3, and 5 years.

Full size image

Discussion

Hepatocellular carcinoma is a common malignant tumor of the liver. Despite the availability of multiple treatment modalities, including ablation, surgery, liver transplantation, and chemotherapy, the survival rate of HCC patients remains low due to its insidious onset, molecular heterogeneity, and resistance to treatment10,11. Investigating the key molecular features driving the initiation and progression of HCC may help improve this situation. Studies have shown that NQO1 is overexpressed in HCC tissues and may promote malignant progression and poor prognosis by remodeling the tumor microenvironment. Therefore, NQO1 holds promise as a novel immunotherapeutic target for HCC.

Although multiple studies have confirmed the overexpression of NQO1 in HCC, these studies primarily employed immunohistochemistry or real-time quantitative PCR techniques, which are limited to detecting NQO1 expression at the tissue or RNA level59,60. In this study, we combined multiplex immunohistochemistry, single-cell transcriptomics, and spatial structural quantification to analyze NQO1 expression from multiple dimensions, including tissue, cellular, and spatial levels. Our findings revealed that NQO1 is significantly enriched in HCC cells within tumor tissues, and NQO1+ HCC cells are associated with poorer prognosis. These results suggest that NQO1 may serve as a potential therapeutic target and prognostic biomarker for HCC. However, the specific functions of NQO1 in HCC require further investigation to fully elucidate its critical role in tumor initiation and progression.

TME is a dynamic system composed of many cellular and non-cellular components, and its complexity and diversity are important factors influencing the efficacy of immunotherapy. In the TME interaction network analysis of this study, frequent interactions between HCC cells and TAMs attracted our attention. TAMs are the most abundant infiltrating immune cells in the tumor stroma and are primarily classified into classically activated M1 macrophages and alternatively activated M2 macrophages, playing roles throughout the stages of HCC initiation, progression, and treatment61. To date, research on the correlation between NQO1 and macrophages remains unexplored. However, evidence suggests a subtle relationship between liver cancer cells and TAMs. Using multiplex immunofluorescence staining, we analyzed the spatial distribution of TAMs around NQO1+ HCC cells within the HCC microenvironment. The results showed that, compared to M1 macrophages, M2 macrophages were more abundantly clustered near NQO1+ HCC cells, and their abundance was inversely correlated with the survival time of patients. In the early stages of HCC, TAMs exhibit an inflammatory M1 classical activation phenotype, triggering immune responses that inhibit tumor growth. However, as the tumor progresses and expands, signals such as hypoxia, inflammation, and metabolism within the TME gradually induce the polarization of TAMs into M2 macrophages, which subsequently promote tumor angiogenesis and cancer cell proliferation62,63,64. Although direct evidence linking NQO1 to TAM polarization is currently lacking, its upstream regulator, nuclear factor erythroid 2-related factor 2 (Nrf2), has been shown to drive the polarization of TAMs into the M2 phenotype, conferring TAMs immunosuppressive properties that facilitate tumor growth65,66. Given that NQO1 is upregulated by Nrf2, it is highly likely to play a role in M2 polarization as well. Therefore, we hypothesize that NQO1 may contribute to driving TAMs toward M2 polarization, making it a potential immunotherapeutic target. However, further experimental validation and mechanistic studies are needed to confirm this role.

Machine learning has been widely applied in clinical predictive models, offering significant advantages over traditional methods by overcoming nonlinear fitting challenges and efficiently handling large-scale datasets with enhanced predictive power. To the best of our knowledge, this study is the first to employ a machine learning approach to develop and evaluate a clinical predictive model for HCC associated with NRGs. To reduce redundant information and improve model performance, we applied multiple rounds of feature dimensionality reduction to NRGs, ultimately identifying TUBG1, CAD, and BSG as the most critical NRGs in HCC. These three genes are closely associated with HCC and may promote tumor progression67,68,69. Based on this, an NRGs-related prognostic risk model for HCC, NRSHC, was developed. Assessed by the ROC curve, the model demonstrated good predictive accuracy in the training set and internal validation set but performed less effectively in the external validation cohorts. We believe this discrepancy may be attributed to dataset heterogeneity, including differences in data processing methods, technical platforms, the source and feature of sample. The GSE14520 dataset utilized two distinct sequencing platforms, GPL571 and GPL3921, which may have introduced systematic bias. In addition to the GSE14520 cohort, the ICGC cohort was also used for external validation. While the cohort demonstrated excellent predictive performance at 1 and 3 years, the limited number of patients surviving beyond 5 years resulted in an extremely low AUC at 5 years (Suppl1F). Consequently, we did not use its predictive results. Furthermore, compared to the training and internal validation set, GSE1450 cohorts exhibited more pronounced gender imbalances, which may have contributed to the performance differences of the NRSHC model across datasets. Therefore, we integrated clinical factors with the NRSHC model to construct a nomogram, further enhancing the predictive performance of the model.

NRSHC demonstrates predictive capabilities comparable to other HCC models constructed using machine learning, such as the prognostic features developed by Wang et al., Xing et al., and Xu et al. Compared to these models, NRSHC and its corresponding nomogram achieve higher AUC values, C-index, and clinical benefits in the TCGA cohort, showcasing superior predictive performance. Notably, NRSHC effectively stratifies HCC patients, identifying a high-risk group associated with poorer prognosis. The high-risk group also exhibits higher levels of immune infiltration and shows greater responsiveness to immune checkpoint inhibitors such as CTLA-4 and PD-1. In contrast, the low-risk group demonstrates higher sensitivity to HCC chemotherapy drugs. Overall, we propose that high-risk patients may derive greater benefit from immunotherapy. Additionally, we have integrated our findings into an online predictive tool, significantly enhancing the clinical value of this study.

This study also has certain limitations. First, we utilized multiple publicly available datasets from different databases, and the heterogeneity among these datasets may have introduced biases into the results. Although we have taken several measures to mitigate such effects, such as data standardization and normalization, feature dimensionality reduction, regularization techniques, the issue persists to some extent. Second, the NRSHC model performed less effectively in the external validation cohorts. Expanding the source and number of samples in the validation cohort, using multiple validation sets, and employing multiple rounds of cross-validation approach may help to improve the effectiveness of the model. Lastly, more experimental works are still needed to validate the specific functions of NQO1 in the TME. Due to limitations in funding, technology, and time, we are currently unable to address these issues. However, we will strive to overcome and address these limitations in future research.

In conclusion, this study reveals that NQO1 is overexpressed in HCC cells and may play a role in promoting the polarization of TAMs toward the M2 phenotype. We constructed an NRGs-associated HCC prognostic risk models a nomogram, and an online prediction tool at the first time. These tools can effectively predict the survival outcomes and immunotherapy responses of HCC patients, offering significant value in the clinical personalized treatment of HCC.