Introduction

Hepatocellular carcinoma (HCC) is a highly invasive malignancy that is challenging to treat and ranks fifth as the most prevalent malignant tumor globally.^1,2 It presents a significant global burden, with a particularly high incidence in Asia, notably in China.^3,4 HCC treatment options currently available encompass surgical resection, liver transplantation, transcatheter arterial chemoembolization (TACE),⁵ radiation therapy,⁶ targeted drug therapy,⁷ and immunotherapy.⁸ Surgical resection is considered highly effective, especially for patients diagnosed at an early stage, leading to improved survival rates. For patient’s ineligible for surgery, treatments such as TACE and radiation therapy are used to manage disease progression and alleviate symptoms.

Immunotherapy and targeted therapy have proposed recently promising results in HCC treatment. However, there are several limitations in the HCC diagnosis and treatment. First, patients are commonly diagnosed at advanced or intermediate-advanced stages, posing treatment challenges that contribute to poor prognoses. Second, some patients are unsuitable or unable to tolerate traditional treatment methods. Third, liver function impairment can impact cognitive function and the immune system, necessitating more precise and individualized treatment approaches.⁹ Thus, long-term follow-up and subsequent treatment are necessary due to high malignancy and recurrence rate of HCC. To address these challenges, future strategies for HCC treatment should focus on early screening, personalized diagnosis, and treatment.¹⁰ Combining innovative diagnostic and therapeutic methods could be key to improving treatment outcomes and patient survival rates.

Heat shock proteins (HSPs) are a diverse protein’s group that respond to changes across the cellular and extracellular environment.¹¹ They can be categorized into seven subfamilies in humans: Large HSP, HSP90, HSP70, HSP60, HSP40, HSP20, and HSP10. The expression levels and molecular weights of HSPs vary across different cell types and states.¹² HSPs have been extensively studied due to their involvement in various diseases pathogenesis, as autoimmune diseases, neurodegenerative diseases, and cancers. In neurodegenerative diseases, HSPs are believed to contribute to disease occurrence and progression by affecting peptide folding and protein degradation, leading to misfolding and accumulation of abnormal proteins.¹³ However, therapeutic approaches targeting HSPs in neurodegenerative diseases are still under development.¹⁴ In autoimmune diseases, research suggests that the interaction between HSPs and microbial antigens can modulate immune responses, influencing the likelihood and intensity of autoimmune reactions,¹⁵ which may contribute to improved treatment and prevention strategies for autoimmune diseases. In the context of tumors, HSPs play crucial roles in tumor development, affecting cancer cell proliferation, apoptosis, invasion and metastasis. In particular, HSP70 is considered a representative HSP associated with cancer.¹⁶ Our previous research showed that HSPA5 (GRP78) influences invasive capacity and liver tumor cells proliferation through its effects on NF-κB pathway¹⁷ and Wnt/HOXB9 pathway,¹⁸ thus impacting liver cancer development. The significance of HSPs in tumorigenesis positions them as potential therapeutic targets and prognostic indicators. Recently, novel HSP70 inhibitors like 2-phenylethynesulfonamide (PES),¹⁹ HSP90 inhibitor pimitespib,²⁰ and NCT-547²¹ have shown promising effects in improving the prognostic endpoints of patients when used alone or in combination therapy in several tumor types. However, their clinical usage still faces challenges.²²

In the study, we conducted a comprehensive multi-omics data analysis from TCGA and ICGC to investigate the correlation between HSPs expression and HCC patient’s prognosis. Our purpose was to develop a prognostic model based on significant HSPs and assess its accuracy in predicting patient outcomes, evaluating the immune microenvironment, and guiding treatment strategies, ultimately providing valuable insights to enhance the survival outcomes of HCC patients.

Materials and Methods

Data Acquisition and Processing

A list of seven subfamily members of HSPs, namely (Heat Shock 10 kDa, 20 Da, 40 kDa, 60 kDa, 70 kDa, and 90 kDa proteins, and Large HSP) was acquired from Genecard database (https://www.genecards.org/) (Table 1). Transcriptomic data and clinical datasets for HCC patients were acquired from The Cancer Genome Atlas (TCGA)²³ and the International Cancer Genome Consortium (ICGC).²⁴ The GSE104580 dataset included 147 hCC patients who underwent transarterial chemoembolization (TACE) treatment, comprising 81 responders and 66 non-responders. The GSE109211 dataset originated from Phase III STORM clinical trial (NCT00692770) and comprised 67 hCC patients who received sorafenib, with 21 responders and 46 non-responders.

Table 1 List of Heat Shock Proteins Family

Copy Number Alteration and Somatic Mutation Analysis

DNA alterations included mutations (truncation and missense) and copy number alterations (amplifications and deep deletions). Somatic mutation and copy number variation data of HCC patients were obtained from TCGA and ICGC, retaining non-silent mutations including Nonsense Mutation, Nonstop Mutation, Frame Shift Del, Frame Shift Ins, In Frame Del, and In Frame Ins. The distribution of HSPs mutations was analyzed using GenVisR, and connection between HSPs expression levels and copy number alterations was calculated using R software.

Functional Enrichment and Protein–Protein Interaction Network Construction Analysis

We utilized the online Search Tool for Retrieval of Interacting Genes (STRING) to construct a protein–protein interaction (PPI) network of HSPs.²⁵ To further explore the functional implications, enrichment analysis (FEA) was conducted utilizing Metascape.²⁶ HSP’s PPI network gene lists were imported, and the genes were clustered based on similarity using the cumulative hypergeometric distribution to calculate the P-value. Subsequently, functional enrichment analysis (FEA) of these genes was conducted utilizing Kyoto Encyclopedia of Genes and Genomes (KEGG) database.²⁷ To visualize PPI network, the network was imported into Cytoscape,²⁸ a software tool for visualizing and analyzing complex networks. To identify central genes within the network, the cyto-Hubba plugin was employed, which utilizes degree and maximum neighborhood component (MNC) topology algorithms to extract highly interconnected genes.

Construction of Heat Shock Protein Expression Risk Score (HRRS)

HSPs related to prognosis in both TCGA and ICGC cohorts, as well as HSPs exhibiting differential expression in both cohorts, a subset of 14 hSPs was identified, based on which the Heat shock protein–Related Risk Score (HRRS) was created utilizing LASSO Cox regression algorithm. The calculation of risk score formula is as follows:

The gene coefficient (ɛ) denotes the regression coefficient of each gene (ε), and gene expression (ε) corresponds to its expression value in each patient. HRRS was computed for each patient, and based on the median score, they were categorized into high-risk and low-risk groups. To enhance prediction accuracy and minimize the influence of confounding factors, patients across TCGA cohort (Table 2) who met following criteria were excluded: overall survival (OS) less than 3 months, inadequate clinical information, and non-R0 surgery. Univariate and multivariate Cox regression analyses were conducted to assess whether HRRS could serve as an independent prognostic factor. A forest plot was generated, incorporating the M stage, neoplasm status, and HRRS. The C-index and receiver operating characteristic (ROC) curve assessed HRRS and Cox regression model accuracy, respectively. Utilizing HRRS, a Kaplan–Meier analysis of survival was conducted, and the time ROC package was applied to construct 1-, 2-, and 3-year ROC curves.

Table 2 Clinical Characteristics of HCC Patients in TCGA Database

Relationship Between HRRS and HCC Immune Microenvironment

The differentially expressed genes (DEGs) analysis between two-risk groups was conducted utilizing “limma” R package. DEGs were subsequently subjected to gene set enrichment analysis (GSEA).²⁹ GSEA based on c2 (KEGG gene sets: c2.cp.kegg.v7.0.symbols.gmt) and h (hallmarks gene sets: h.all.v7.0.symbols.gmt). To further investigate immune characteristic differences across two-risk groups, comparisons were made regarding immune responses and major histocompatibility complex (MHC) genes’ expression. Additionally, relationship between HRRS and immune cell levels, including CD8 T cells, CD4 T cells, B cells, dendritic cells, macrophages, and neutrophils, was investigated.

Role of HRRS in Predicting the Therapeutic Effects of HCC

HRRS’s ability to predict TACE treatment response in the GSE104580 dataset and the efficacy of sorafenib in the GSE109211 dataset was assessed. In addition, immune checkpoint, Tumor Immune Dysfunction and Exclusion (TIDE) scores, a novel approach for assessing immune checkpoint blockade (ICB) efficacy, were obtained from the TIDE website (http://tide.dfci.harvard.edu/), following which the relationship between HRRS and the efficacy of immune therapy was evaluated using the TIDE prediction index. Furthermore, cell line data from Genomics of Drug Sensitivity in Cancer (GDSC) database³⁰ (https://www.cancerrxgene.org/) was obtained to predict chemotherapy response in cancer. The pRRopheticR package was employed to predict half-maximal inhibitory concentration (IC50) values of frequently utilized chemotherapeutic drugs among HCC patients from TCGA set, including drugs like paclitaxel and sunitinib, and the effectiveness of these drugs was evaluated between the two groups.

Tissue Specimens and Immunohistochemical Staining

Human HCC tissue and corresponding adjacent tissues were collected from patients diagnosed with hepatocellular carcinoma at No. 215 hospital of Shaanxi Nuclear Industry from June 2014 to January 2020. All patients signed informed consent forms, and the project has been approved by the Ethics Committee of No. 215 hospital of Shaanxi Nuclear Industry (No. 2024(039)). Samples were collected at 12, 3, 6 and 9 points, respectively, and at least one piece of tumor tissue was collected for molecular pathology examination. Liver tissue within <1 cm (para-cancer liver tissue or cut edge) and >1 cm (para-cancer liver tissue or cut edge) from the tumor edge were collected respectively. Immunohistochemical staining (IHC) was performed as previously described.¹⁸

Statistical Analysis

Statistical analyses were conducted utilizing R software, V. 4.1.3 (https://www.r-project.org/), and SPSS software, V. 26.0 (SPSS Inc., Chicago, IL). Univariate Cox analysis evaluated the impact of clinical parameters and HSP expression on HCC patient survival. LASSO and multivariate Cox regression analyses constructed the HRRS model utilizing HSPs. Additionally, multivariate Cox regression analysis assessed HRRS independence and other clinical indicators as prognostic factors. Statistical significance was denoted as * P < 0.05, ** P < 0.01, and *** P < 0.001.

Result

Landscape of Somatic Alterations in HSPs

To investigate HSPs genomic characteristics among HCC, we analyzed the mRNA expression and HSPs mutation frequency in HCC patients from ICGC and TCGA cohorts. Our analysis revealed that most HSPs exhibited elevated expression levels in liver cancer tissues, as observed in both large-scale databases (Figure 1A and B). Somatic DNA mutations in HSPs were detected in around 20% of patients, with an overall mutation frequency ranging from 2% to 7% (Figure 1C and D). Missense mutations were found to be the most common type of mutation. Furthermore, HSPs were positively correlated with copy number variations (CNVs) (Figure 1E). Further correlation analysis between HSP expression levels and CNVs identified 20 hSPs with a significant positive correlation (Figure S1).

Figure 1 The genomic characteristics of HSPs in HCC. The expression of HSPs in ICGC (A) and TCGA (B) databases. Somatic DNA mutation of HSPs in ICGC (C) and TCGA (D) databases. Copy number variation of HSPs in hepatocellular carcinoma (E). *p < 0.05,**p < 0.01 and ***p < 0.001.

Functional Enrichment Analysis (FEA) of HSPs

After analyzing the genomic differences in heat shock proteins (HSPs), we conducted functional enrichment analysis using the STRING database to identify significantly associated neighboring genes of the HSPs (Figure 2A). Gene functional enrichment analysis was performed utilizing Metascape database, revealing enrichment in functions related to protein folding and cellular stress response (Figure 2B). Furthermore, we conducted GO and KEGG pathway enrichment analyses utilizing Database for Annotation, Visualization, and Integrated Discovery (DAVID). The biological process analysis demonstrated significant enrichment in protein folding and response to topologically incorrect proteins (Figure 2C). The cellular component analysis revealed significant enrichment in the chaperone complex and the serine/threonine protein kinase complex (Figure 2D). Additionally, molecular function analysis indicated significant enrichment in heat shock protein binding and chaperone binding (Figure 2E). The KEGG pathway analyses demonstrated significant enrichment in protein processing in the endoplasmic reticulum, lipid and atherosclerosis (Figure 2F). Through degree and MNC calculations using Cytoscape 3.7.1, we identified the top 15 genes, which included HSPs, TP53, and AKT1 (Figure 2G and H), among which TP53 and AKT1 have been widely studied.

Figure 2 Functional enrichment analysis of HSPs. Neighbor genes significantly related to HSPs identified through the STRING database (A), HSP-related genes were analyzed by Metascape database (B), HSP-related genes were analyzed by DAVID database (C-F), Using Cytoscape 3.7.1 to display PPI results through degree and maximum neighborhood component (MNC) calculations (G&H).

Construction of an HSPs-Related Prognostic Model for HCC Patients

By analyzing the association between HSPs and liver cancer patient’s prognosis, we found 26 hSPs associated with the prognosis in the ICGC database (Figure 3A), and 24 hSPs in the TCGA database (Figure 3B). Then, we selected 14 hSPs (Figure 3C) that showed statistically significant expression and prognosis in both databases for further investigation. A prognostic model for HCC patients with overall survival (OS) data was developed in the ICGC cohort. Multivariate and LASSO Cox regression analyses were utilized to assess optimal prognostic features from the 14 hSP genes. After minimizing lambda in the LASSO Cox regression model, four HSPs were chosen to develop HSPs-Related Prognostic Model (Figure 3D and E). The HRRS was calculated using the formula:

Figure 3 Construction of an HSPs-Related Prognostic Model for HCC Patients. A total of 26 hSPs were associated with the prognosis of patients with HCC in ICGC database (A), 24 hSPs in TCGA database were related to the prognosis of patients with HCC (B), The expression levels and prognostic significance of HSPs were statistically significant in both the ICGC and TCGA databases (C), The LASSO regression model was employed in a 10-fold cross-validation to evaluate the partial likelihood deviance, the optimal values were determined based on both the minimum criterion and the 1-SE criterion, and these values are depicted by vertical dotted lines in the plot (D), LASSO coefficient profiles of four HSPs in the 10-fold cross-validation. The vertical dotted lines were drawn at the optimal values by using the minimum criteria and 1-SE criteria (E), Coefficients of the four HSPs in the cox regression model (F), The HR and P-value of four HSPs in the cox regression model (G), Patients were divided into high-risk and low-risk groups based on median level of HRRS in ICGC (H), Kaplan–Meier survival curve of OS between high-risk and low-risk groups in the ICGC (I). Patients were divided into high-risk and low-risk groups based on median level of HRRS in TCGA (J), Kaplan-Meier survival curve of OS between high-risk and low-risk groups in the TCGA (K).

HRRS = (0.00314 * HSPD1) + (0.03620 * DNAJC1) + (0.06078 * DNAJC5) + (0.01983 * DNAJC8) (Figure 3F).

Notably, the HSPD1 and DNAJC5 expression levels significantly contributed to poor prognosis (Figure 3G). The distribution of HRRS and survival status in the ICGC database were visually represented (Figure 3H). Kaplan–Meier survival analysis revealed that the overall survival (OS) was lower in the high-risk group than in the low-risk group in the ICGC database (Figure 3I). Similar results were observed in the TCGA database (Figure 3J and K).

HRRS as an HCC Patients Independent Prognostic Factor

To assess independent prognostic significance of HRRS risk score, which consists of HSPs, for HCC patients, univariate and multivariate Cox regression analyses were conducted considering various covariates, including Age, Gender, Neoplasm status, Family history, Pharmaceutical therapy, Ablation embolization, Neoplasm histologic Grading, T, N, M staging, Neoplasm disease stage, Vascular invasion, Child–Pugh classification grade, and Adjacent hepatic tissue inflammation extent type. The analysis revealed that Neoplasm status and HRRS were independent factors capable of predicting HCC patient’s prognosis (Figure 4A and B).

Figure 4 HRRS as an independent prognostic factor for HCC patients. Univariate and multivariate Cox regression analyses of the association between clinicopathological variables (including HRRS) and overall survival of patients (A and B), predicted the probability of HCC patients for 1-, 3- and 5-year OS in the TCGA cohort (C), The calibration chart of HCC patients for the 1-, 3- and 5-year OS in the TCGA cohort (D), The time-dependent c-index graph of HRRS, Neoplasm Status and M Stage (E), HRRS decision curve analysis of time-related risks in ICGC and TCGA (F&G).

To visually represent the prognostic parameters, a nomogram was constructed with variables having a p-value less than 0.1 from multivariate Cox analyses. Each patient was assigned a total score by summing each prognostic parameter scores. Higher total scores were associated with poorer patient outcomes (Figure 4C). Calibration plots for 1-, 3-, and 5-year survival probabilities demonstrated strong concordance between column chart predictions and actual observations (Figure 4D). Furthermore, TCGA-based time-dependent C-index curves showed that HRRS exhibited better predictive performance compared to other single factors (Figure 4E). The average AUC values for 1-, 2-, and 3-year prognosis predictions in ICGC cohort were 0.858, 0.751, and 0.711, respectively (Figure 4F). For TCGA cohort predictions, the average AUC values for 1-, 2-, and 3-year survival were 0.740, 0.792, and 0.786, respectively (Figure 4G).

Comparison of Clinical Pathological Features and Immune Microenvironment Between High-Risk and Low-Risk Groups

Next, we compared clinical pathological features differences between two-risk groups in TCGA cohort. No significant differences were found in diagnostic age, gender, N, M Stages, and Neoplasm Status between two-risk groups. However, statistically significant differences between two groups were observed in T stage, Neoplasm Histologic Grade, and Neoplasm disease stage (Figure 5A).

Figure 5 Comparison of clinical pathological features and immune microenvironment between high-risk and low-risk groups. The difference in clinicopathological characteristics between high- and low-risk groups (A), The biological process of enrichment of DEGs between the high- and low-risk groups analyzed by GSEA (B and C), GSVA was performed to analyze the immune-related pathways of the high- and low-risk groups (D), Gene expression of HLA genome between the high- and low-risk groups (E), The relationship between HRRS and the infiltration level of different types of immune cells (F-K).

GSEA was conducted on DEGs between two-risk groups to explore molecular mechanisms underlying these differences. GSEA analysis of HALLMARK gene sets demonstrated significant enrichment in Bile Acid Metabolism, Epithelial-Mesenchymal Transition, and Angiogenesis (Figure 5B). Notably, differences in Inflammatory Response were also assessed across two groups (Figure 5C).

To investigate the differences in immune infiltration patterns, Gene Set Variation Analysis (GSVA) was done, revealing distinct immune infiltration patterns between two-risk groups in TCGA cohort (Figure 5D). Additionally, the relationship between the high-risk and low-risk groups and MHC was examined, indicating higher expression levels of MHC in the high-risk group (Figure 5E). Finally, we assessed the relationship between HRRS and various immune cell infiltrate levels. The results indicated that as the HRRS increased, the infiltration levels of B cells, CD4 T-cells, CD8 T-cells, neutrophils, macrophages, and dendritic cells also increased (Figure 5F–K).

HCC Treatment Strategies Based on HRRS

To assess the utility of HRRS in guiding treatment strategies for HCC, we investigated its predictive value in HCC patients undergoing TACE using the GSE104580 dataset. This dataset comprised HCC patient’s mRNA expression arrays categorized into responders and non-responders TACE. Our analysis revealed that HRRS was significantly higher among TACE non-responder group than TACE responder group (P-value <0.0001, Figure 6A). Additionally, ROC curve analysis demonstrated that HSPs exhibited significant diagnostic value in predicting TACE response, with an AUC of 0.791 (Figure 6B). Furthermore, we examined another independent dataset (GSE109211) and observed that HRRS was significantly decreased in HCC patients who responded to sorafenib treatment (P = 0.0013, Figure 6C). ROC curve analysis indicated that HRRS had predictive value in identifying the effectiveness of sorafenib treatment, with an AUC of 0.649 (Figure 6D).

Figure 6 HCC treatment strategies based on HRRS. The box diagram of the HRRS between the TACE response group and the TACE non-response group in the GSE104580; ROC analysis of the HRRS used to predict the TACE response (A and B), The box diagram of HRRS between sorafenib treatment response group and sorafenib treatment non-response group in GSE109211; ROC analysis of HRRS to predict sorafenib treatment response (C and D), Box diagram of the TIDE score of the high- and low-risk group (E). Expression of five immune checkpoint molecules (PDCD1, HAVCR2, CD274, LAG3, and CTLA4) in the high and low groups (F), Chemotherapeutic responses to four chemotherapeutic drugs in the high- and low-risk groups (G–J).

To assess the HRRS potentiality as a predictor biomarker for clinical response to immune checkpoint inhibitor therapy, we conducted a comparative analysis of TIDE score distribution among different risk groups. Our findings demonstrated that high-risk group had lower TIDE scores (Figure 6E). Additionally, both two-risk groups were analyzed for five immune checkpoint molecules (PD1, HAVCR2, PD-L1, LAG3, and CTLA4). The findings revealed a significant upregulation of (PD1, HAVCR2, and CTLA4) among high-risk group (Figure 6F).

Moreover, GDSC data were utilized to predict drug treatment IC50 values for TCGA cohort HCC patients. Our analysis revealed that four drugs (Paclitaxel, Sunitinib, Pyrimethamine, Epothilone B) exhibited significantly reduced IC50 values across high-risk group, suggesting that HCC patients among high-risk group might be more sensitive to treatment regimens containing these drugs (Figures 6G–J & S2).

Finally, we retrieved protein expression data of the four HSPs comprising HRRS from Human Protein Atlas database. The analysis showed increased expression levels of HSPD1, DNAJC1, DNAJC5, and DNAJC8 in HCC tissues (Figure S3A-D). Further analysis of the protein expression levels of HSPD1 and DNAJC5 in HCC tissues and adjacent tissues by IHC also showed that there were significant differences in the expression of HSPD1 (Figure 7A and B) and DNAJC5 (Figure 7C and D) in HCC tissues and adjacent tissues.

Figure 7 Protein expression of HSPD1 and DNAJC5 in normal liver tissues and hepatocellular carcinoma tissues. The expression level of HSPD1 in hepatocellular carcinoma tissues was significantly higher than that in normal liver tissues (A), and the IHC score was significantly different (B); the expression level of DNAJC5 in hepatocellular carcinoma tissues was significantly higher than that in normal liver tissues (C), and the IHC score was also significantly different (D).*p < 0.05,**p < 0.01 and ***p < 0.001.

Discussion

HCC remains a major burden on human lives and resources. HSPs, a protein’s group implicated in cellular protection and repair, play crucial roles in various cellular stress responses. They can assist in refolding denatured proteins under heat stress and other forms of cellular stress, including oxidative stress, chemical stress, and pathogen infections.^31,32 HSPs are highly expressed in various human tumors, including HCC. Their elevated expression in cancer cells helps them adapt to adverse microenvironmental conditions and develop resistance to chemotherapy. HSPs can impact HCC occurrence and development via multiple pathways.¹⁶ However, the use of drugs targeting HSPs in cancer patients is currently hindered due to several challenges.²² Therefore, this study aimed to comprehensively investigate HSPs; impact on HCC tumor development, prognosis assessment, immune microenvironment alterations, tumor immune evasion, and treatment response.

We comprehensively analyzed transcriptomic changes in 77 hSPs in HCC using multiple large-scale transcriptomic datasets. Our findings revealed that most HSPs exhibited elevated expression levels in liver cancer tissues. However, the analysis of HSP mutations indicated that the overall mutation frequency of HSPs in HCC was relatively low, which aligns with the highly conserved nature of HSPs in evolutionary processes. Interestingly, we observed a close association between HSP expression levels and CNVs, suggesting that imbalanced HSP expression could be influenced by CNVs. Furthermore, by employing topological algorithms and bioinformatics data filtering, we identified two genes, AKT1 and TP53, which displayed close relationships with HSPs within the complex pathway network. These genes have significant functions in the tumors onset and progression, as alterations in AKT1 and TP53 have been implicated in HCC occurrence. Therefore, further investigation should concentrate on investigating targeted relationship between HSPs and AKT1³³ with TP53³⁴ to gain deeper insights into their interplay in HCC.

Then, we examined the association between HSPs and HCC patient’s prognosis. We identified 14 hSPs significantly linked to patient prognosis. According to LASSO and multivariate Cox regression analyses, we constructed HRRS, which is composed of HSPD1, DNAJC1, DNAJC5 and DNAJC8. HSPD1 belongs to the HSP60 subfamily. Kim Seon-Kyu et al³⁵ have shown that the downregulation and relocalization of HSPD1 in breast cancer affected apoptosis across cells. Parma Beatrice et al³⁶ found a close association between HSPD1 and poor prognosis in NSCLC patients and reported that reducing HSPD1 inhibited tumor cell proliferation in-vivo and in-vitro. Additionally, HSPD1 has been found to have a role in tumor cell survival, proliferation, migration and metabolism in esophageal cancer³⁷ and oral squamous cell carcinoma.³⁸ However, studies examining the role of HSPD1 in HCC tumorigenesis are scarce.

DNAJC1, DNAJC5 and DNAJC8 all belong to the HSP40 subfamily. In the present study, we found that DNAJC1, DNAJC5 and DNAJC8 were extremely expressed across liver cancer tissues and closely linked to patient prognosis. DNAJC1, also known as MTJ1, can interact with GRP78,³⁹ and studies have indicated that the binding of DNAJC1 to cell surface GRP78 stimulates cell proliferation.⁴⁰ DNAJC5, cysteine string protein alpha (CSPα), has been identified as a crucial mediator of misfolded protein secretion in MAPS. It can accompany misfolded proteins in MAPS to the extracellular space.⁴¹ In HCC research, Wang Hailong et al⁴² found that DNAJC5 boosts SKP2-mediated p27 degradation, a cyclin-dependent kinase inhibitor 1B, via promoting SKP2-p27 complex formation. Conversely, DNAJC5 Knockdown rescues SKP2-induced p27 level reduction. DNAJC8 regulates glycolysis and serves as a TIG1target. Wang Chun-Hua et al⁴³ discovered that silencing DNAJC8 reduces GLUT1 expression upregulation and glucose uptake induced by ectopic DNAJC8 expression, thus affecting the proliferation of cervical cancer cells.

Furthermore, we assessed each patient’s HRRS and stratified them into high and low HRRS groups. In the ICGC cohort, we observed a strong association between HRRS and OS, where low-risk patients had a significantly more favorable prognosis. Utilizing TCGA cohort, we conducted a survival study to further verify HRRS’s efficacy. Similar to ICGC cohort results, OS median for low HRRS patients was considerably higher than that high HRRS patients. Univariate and multivariate Cox regression analysis indicated that HRRS served as an independent HCC prognostic factor. Further HRRS analysis and clinical pathological parameters demonstrated that high HRRS patients had significantly higher T stage, Neoplasm Histologic Grade, and Neoplasm disease stage than low HRRS patients. To investigate underlying reasons for these differences, we employed GSEA. The results demonstrated significant enrichment of pathways related to Epithelial-Mesenchymal Transition and Angiogenesis in the high HRRS group, which are strongly related to HCC malignancy. Additionally, we demonstrated significant differences across Inflammatory Response between two-HRRS groups, following which HRRS and tumor immune cell infiltrate association were examined. The results revealed that higher HRRS patients had elevated levels of B cells, CD4 T cells, CD8 T cells, dendritic cells, macrophages and neutrophils. Among them, HRRS showed highest correlation with macrophages, the tumor microenvironment’s key component (TME). Considering that increased macrophage infiltration promotes HCC growth, progression and resistance to sorafenib, further exploration of how HSPs influence the enrichment of various immune cell populations across HCC TME would contribute to enhancing immunotherapy efficacy and improving survival outcomes of HCC patients.⁴⁴

Various strategies such as TACE, immunotherapy, targeted drugs and chemotherapy are used to treat advanced HCC.^45,46 In our study, we observed a significant increase in HRRS in HCC patients who did not respond to TACE treatment, demonstrating an area under the curve (AUC) greater than 0.7. Sorafenib, an oral multi-target tyrosine kinase inhibitor, is currently a first-line treatment for HCC.⁴⁷ Similarly, we also found a significant increase in HRRS in HCC patients who did not respond to sorafenib treatment, indicating that HRRS can serve as a predictor marker for HCC patients’ response to TACE and sorafenib. Furthermore, with the advancement of numerous studies on TME, co-administration of atezolizumab, an inhibitor of programmed death-ligand 1 (PD-L1), and bevacizumab, an antibody targeting vascular endothelial growth factor (VEGF), has demonstrated encouraging outcomes in enhancing prognosis of unresectable or metastatic HCC patients who have not received systemic treatment.⁴⁸ Our study revealed a strong correlation between HRRS and checkpoint molecules (PD-1, HAVCR2 and CTLA4) expression, suggesting that HSPs may have a function in immunotherapy prediction efficacy, emphasizing the importance of HRRS in potentially guiding immunotherapeutic approaches. Finally, we found a negative correlation between the IC50 values of certain drugs (ie, Paclitaxel, Sunitinib, Pyrimethamine, Epothilone B) and HRRS, suggesting that HCC patients with high HRRS may have a better response to these drugs. Thus, combining these drugs with existing treatment regimens may further enhance HCC patient’s prognosis.

Our study has a few limitations that merit consideration. Initially, analysis was relied on publicly available high-throughput sequencing datasets, and the lack of experimental validation restricts the robustness and reliability of our findings. Therefore, further, in vitro or in vivo experiments are necessary to verify HSPs role across HCC development and verify the prognostic significance of HRRS. Secondly, our study primarily focused on transcriptomic analysis, and integrating other omics data, such as proteomics and metabolomics, would provide a full comprehensive knowledge of immunological changes across HCC. Exploring the interplay between different layers of biological information could potentially enhance the immunotherapy efficacy among HCC patients.

In conclusion, our study provides evidence supporting the significant HSPs role in the development and progression of HCC. We have demonstrated that the HRRS, comprising specific HSP members, is a robust prognostic indicator for HCC. Furthermore, our findings reveal the correlation between HRRS and immune infiltrate patterns and for potential HCC treatment implications. The HRRS holds promise as a valuable tool for predicting HCC survival outcomes and guiding treatment decisions.