The mRNA expression data and clinical profiles were provided by the TCGA (The Cancer Genome Atlas) database (https://portal.gdc.cancer.gov), ICGC (International Cancer Gene Consortium) database (https://dcc.icgc.org/releases/current/Projects) and GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/). In the TCGA database, we collected 371 HCC samples and 50 normal samples. In the ICGC database, we collected 240 HCC samples and 202 normal samples. GSE14520 and GSE36376 datasets were used to analyze the expression changes of four key URGs in normal and HCC tissues. The TNMplot database (https://tnmplot.com/analysis/) was used to compare the expression of UPR-related signature in HCC tissues and normal tissues.

In the TCGA dataset, |Log₂(fold change)| > 0.58 and P < 0.05 were set as a threshold to screen the differentially expressed genes (DEGs) between HCC and normal samples. Then, we obtained genes associated with UPR in the MSigDB (http://software.broadinstitute.org/gsea/msigdb/index.jsp) and UALCAN (http://ualcan.path.uab.edu) databases, and intersected them with DEGs to finally obtain 31 URGs for establishing a UPR-related gene set. The R language packages “ggplot2”, “ggunchained” and “pheatmap” were used. Signaling pathway enrichment analysis was completed by the Metascape database (https://metascape.org/gp/index.html#/main/step1).

The R language package “ConsensusClusterPlus” was used to investigate the underlying molecular clusters based on the expression of 31 URGs. Principal component analysis (PCA) was used to evaluate the expression patterns of URGs. The R package “pheatmap” was used to draw a clustering heatmap. The R packages “survival” and “survminer” were used for Kaplan–Meier (KM) analysis.

Univariate Cox analysis was utilized to screen prognosis-related UPRs based on the TCGA and ICGC databases. Least absolute shrinkage and selection operator (LASSO) regional analysis was further carried out to build a risk scoring model with the R language package “glmnet”. In the risk scoring model, risk score = sum of coefficients × prognostic URG expression. The HCC samples were separated into high-risk and low-risk groups by means of the median risk score as the cutoff point. KM analysis was performed to assess the differences in prognosis using the log-rank test. Receiver operating characteristic (ROC) curves were generated using the R package “timeROC” to evaluate the accuracy of the risk scoring model.

The R language package “forestplot” was used to construct the forest plot and display variables such as P value, HR and 95% CI in regression analysis. The “rms” package was used to make nomogram diagrams based on URG expression, age tumor (T), node (N), metastasis (M) and grade. The predictive and discriminative ability of the nomogram for 1-, 3-, and 5-year OS was assessed using ROC and concordance index (C-index) using the R package “survival”. Calibration curves were used to evaluate the difference between predicted results and actual survival rate of patients.

The R language package “maftools” was applied to visualize the mutation maps of two clusters and generate the waterfall plots.

The R language package “limma” was used to select DEGs between cluster 1 and cluster 2. |Log₂(fold change)| > 0.58 and P <0.05 were used as thresholds. The “ClusterProfiler” package was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs between cluster 1 and cluster 2.

The R language package “immunedecov” was utilized to calculate the infiltrated abundance of immune cells in two clusters according to the XCell and TIMER algorithms. The “pheatmap” package was used to display the abundance of immune cells in each HCC sample.

We analyzed the difference in the expression of common immune checkpoint genes between cluster 1 and cluster 2, and used the “ggplot2” and “ggpubr” packages for visualization. The tumor immune dysfunction and exclusion (TIDE) algorithm (https://tide.dfci.harvard.edu) was used to evaluate the difference in immunotherapy response between two clusters.

Immunohistochemical staining (IHC) results of normal liver tissues and HCC tissues were obtained from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). To verify the expression of core URG expression in different cells in HCC tissues, single-cell RNA-sequencing results were obtained from the Human Liver Browser (https://itzkovitzwebapps.weizmann.ac.il/webapps/home/session.html?app=HumanLiverBrowser).

The Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) and R language package “pRRophetic” were used to evaluate the sensitivity of two clusters to chemotherapeutic drugs, and the half maximum inhibitory concentration (IC50) was calculated by ridge regression. The correlations between IC50 and the risk score and expression of core URGs were evaluated.

Human hepatocellular carcinoma cell lines HepG2 and Huh7 were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin in a humid atmosphere (37°C, 5% CO₂) as described previously (11, 23, 24). HepG2 and Huh7 cells were treated with different concentrations of sorafenib for 24 h. Total RNA was extracted using the RNA-easy kit (Vazyme, Nanjing, China). Reverse transcription kit (Biosharp, Beijing, China) was used to reverse transcribe RNA into cDNA. RT-PCR was performed using SYBR qPCR Master Mix (Biosharp, Beijing, China) to examine the expression levels of four target genes according to the manufacturer’s instruction. The primer sequences are shown in Table 1 . Relative mRNA level was calculated by 2^-ΔΔCT method normalized with internal reference gene S18.

All statistical analyses were completed by R software (version 4.2.2). The Wilcoxon test was used to compare the significant differences between two groups. The Kruskal–Wallis test was applied to estimate the differences among three groups. The log-rank test was used for KM analysis. Pearson correlation analysis was used to measure the correlation between the two groups. P < 0.05 was defined as a statistically significant difference.

To obtain differentially expressed UPR-related genes (URGs), we comprehensively considered the following three aspects. First, we collected mRNA expression data from 371 HCC samples and 50 normal samples in the TCGA database. A total of 7509 DEGs were selected according to |Log₂(fold change)| > 0.58 and P < 0.05 ( Figure 1A , S2A ). Second, we retrieved 113 and 250 URGs from the MSigDB and UALCAN databases, respectively. Finally, we selected the intersection of three gene sets as differentially expressed URGs, and a total of 31 URGs were collected ( Figure 1B ).

We investigated the expression of 31 URGs in HCC and normal samples using the TCGA database. Among 31 URGs, ARFGAP1, ATF4, ATF6, ATP6V0D1, BAG3, CALR, CXXC1, DCTN1, FKBP14, GOSR2, HSP90B1, HSPA5, HSPA9, HYOU1, KDELR3, NPM1, PDIA6, PREB, SHC1, SRPRB, SSR1, STC2, TATDN2, TSPYL2, WFS1 and WIPI1 were significantly upregulated in HCC samples ( Figures 1C, D and Figure S2B ). The expression of ATF3, CCL2 and IGFBP1 was obviously downregulated in HCC compared with normal samples ( Figures 1C, D and Figure S2B ).

To further validate the reliability of the above conclusion, we performed similar analysis using the ICGC database. ARFGAP1, ATF4, ATF6, ATP6V0D1, HYOU1, KDELR3, NPM1, PDIA6, PREB, SHC1, SRPRB, SSR1, STC2, TATDN2, WFS1, WIPI1, YIF1A and ZBTB17 expression was higher in HCC samples, while ATF3, CCL2, IGFBP1 and TSPYL2 expression was lower in HCC samples than normal samples ( Figures S3A , S3B ). Collectively, most URGs were upregulated in HCC samples.

To further understand the relationships, signaling pathways and molecular mechanisms of URGs, we used Metascape tool to perform signaling pathway enrichment analysis based on 31 URGs. Consistent with our speculation, these URGs were mainly enriched in a series of pathways related to ER stress, such as unfolded protein response (UPR) and cellular response to unfolded protein ( Figure S4 ). URGs were also associated with immunity-related pathway, such as IL-18 signaling pathway ( Figure S4 ).

The correlations among URGs were further calculated based on the TCGA and ICGC databases, and we observed a large number of significant correlations between URGs, of which positive correlations accounted for the majority ( Figures 1E, F ).

To apply URGs to personalized treatment of HCC patients, we classified the molecular subtypes of patients based on the expression levels of 31 URGs. We increased the number of clusters (k) from 2 to 6, and found that the well clustering effect was obtained when k = 2 ( Figure 2A ). Therefore, HCC samples were divided into two clusters ( Figure 2B ). The PCA analysis shows that the two clusters are well separated ( Figure 2C ). The expression of URGs in cluster 1 (C1) was generally higher than that in cluster 2 (C2) ( Figure 2D ). To determine the prognostic value of molecular classification, we then used KM analysis to analyze the overall survival (OS), progression-free survival (PFS), disease-free survival (DFS) and disease-specific survival (DSS) of HCC patients. Patients in C1 had worse prognosis than those in C2 ( Figures 2E–H ). Furthermore, HCC patients in C1 and C2 exhibited significant differences in terms of T stage, TNM stage and grade ( Table 2 ).

To provide a more comprehensive comparison between two molecular subtypes, we performed differential expression analysis, GO and KEGG analyses between C1 and C2. Upregulated genes and downregulated genes were identified based on |Log2(fold change)| > 0.58 and P < 0.05 ( Figures 3A, B ). The top five KEGG pathways for upregulated genes were Salmonella infection, human T−cell leukemia virus 1 infection, endocytosis, focal adhesion, and pathogenic Escherichia coli infection ( Figure 3C ). The top five GO pathways for upregulated genes were histone modification, covalent chromatin modification, RNA splicing, regulation of mitotic cell cycle phase transition, and regulation of cell cycle phase transition ( Figure 3D ). The top five KEGG pathways for downregulated genes were complement and coagulation cascades, metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450, chemical carcinogenesis-DNA adducts, and bile secretion ( Figure 3E ). Additionally, the top five GO pathways for downregulated genes were small molecule catabolic process, organic acid catabolic process, organic acid biosynthetic process, carboxylic acid catabolic process, and carboxylic acid biosynthetic process ( Figure 3F ). Of note, several immunity-related pathways were also enriched in the KEGG results, including TNF signaling pathway, Salmonella infection, pathogenic Escherichia coli infection, human T-cell leukemia virus 1 infection, and chemokine signaling pathway ( Figure 3C ).

Ferroptosis is a new type of programmed cell death and N6-methyladenosine (m6A) methylation is a common modification in eukaryotic mRNA. Both ferroptosis and m⁶A methylation play key roles in tumorigenesis, metastasis and immune escape in HCC. We then analyzed the expression of ferroptosis-related genes and m⁶A-related genes between C1 and C2. Most ferroptosis-related genes and m⁶A-related genes were upregulated in C1 compared with C2 ( Figures S5A , S5B ).

Because identified 31 URGs and DEGs between C1 and C2 were involved in immunity-related pathways ( Figure 3C ), we further evaluated the differences in the immune microenvironment between two subgroups. The TIMER algorithm was used to assess the abundance of infiltrating immune cells in C1 and C2 subgroups. Except for CD8+ T cells, the abundance of five major immune cells in C1 was significantly higher than that in C2 ( Figures 4A, B ). Among the immune microenvironments of all HCC samples, dendritic cells (DCs) accounted for the largest proportion ( Figure 4C ). To further verify the accuracy of the conclusion, we performed similar analysis using the XCell algorithm. Among the 28 types of immune cells, 16 types exhibited significant differences in infiltrated abundance of immune cells between C1 and C2 ( Figures S6A–C ). Additionally, the expression of immunosuppressive factors, including CD96, CSF1R, HAVCR2, IDO1, LAG3, LGALS9, NECTIN2, PDCD1, PDCD1LG2, TGFB1, TGFBR1, BTLA, CD160, CD244, CD274, CTLA4, IL10, IL10RB, TIGIT and VTCN1, was significantly higher in C1 than that in C2 ( Figure 4D ). To further estimate the responses to immunotherapy, TIDE algorithm was used. HCC patients with lower TIDE score may respond better to immunotherapy. Notably, the TIDE score of C1 was significantly higher than C2, indicating that HCC patients in C2 were more sensitive to immunotherapy compared with C1 ( Figure 4E ).

We obtained the mutation information of C1 and C2 from the TCGA database. The waterfall plot results demonstrated that the top five genes with the high frequency of alteration in C1 were TP53, TTN, CTNNBI, MUC16 and RYR2 ( Figure 5A ), and the top five genes with high frequency of mutation in C2 were CTNNB1, TTN, TP53, MUC16 and PCLO ( Figure 5B ). Missense mutation was the most frequent mutation type, in which single-nucleotide polymorphism (SNP) was the main part. SNV class showed that the most common variation was C > T in both subgroups.

Considering the clinical importance of chemotherapy, we then evaluated the sensitivity of 9 common chemotherapeutic drugs between two subgroups based on the GDSC database. The IC50 of 9 chemotherapeutic drugs, including sorafenib, sunitinib, vinblastine, doxorubicin, gemcitabine, 5-fluorouracil, etoposide, docetaxel and paclitaxel, in C1 was significantly lower than that in C2 ( Figure 6 ). These results imply that although HCC patients in C1 have a lower survival rate than patients in C2, fortunately, these patients were more sensitive to chemotherapy.

To further screen for URGs with prognostic value and better apply URGs to the clinical diagnosis of HCC patients, we performed a univariate Cox regression analysis and selected genes with prognostic potential (P < 0.05) ( Figures 7A, B ). 22 genes were selected from the TCGA dataset and 9 genes were selected from the ICGC dataset. Finally, 8 genes, including SRPRB, PDIA6, GOSR2, ATF4, WIPI1, HSPA5 and HYOU1, were selected from the intersection for subsequent analysis ( Figure 7C ).

We then used the optimum λ value to reduce the dimension of 8 genes by LASSO and established a prognostic risk model. According to the LASSO algorithm, we identified 4 core prognostic genes: ATF4, PDIA6, GOSR2 and SRPRB ( Figures 7D, E ). Based on the expression levels of 4 URGs, the risk score was calculated through the following formula: risk score = (0.0196×ATF4) + (0.1309×GOSR2) + (0.0889×PDIA6) + (0.4707×SRPRB). According to the median risk score, we divided HCC patients into high-risk group and low-risk group. The risk score, survival time and expression of four URGs in each patient were shown ( Figure 7F ). KM analysis results indicated that the survival rate of high-risk group was significantly lower than low-risk group ( Figure 7G ). We then used receiver operating characteristic (ROC) analysis to test the reliability and sensitivity of this prognostic model. The 1-year survival rate of areas under the ROC curve (AUC) was 0.739, the 3-year survival rate of AUC was 0.696, and the 5-year survival rate of AUC was 0.658 ( Figure 7H ). The above results indicate that the prognostic model we established is reliable.

To make our model more convincing, we further established a multivariate Cox regression analysis model using the four core genes using the ICGC database ( Figure S7A ). The survival rate of patients with high risk scores was worse than that of patients with low risk scores ( Figure S7B ). The AUC scores in the validation cohort were 0.815, 0.693 and 0.678 for 1-, 2-, and 3 years, respectively ( Figure S7C ).

To further verify the prediction accuracy of the 4-URG prognostic model, we carried out univariate and multivariate Cox regression analyses. Univariate Cox regression analysis indicated that ATF4, GOSR2, PDIA6, SRPRB, T stage, M stage and grade were significantly correlated with OS ( Figure S8A ). Multivariate Cox regression analysis suggested that SRPRB and T stage were markedly correlated with OS ( Figure S8B ). According to the results of multivariate Cox regression analysis, we used SRPRB and T stage to build a nomogram model to predict the OS of HCC patients at 1-, 3-, and 5 years ( Figure S8C ). Calibration plots showed agreement between actual and predicted survival rates ( Figure S8D ).

To assess the prognostic value of the 4-URG risk model in clinical practice, we analyzed the relationships between risk scores and various clinicopathological parameters, including age, sex, early grade (G1+G2) and advanced grade (G3+G4), early stage (T1+T2) and advanced stage (T3+T4), M0 and N0, TNM stage I+II and TNM stage III+IV. In the above clinicopathological parameters, the OS of high-risk group was worse than that of low-risk group ( Figure S9 ). These results proved that 4-URGs have good potential application value in the clinical prognosis of HCC.

To further verify the expression pattern of four URGs in liver cells, we collected single-cell RNA-sequencing data from the Human Life Browser database. Four URGs were expressed not only in HCC cells but also in immune cells ( Figure 8A ). The expression of 4-URG signature was not only upregulated in HCC tissues, but also much higher in metastatic tissues ( Figure S10A , S10B ). We also analyzed the expression differences of four URGs in the GSE36376 and GSE14520 datasets. As expected, the expression of four URGs was markedly elevated in HCC samples ( Figure S10C , S10D ). The HPA database was then used to compare the protein expression of four URGs. IHC results showed that the protein levels of ATF4, GOSR2, PDIA6 and SRPRB in HCC tissues were obviously higher than normal tissues ( Figure 8B ). We also collected immunofluorescence results in the HPA database to analyze subcellular localization. The results demonstrated that ATF4 was mainly located in the nucleus, GOSR2 was located in the Golgi apparatus, PDIA6 was located in the ER, and SRPRB was mainly located in the ER membrane ( Figure 8C ).

Sorafenib is the first-line drug for the treatment of advanced liver cancer. We then analyzed the relationship between sorafenib sensitivity and expression of 4 URGs. The IC50 of sorafenib was significantly negatively correlated with the expression of 4 URGs ( Figure 9A ). The HCC patients with higher expression of 4 URGs had lower IC50 values, indicating that patients with high expression of 4 URGs were more sensitive to sorafenib treatment ( Figure 9B ).

To further confirm whether sorafenib could alter the expression level of four URGs, HepG2 and Huh7 cells were treated with different concentrations of sorafenib for 24 h. Consistent with previous results, cell viability was decreased in a dose-dependent manner ( Figures 9C, D ). Additionally, the expression of GOSR2, PDIA6 and SRPRB was gradually decreased, but the expression ATF4 was significantly increased in HepG2 and Huh7 cells following sorafenib treatment ( Figures 9E, F ).