Integrative Analysis of scRNA-Seq and Bulk RNA-Seq Identifies Plasma Cell Related Genes and Constructs a Prognostic Model for Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the most common malignant tumor, posing a severe threat to human health.^1,2 Due to its subtle onset, most patients are typically diagnosed at an advanced stage, making them ineligible for surgical intervention and leading to a very poor prognosis.³ Recently, immunotherapy has emerged as a promising treatment for HCC, with rapid advancements in the field.⁴ The IMbrave150 study demonstrated that the combination of atezolizumab and bevacizumab increased the median overall survival (OS) to 19.2 months, compared to 13.4 months with sorafenib.⁵ Despite these advancements, the prognosis for HCC remains poor, and the mechanisms of tumorigenesis remain unclear. Therefore, studying the molecular mechanisms underlying HCC can help improve the poor prognosis of patients. Recent studies have shown that the complexity of the tumor immune microenvironment (TIME) and intra-tumor heterogeneity contribute to immunotherapy resistance and poor prognosis in HCC patients.⁶ Moreover, biomarkers can predict patient prognosis and guide clinical treatment by identifying heterogeneity in the TIME.^7,8 Therefore, developing new biomarkers can help understand the complexity of the TIME and tumor heterogeneity, ultimately enabling precise prognostic evaluation and guiding personalized treatment for HCC.

In recent years, the roles of various immune cells within the TIME, such as T cells, macrophages, and neutrophils, have been extensively explored in HCC, highlighting their potential as innovative therapeutic targets.^9–11 However, the study of plasma cells remains insufficiently understood. Plasma cells, recognized primarily for secreting antibodies, play a multifaceted role in regulating the TIME. For instance, antibodies produced by plasma cells can kill malignant cells through antibody-dependent cell-mediated cytotoxicity (ADCC).¹² Meanwhile, studies have shown that antibodies derived from the TIME are associated with a favorable prognosis in many types of cancer.¹³ Recently, Ma et al elucidated the heterogeneity of plasma cells within tumors, identifying ten distinct subgroups, each potentially playing unique roles in various types of tumors.¹⁴ Therefore, understanding the function and role of plasma cells may help develop new therapeutic strategies to improve the prognosis of cancer patients.

Bulk RNA sequencing (bulk RNA-seq) technology has been the primary tool for exploring the heterogeneity of the TIME and the mechanisms of tumorigenesis over the past decade. Numerous biomarkers have been identified as prognostic factors in patients with cancer, aiding clinical risk stratification and therapy guidance. For example, dynamic changes in the TIME are analyzed based on the overall expression levels of specific genes, and prognostic biomarkers have been developed to predict bladder cancer outcomes.¹⁵ However, bulk RNA-seq measures the average gene expression across the entire transcriptome of all cells, which overlooks potential cellular heterogeneity.¹⁶ Single-cell RNA sequencing (scRNA-seq) is a revolutionary genomics technology that enables comprehensive analysis of gene expression profiles at the single-cell level, revealing the roles of different cell subsets and intercellular communication within tumors.¹⁷ For example, Meng et al found that tumor cells can reprogram CD10+ALPL+ neutrophils through the NAMPT-NTRK1 signaling axis, leading to immune resistance and HCC progression.¹⁸ Therefore, it is necessary to develop biomarkers that guide clinical decision-making and provide personalized therapy by integrating scRNA-seq and bulk RNA-seq.

In this study, we integrated scRNA-seq and bulk RNA-seq analyses to identify plasma cells as the main cell subset involved in hepatocarcinogenesis and constructed a prognostic model based on PCRGs. The stability and reliability of the model were validated using external datasets and a nomogram. Additionally, we explored the molecular signaling pathways of the PCRGs model in HCC pathogenesis and their relationship with immune infiltration to better understand the mechanisms underlying the poor prognosis of HCC. We also elucidated the biological function of the key model gene, SSR3, in HCC in vitro. In summary, our results identified plasma cells as key clusters in HCC development. In addition, a prognostic model based on PCRGs can accurately predict the prognosis of patients with HCC and guide clinical treatment.

Materials and Methods

The overall design of the study is shown (Figure 1).

Figure 1 The flow chart of this study.

Data Sources and Access

The GSE189935 dataset, which included single-cell expression profiles for nine tumor samples and nine adjacent liver samples, was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov). The GSE14520 dataset, containing the expression profile data for 221 patients, was also obtained from the GEO database. Additionally, TCGA-LIHC expression data for 424 patients were acquired from TCGA database (https://portal.gdc.cancer.gov/).

ScRNA-Seq Analysis

First, the scRNA-seq data were processed using the “Seurat” R package, with abnormal samples filtered based on Unique Molecular Identifier(UMI) counts, gene counts, and mitochondrial gene ratios for each cell.¹⁹ The data were then standardized, normalized, and reduced using principal component analysis (PCA). The optimal number of principal components was determined via an elbow plot and unified manifold approximation and projection (UMAP) was used for nonlinear dimensionality reduction to delineate the positional relationships between clusters. Cell types and corresponding marker genes in the tissue were annotated using CellMarker and PanglaoDB databases.^20,21

The Disease Contribution of Distinct Cell Subsets

We characterized the contribution of different cell subsets to HCC by analyzing changes in cell numbers and gene expression.²² The signature genes for each group were initially identified. Specifically, the bulk differential gene expression analysis was performed, and the 231 differentially expressed genes (|Fold Change| > 1.5) were selected to characterize the changes between the “Normal” and “Tumor” groups. In addition, we independently analyzed differentially expressed genes between the “Normal” and “Tumor” groups within each subset. Genes that overlapped with 231 differentially expressed genes were selected as signature genes for each subset. We then established the FCscore to quantify both the quantity and expression level changes of signature genes throughout biological processes. . FCexp(i, j) represents the fold change of the ith signature gene in the jth cluster and FCprop(j) is the proportion of the fold change of the jth cluster. Finally, the average FCscore of all feature genes within each subset was used to assess the contribution of different cell subsets to HCC.

Construction and Validation of the Prognostic Model

Prognostic genes were used to create a prognostic model using the LASSO regression. A risk score formula for individual patients was devised by integrating the expression levels of each gene adjusted by the regression coefficients from the LASSO analysis. The risk score was calculated as follows: . Patients were stratified into low-risk (LRisk) and high-risk (HRisk) groups based on the median risk score Survival differences between the two groups were evaluated using Kaplan-Meier analysis, and statistical comparisons were made using the Log rank test. The predictive capacity of the risk score was investigated through lasso regression and stratified analysis. The efficacy of the prognostic model was assessed using receiver operating characteristic (ROC) curves.

Nomogram Construction and Clinical Relevance Analysis

A nomogram was used to evaluate the clinical characteristics and risk scores of patients. We constructed a nomogram based on the risk scores from the PCGR features and clinical variables to estimate the 3-year and 5-year OS probabilities for patients with HCC. The predictive performance of the nomogram was evaluated using concordance index and calibration curves. In our study, we combined the risk scores with sample characteristics to explore their clinical relevance.

GSVA Analysis and GSEA Analysis

Gene Set Variation Analysis (GSVA) is a non-parametric unsupervised method used to assess the enrichment of gene sets within transcriptomes. In this study, gene sets from the Molecular Signatures Database were used, and the GSVA algorithm was applied to score each gene set systematically to assess potential alterations in biological functions across various samples. Gene Set Enrichment Analysis (GSEA) was employed to conduct a more in-depth examination of disparities in signaling pathways between the high-risk and low-risk groups. Reference gene sets were obtained from the Molecular Signatures Database (https://www.gseamsigdb.org/gsea/msigdb/).

miRNA Network Construction

MicroRNAs(miRNAs) are small non-coding RNAs that regulate gene expression by promoting mRNA degradation or inhibiting mRNA translation.²³ Thus, we examined whether miRNAs regulate the transcription or degradation of some risk genes by targeting important genes. In this study, we obtained miRNAs associated with important genes from the miRcode database and visualized the miRNA-gene network using Cytoscape (version 3.8.0).

Transcriptional Regulation Analysis of Important Genes

This study used the R package “RcisTarget” to predict transcription factors based on motifs. The normalized enrichment score (NES) of a motif depended on the total number of motifs in the database. In addition to the motifs annotated by the source data, we inferred additional annotations based on the motif similarity and gene sequences. To estimate the overexpression of each motif in the gene set, the area under the curve (AUC) for each motif-gene set pair was calculated by computing the recovery curve of the gene set sorted by motifs. The NES of each motif was calculated based on the AUC distribution of all motifs in the gene set.

Tumor Immune Cell Infiltration Analysis

The CIBERSORT method is commonly used to assess immune cell populations within the TIME.²⁴ In this study, the CIBERSORT algorithm was applied to analyze patient data, estimate the relative proportions of 22 distinct types of immune infiltrating cells, and investigate the relationship between gene expression and immune cell composition.

Drug Sensitivity Analysis and the Response to Immunotherapy

Using the extensive pharmacogenomics database GDSC (Genomics of Drug Sensitivity in Cancer, https://www.cancerrxgene.org/), we employed the R package “pRRophetic“ to predict the sensitivity of individual tumor samples to chemotherapy. Regression methods were applied to derive IC50 estimates for distinct chemotherapy agents, and the accuracy of regression and prediction was validated through a 10-fold cross-validation using the GDSC training dataset.²⁵ All parameters were set to default values, including the “combat” parameter for batch effect removal and the average value for repeated gene expressions. Additionally, we utilized the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) to predict the impact of immunotherapy on the low-risk (LRisk) and high-risk (HRisk) groups.²⁶

Cell Culture

HCC cell lines (Hep3B, HepG2, Huh-7, MHCC97L, MHCC97H, SUN449, and HCCLM3) and the normal human liver cell line (LO2) were obtained from the National Collection of Authenticated Cell Cultures (Shanghai, China) with short tandem repeat certifications. All HCC cell lines were cultured in DMEM (HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (ExCell Bio, New Zealand) at 37 °C in a humidified atmosphere of 5% CO2.

Clinical Specimens

The clinical specimens used in this study were approved by the Ethics Committee of Zhujiang Hospital, Southern Medical University (Approval Number: 2023-KY-199-01). All specimens were collected in accordance with relevant ethical guidelines and with informed consent. In this study, they were used solely for their intended scientific purposes, without additional collection or any impact on patient privacy.

RNA Isolation and Quantitative Real-Time PCR (qRT-PCR)

Total RNA was isolated from cells using a reagent kit (GOONIE, China, Cat#400-100) according to the manufacturer’s instructions. RNA was reverse-transcribed into cDNA using a reverse transcription kit (TIANYA BIO, China, Cat#P1504). qRT-PCR was performed using SYBR Green Master Mix on a real-time PCR system (Roche LightCycler 480; Basel, Switzerland). Gene expression levels were normalized to GAPDH as an internal control and calculated using the 2 − ΔΔCt method. Primer sequences are listed in Table S1.

Transfection

The designed sgRNAs and control-sgRNAs were acquired from Qingke (Beijing, China) to knock out the SSR3 gene. The sgRNA sequences for SSR3 were as follows: SSR3-1F#CACCGAGGACTTGGCCGAGAGATTG, SSR3-1R#AAACCAATCTCTCGGCCAAGTCCTC, SSR3-2F#CACCGCTCCAAACAGCAGTCTGAGG, and SSR3-2R#AAACCCTCAGACTGCTGTTTGGAGC. Cells were transfected with plasmids or sgRNAs using Lipofectamine 3000 (GLPBIO, USA, Cat#GK20006) according to the manufacturer’s protocol. Transfection efficiency was assessed using qRT-PCR.

CCK8 Assay, Colony Formation Assay, and Transwell Migration Assay

Cell proliferation was assessed using the CCK8 assay. Cells were plated in 96-well plates and exposed to the CCK8 reagent for 2 h. Absorbance was recorded at 450 nm. For the colony formation assay, cells were seeded in 6-well plates and allowed to grow for two weeks. The colonies were then fixed with methanol, stained with crystal violet, and quantified. Transwell migration assays were conducted using Transwell chambers with an 8 µm pore size membrane. Cells were introduced into the upper chamber in serum-free medium, and medium containing 20% FBS was introduced into the lower chamber. After a 24-hour incubation, migrated cells were fixed, stained, and quantified.

Statistical Analysis

Survival curves were constructed using the Kaplan-Meier method and compared using the Log rank test. Multivariate analysis was performed using the Cox proportional hazards model. All statistical analyses were performed using R software (version 4.3.0), with the significance level set at p < 0.05.

Results

ScRNA-Seq Analysis Reveals the Contribution of Plasma Cells to HCC Development

In the present study, we obtained scRNA-seq data with complete single-cell expression profiles from the GSE189935 dataset. We used the “Seurat” R package to read the expression profiles. After quality control, we obtained a total of 5095 cells (Figure 2A–D). The data were processed through normalization, PCA, and harmony analysis for subsequent analysis (Figure 2E–G). UMAP analysis was used for cell clustering and visualization, and each cell cluster was annotated using databases such as CellMarker and Panglaodb,^20,21 resulting in the identification of 17 cell clusters and 14 different cell types (Figure 3A–C). A bar graph was used to display the proportion of cells in different groups (Figure 3D), revealing significant differences in the composition of immune cells between normal liver and tumor tissues. We then characterized the contribution of different cell subsets to HCC by considering changes in cell numbers and gene expression.²² Meanwhile, we selected 231 differentially expressed genes (|FC| > 1.5) to describe the changes in this process, and the results showed that plasma cells were found to contribute the most to HCC (Figure 3E). This finding was further validated using the TCGA-LIHC dataset (Figure 3F).

Figure 2 scRNA-seq analysis. (A–C) Following the quality control of scRNA-seq, 5,095 core cells were recognized. (D) The variance diagram shows the variation of gene expression in all cells of HCC. The red dots represent highly variable genes and the black dots represent non-variable genes. (E–G) PCA and harmony analysis.

Figure 3 Identification of plasma cells. (A) UMAP plot colored by various cell clusters. (B) UMAP plot colored by subpopulation of cells after annotation. (C) The image presents the expression levels of marker genes across various cell clusters. (D) The bar chart shows the proportion of different cell subpopulations in normal versus tumor tissue. (E) This image shows the contribution of different cell subpopulations to HCC. (F) Relative levels of plasma cell immune infiltration in the TCGA-LIHC dataset. **** represents p < 0.0001, ns, no significance in comparison with the control group.

Construction and Validation of PCRGs Prognostic Model

First, we extracted marker genes of the plasma cell population from scRNA-seq data using the “FindAllMarkers” function (p-value adj < 0.05 and |FC| > 1.5). Clinical information of HCC patients was collected from TCGA database, and 13 prognostic-related genes were identified using univariate Cox regression analysis (p < 0.05) (Figure 4A). LASSO regression was further screened for the characteristic genes in HCC (Figure 4B and C), and identified eight characteristic genes (Figure 4D). The risk score formula for the model is: Risk Score =CD79A×(−0.108603890133031)+SPCS2×0.00660846996462222+TRAM1×0.0148968367803627+PRDX4×0.0387216010891268+RPN2×0.0510535584736123+SSR3×0.0808719744617527+SEC61A1×0.136991700234778+ST6GALNAC4×0.158734118556289. The TCGA-LIHC dataset with survival data from TCGA database was randomly divided into a training set and a testing set at a 4:1 ratio. The patients were divided into HRisk and LRisk groups based on the median risk score Kaplan-Meier survival analysis indicated that the OS of the HRisk group was significantly lower than that of the LRisk group in both the training and testing sets (Figure 4E and F). Similar results were obtained in an independent external cohort, GSE14520. (Figure 4I) Furthermore, the ROC curve for OS in the TCGA training set and internal testing set showed AUC values greater than 0.70 at 1, 3, and 5 years (Figure 4G and H), indicating the good predictive performance of the risk model. The GSE14520 cohort was used as an independent group to validate the prognostic accuracy, with AUC values of 0.62, 0.61, and 0.63 at 1, 3, and 5 years, respectively (Figure 4J). Overall, our PCRGs prognostic model showed a good predictive ability for HCC prognosis.

Figure 4 Construction and validation of PCRGs prognostic model. (A)) Univariate Cox regression analysis showed that 13 PCRGs were associated with HCC prognosis. (B and C). LASSO regression of 13 PCRGs. (D) This image shows the coefficients of different genes in the model. (E–J) KM analysis and ROC curves in different datasets.

Clinical Relevance Analysis and Nomogram Construction

By analyzing the relationship between risk scores and clinical characteristics such as stage, T stage, and N stage (Figure S1A–H), we found that the risk score effectively classified HCC patient samples. Additionally, both univariate and multivariate analyses demonstrated that the risk score serves as an independent prognostic factor for patients with HCC. (Figure 5A and B). To confirm the precision of the risk score model, we performed regression analysis by dividing the TCGA-LIHC dataset into HRisk and LRisk groups based on the median risk score, showing significant contributions to the nomogram predictive model (Figure 5C). The 1-year, 3-year, and 5-year OS predictions demonstrated the high predictive capability of the risk score model for HCC patients with OS (Figure 5D–F).

Figure 5 Construction of nomogram model. (A)Univariate Cox analysis of risk scores and clinical characteristics. (B) Multifactorial Cox analysis. (C) Construction of the nomogram model. (D) The calibration curve of the nomogram. (E) The AUC of the prediction of 1, 3, and 5-year survival rates of HCC. (F) The DCA curve of the nomogram.

GSVA and GSEA Analysis

Our study examined the distinct signaling pathways associated with the HRisk and LRisk groups to elucidate the potential molecular mechanisms underlying the impact of risk score on tumor progression. GSVA findings revealed that the differential pathways between the two groups were predominantly enriched in mTORC1 signaling and PI3K/Akt/mTOR signaling (Figure 6A). GSEA results showed involvement of pathways such as the Fanconi anemia pathway, DNA replication, and cell cycle (Figure 6B).

Figure 6 Exploration of the mechanism behind the PCRGs prognostic model. (A) GSVA analysis. (B) GSEA analysis. (C)miRNA network construction. (D) Transcriptional regulation analysis.

miRNA Network Construction and Transcriptional Regulation Analysis

MicroRNAs (miRNAs) are small non-coding RNAs that regulate gene expression by promoting mRNA degradation or inhibiting mRNA translation.²³ We analyzed to determine if specific miRNAs play a role in regulating the transcription or degradation of key risk genes. Using the miRcode database to perform reverse prediction on eight model genes, we identified 79 miRNAs and 225 mRNA-miRNA interactions. These relationships were visualized using Cytoscape (Figure 6C). Further analysis revealed that these eight model genes were regulated by multiple transcription factors. Enrichment analysis of these transcription factors, motif-TF annotation, and selection analysis of important genes showed that cisbp _ M1670 had the highest normalized enrichment score (NES: 7.73) among the motifs. We displayed all the enriched motifs and the corresponding transcription factors for the model genes (Figure 6D).

Immune Characteristics of Risk Score and Immunotherapy Response

To evaluate immune cell infiltration in the TIME, we analyzed the proportion of immune cells among the different risk groups (Figure 7A). The findings revealed that naive B cells, M1 macrophages, resting mast cells, monocytes, resting NK cells, T cells CD4 memory resting, and others were significantly reduced in the HRisk group, while M0 Macrophages and regulatory T cells (Tregs) were significantly increased (Figure 7B). Studies indicate that Treg cells predominantly create an immunosuppressive microenvironment by inhibiting anti-tumor immune responses, suggesting that TIME in the high-risk group may be in an immunosuppressive state. The tumor Immune Dysfunction and Exclusion (TIDE) score is used to evaluate the response of cancer patients to immunotherapy. These results indicated that the HRisk group was more likely to undergo immune escape, reducing the effectiveness of immunotherapy. (Figure 7C). Additionally, we presented the immune cell infiltration and TIDE scores of the PCRGs (Figure S2A–H). In summary, this finding suggests that the HRisk group is in an immunosuppressive microenvironment and has poor immunotherapy efficacy.

Figure 7 Immunoinfiltration analysis, immunotherapy response prediction and drug sensitivity analysis. (A and B): Immunoinfiltration analysis in the HRisk and LRisk group. (C) TIDE score in the HRisk and LRisk group. (D–K) Sensitivity analysis of key drugs in HRisk and LRisk groups. * represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001, ns, no significance in comparison with the control group.

Drug Sensitivity Analysis

In addition to immunotherapy, targeted therapy and chemotherapy are still important treatments for HCC.²⁷ Therefore, we used the GDSC database to screen for drugs that are sensitive to the HRisk group. (Figure 7D–K) Our results showed that patients in the high-risk group were more sensitive to Sorafenib, Pazopanib, ABT-888, AZD-2281, Bleomycin, Doxorubicin, Cytarabine, and Gemcitabine, providing new insights for personalized treatment of patients with HCC. These chemotherapy drugs have therapeutic value for the high-risk group when immunotherapy is ineffective.

Expression Levels of PCRGs

Violin plots illustrate the expression of the model genes across different cell populations. (Figure 8A–H and Figure S3) To further investigate the mRNA expression levels of the model-related genes, we analyzed 10 paired fresh tissue samples of normal and HCC origin previously collected by our group. The qPCR results indicated that, except for CD79A, all other genes were significantly upregulated in tumor tissues (Figure S4). Interestingly, we found that these genes were also differentially expressed in hepatic cells. Consequently, we employed qRT-PCR to confirm the differential expression of model genes between HCC cell lines and normal liver cell lines. Notably, SSR3 was significantly upregulated in the HCC cell lines (Figure 8I–O). To further investigate the role of SSR3 in HCC, we selected SSR3 for subsequent experiments.

Figure 8 Verifying the expression of genes that constitute the risk model in single-cell datasets and HCC cell lines. (A–H) The violin plots show the major differentially expressed genes in different cell types. (I–O) Relative expression levels of HepaG2, Huh-7, Hepa3B, MHCC97H, MHCC97L, SUN449, and HCCLM3 in normal cell lines (LO2) and HCC cell lines via qRT-PCR. For bar plots, data are presented as means ± SD. One-way ANOVA with multiple comparisons was used to compare the mean of the control column with the mean of every other column.* represents p < 0.05, ** represents p < 0.01,*** represents p < 0.001, **** represents p < 0.0001, ns, no significance in comparison with the control group.

Knockdown of SSR3 Improves HCC Cell Proliferation and Migration in Vitro

We investigated the biological function of SSR3 in the proliferation and migration of MHCC97H and HCCLM3 cell lines. First, SSR3 was knocked down in Huh-7 and HCCLM3 cell lines(Figure 9A and B). Cell proliferation was assessed using the CCK8 assays (Figure 9C and D), and colony formation assays (Figure 9E). We observed that SSR3 knockdown reduced cell proliferation. Meanwhile, transwell assays (Figure 9F) revealed that the knockdown of SSR3 markedly weakened migration ability. In summary, our results indicated that SSR3 knockdown significantly impaired HCC progression.

Figure 9 CCK8 assay, colony formation assay, and transwell migration assay revealed that SSR3 promotes proliferation and migration of HCC. (A and B)) RT-qPCR validation confirmed that SSR3 was knocked down in Huh-7 and HCCLM3 cell lines. (C and D) CCK8 assay revealed that knocking down SSR3 could inhibit the proliferation of HCC cells. (E) Colony formation assay revealed that knocking down SSR3 could inhibit the colony formation ability of HCC cells. (F) Transwell migration assay revealed that knocking down SSR3 could inhibit the migration ability of HCC cells. For curves and bar plots, data are presented as means ± SD. One-way ANOVA with multiple comparisons was used to compare the mean of the control column with the mean of every other column. ** represents p < 0.01,*** represents p < 0.001, **** represents p < 0.0001, ns, no significance in comparison with the control group.

Discussion

TIME plays a crucial role in the development of HCC by modulating the function and metabolism of immune cells, thereby influencing the overall dynamics of tumor progression.²⁸ In recent years, the roles of various immune cells, such as T cells, macrophages, and neutrophils, have been extensively explored in HCC, highlighting their potential as innovative therapeutic targets.^9–11 However, the function of plasma cells remains insufficiently understood. Recently, scRNA-seq has been regarded as one of the best tools to interpret TIME heterogeneity. In our study, we discovered that plasma cells made the greatest contribution to HCC development through scRNA-seq analysis. Previous studies have reported that plasma cells can regulate TIME by secreting antibodies or cytokines, ultimately affecting tumor development.^12,13 Therefore, we constructed a prognostic model based on PCRGs in HCC to guide clinical decision-making and personalized treatment.

In this study, eight PCRGs were identified as prognostic biomarkers for HCC through comprehensive bioinformatics analysis. These genes included ST6GALNAC4, SEC61A1, SSR3, RPN2, PRDX4, TRAM1, SPCS2, and CD79A. ST6GALNAC4 is linked to poor prognosis in various tumors, such as HCC and uterine corpus endometrial carcinoma.^29,30 It has been reported that ST6GALNAC4 promotes proliferation and invasion in HCC by driving high expression levels of TGFBR2.²⁹ SEC61A1, the main subunit of the SEC61 complex, is located in the endoplasmic reticulum membrane and is associated with the development of tumors, including chronic obstructive pulmonary disease (COPD), liver cancer, and colorectal cancer.^31–33 SSR3, an ER-associated protein, plays a crucial role in protein secretion by mediating the translocation of nascent proteins, and is closely related to many tumor types, including HCC.³⁴ RPN2, a crucial component of the oligosaccharyltransferase complex, is responsible for N-glycosylation of many proteins and has been found to promote malignant progression in glioma, bladder cancer, and colon cancer.³⁵ PRDX4, part of a family of enzymes including PRDX1, PRDX2, PRDX3, PRDX4, PRDX5, and PRDX6, plays significant roles in HCC, breast cancer, and prostate cancer. TRAM1, also known as SRC-3, i is involved in ER-associated degradation of substrate stability. Previous studies have demonstrated that knocking down TRAM1 inhibits tumor growth in pancreatic ductal adenocarcinoma.³⁶ SPCS2, a subunit of the signal peptidase complex, is mainly involved in the removal of signal peptides from newly synthesized proteins in the endoplasmic reticulum, and its expression levels in pediatric acute myeloid leukemia may have diagnostic and prognostic value.³⁷ CD79A, primarily expressed in B cells, forms a heterodimeric B cell receptor (BCR) with CD79B, participating in B cell signal transduction and development.³⁸ All these genes is associated with patient survival in HCC. These findings suggest that our PCRG model is strongly associated with the prognosis of patients with HCC.

To further investigate the role of the model gene in HCC cell lines, we selected SSR3 for subsequent experiments. We found that SSR3 knockdown significantly inhibited the proliferation and migration of HCC cells in vitro. Bioinformatics analysis revealed that high expression of SSR3 is insensitive to immunotherapy. These findings suggest that our PCRGs model is strongly associated with the prognosis of patients with HCC. Additionally, the PCRGs model accurately and stably predicted the prognosis of patients with HCC in the training, internal validation, and external validation cohorts. We evaluated the stability and accuracy of the model in predicting overall survival in HCC using calibration curves and ROC analysis. Overall, our PCRG models can serve as reliable biological indicators for predicting outcomes in patients with HCC.

Another important finding of our study is the analysis of the mechanisms underlying the PCRG model. GSVA analysis indicated that the HRisk group was mainly associated with the activation of mTORC1 and PI3K/Akt/mTOR signaling pathways. The mTORC1 signaling pathway plays a crucial role in tumor development by regulating cell growth, proliferation, and survival.³⁹ Similarly, the PI3K/Akt/mTOR signaling pathway is pivotal in HCC by promoting cell proliferation, inhibiting apoptosis, and regulating metabolism.⁴⁰ In the future, drugs targeting the PI3K/Akt/mTOR signaling pathway may provide new hope for HCC treatment. To further explore the upstream regulatory mechanisms of this prognostic signature, we examined the miRNAs upstream of the eight hub genes in the PCRG model. We identified significant mRNA–miRNA pairs by predicting the targeting relationships between candidate genes and miRNAs, thereby enhancing our understanding of HCC tumorigenesis and progression.

Recently, immunotherapy has become one of the main treatments for advanced HCC.⁴¹ Previous studies have reported that tumor sensitivity to immune checkpoint therapy can be predicted using the TIDE score.²⁶ We observed that an HRisk score was significantly associated with higher immune dysfunction scores, indicating an immune dysfunction status. Additionally, we explored the correlation between risk score and immune cell infiltration in HCC. The results revealed a negative correlation between the risk score and infiltration of four types of immune cells (B cells, M1 macrophages, resting NK cells, and memory CD4 T cells) and a positive correlation with Treg cells. Treg cells are known to suppress anti-tumor immunity and facilitate tumor progression,⁴² which may explain the poor therapeutic efficacy of immunotherapy in the HRisk group. Moreover, we investigated drug sensitivity in the HRisk group by using the GDSC database. Our results showed that patients in the HRisk group were more sensitive to Sorafenib, Pazopanib, ABT-888, AZD-2281, Bleomycin, Doxorubicin, Cytarabine, and Gemcitabine, providing new insights for personalized treatment of patients with HCC. Sorafenib is a first-line treatment for advanced HCC.⁴¹ ABT-888 (veliparib) and AZD-2281 (olaparib) are PARP inhibitors that enhance the effectiveness of cancer treatments by preventing DNA repair in tumor cells, leading to their increased susceptibility to cell death.⁴³ Overall, our findings indicate that although the HRisk group is less responsive to immunotherapy, effective options for targeted therapy and chemotherapy remain.

However, this study had several limitations. Firstly, the analyzed data were acquired from public databases, necessitating external validation. Prospective clinical studies with a broader sample range are crucial for further validation. If validated, our PCRGs model could be a valuable tool to aid clinicians in personalized risk assessment and treatment decision-making. Moreover, the data set in our study was filtered to 5095 cells, which somewhat improved the quality and reliability of the data but was limited by the small number of samples. In future studies, the analysis could be further optimized by increasing the sample size and exploring other approaches. This would allow for a more comprehensive understanding of the various cell types and their functions, thereby enhancing the depth and robustness of the findings.In conclusion, these findings suggest that the PCRGs prognostic model can accurately and reliably predict the prognosis of HCC patients. Additionally, our findings suggest that inhibiting the expression of SSR3 can suppress the proliferation and migration capabilities of HCC cell lines.

Conclusions

In this study, we identified eight PCRGs that predict the prognosis of patients with HCC. We also quantified the expression of these model genes in HCC cell lines. Additionally, our study demonstrated that SSR3 knockdown inhibited the proliferation and migration of HCC cells. In summary, we confirmed the close association between the PCRGs and HCC prognosis through integrative bioinformatics analysis and functional experiments, providing new insights to guide clinical decision-making and personalized treatment.