Hepatocellular carcinoma: An analysis of the expression status of stress granules and their prognostic value

Abstract

BACKGROUND

Hepatocellular carcinoma (HCC) is a global popular malignant tumor, which is difficult to cure, and the current treatment is limited.

AIM

To analyze the impacts of stress granule (SG) genes on overall survival (OS), survival time, and prognosis in HCC.

METHODS

The combined The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC), GSE25097, and GSE36376 datasets were utilized to obtain genetic and clinical information. Optimal hub gene numbers and corresponding coefficients were determined using the Least absolute shrinkage and selection operator model approach, and genes for constructing risk scores and corresponding correlation coefficients were calculated according to multivariate Cox regression, respectively. The prognostic model’s receiver operating characteristic (ROC) curve was produced and plotted utilizing the time ROC software package. Nomogram models were constructed to predict the outcomes at 1, 3, and 5-year OS prognostications with good prediction accuracy.

RESULTS

We identified seven SG genes (DDX1, DKC1, BICC1, HNRNPUL1, CNOT6, DYRK3, CCDC124) having a prognostic significance and developed a risk score model. The findings of Kaplan-Meier analysis indicated that the group with a high risk exhibited significantly reduced OS in comparison with those of the low-risk group (P < 0.001). The nomogram model’s findings indicate a significant enhancement in the accuracy of OS prediction for individuals with HCC in the TCGA-HCC cohort. Gene Ontology and Gene Set Enrichment Analysis suggested that these SGs might be involved in the cell cycle, RNA editing, and other biological processes.

CONCLUSION

Based on the impact of SG genes on HCC prognosis, in the future, it will be used as a biomarker as well as a unique therapeutic target for the identification and treatment of HCC.

Key Words: Stress granule genes, Hepatocellular carcinoma, Gastrointestinal neoplasms, Bioinformatics prognosis, Prognostic value

Core Tip: Hepatocellular carcinoma (HCC) is challenging to cure because the available treatments are ineffective, and this cancer demonstrates a significant degree of resistance to standard chemotherapy and radiotherapy modalities. In cancer cells, the formation of stress granules (SGs) is thought to protect them from apoptosis and induce resistance to anti-cancer drugs or radiotherapy, making SGs a potential target for cancer therapy. Throughout this investigation, we found that signatures of seven SG genes can serve as independent predictors of overall survival in patients with HCC, and survival and prognosis analyses have been conducted for each expression of the hub gene.

INTRODUCTION

Liver cancer ranks as the fifth highest prevalent malignancy globally and the fourth primary contributor to cancer-associated mortality[1], with a calculated number of new cases estimated to be 906000, corresponding to 830000 deaths worldwide each year[2]. There are two primary forms of significant liver disorders, cancer-hepatocellular carcinoma (HCC) and intrahepatic cholangiocarcinoma, and less prevalent tumors, including hemangiosarcoma and hepatoblastoma, of which HCC constitutes over 80% of the global incidence of primary liver cancer[3]. Asian countries exhibit the highest prevalence of HCC, which ranks as the second highest prominent cause of death related to cancer, and it accounts for 72.5% of global incidence and 73.3% of global mortality in 2020[4]. HCC can occur in a category of chronic liver disorders, for instance, and non-alcoholic fatty liver disease, alcoholic liver disease, chronic hepatitis B and C virus infection, most of which begin with an abnormal buildup of fat through liver cells, a status known as steatosis and hence fatty liver disease. Over a period of time, hepatocyte death and inflammation are caused by endoplasmic reticulum stress, lipotoxicity, and mitochondrial dysfunction (steatohepatitis) and result in a buildup of fibrotic tissue throughout the liver cells. The progression of unresolved chronic inflammation and fibrosis is observed over a certain amount of time, and regenerative nodules comprising hepatocytes with low differentiation, arise within the parenchyma, ultimately resulting in cirrhosis[1,5]. Most HCC (80%-90%) is preceded by cirrhosis[6]. HCC is challenging to cure because the available treatments are ineffective, and this cancer demonstrates a significant degree of resistance to standard chemotherapy and radiotherapy modalities. The most efficacious approach continues to be surgical resection or liver transplantation. However, not all patients meet the criteria to undergo surgical intervention, with recurrence rates approaching 70% five years after resection[1,7]. Therefore, further understanding of the HCC molecular mechanisms is needed to construct new effective prevention and treatment approaches.

The pathogenesis of HCC is complex, involving multiple molecular malfunctions, including cell cycle dysregulation, chromosomal instability, DNA methylation changes, immune regulation, epithelial-mesenchymal transition, HCC stem cells progression, and the dysregulation of microRNA (miRNA)[8], among these, aberrant gene expression post-transcriptional regulation is the most crucial element that drives the progression of liver carcinogenesis. Alterations in the expression or function of miRNAs and RNA binding proteins (RBPs), or both, can lead to the reduction of mRNA levels or inhibit the translation of transcripts that regulate liver metabolism, inflammation, and carcinogenesis[9]. In the case of gene mutations, alterations in post-transcriptional mechanisms that regulate gene expression have the potential to result in the oncogenes’ overexpression or the tumor suppressors’ under-expression, thus, promoting the progression of cancer. Specific regulatory processes occur within small cytoplasmic ribonucleoprotein (RNP) foci, including stress granules (SGs) or p-bodies, and the mRNA’s destiny is identified (i.e., degradation and translation)[10]. SGs are RNP granules located in the cytoplasm that emerge as a result of stress. These particles are thought to be evolutionarily conserved cellular defense mechanisms that can resist a variety of stresses[11]. During times of physiological stress, SG selectively isolates specific transcripts and proteins from the soluble fraction of cytoplasm, which may be possible to be a mechanism that serves as a response to stress, regulating both gene expression and cell signaling[12]. SGs are dynamically regulated according to cell type and stress, they can influence cell development, apoptosis, and senescence, and their regulation is frequently linked to age-related human disorders, including cancer and neurodegenerative diseases[13]. The overexpression of various components of SGs has been detected in different tumor types, including sarcoma, pancreatic cancer, HCC, and glioblastoma[14]. Given that cancer cells are naturally subjected to significant stress during both tumorigenesis and anti-cancer therapy, there is growing evidence to suggest that modifications in the formation of SGs may serve as a protective mechanism for cancer cells against apoptosis[15]. Therefore, further investigation is required to determine the prognostic significance of SGs formation in HCC.

Throughout this investigation, a comprehensive analysis of SGs was performed in HCC. The gene set associated with SGs was obtained from the RNA granule database (rnagranuledb.lunenfeld.ca) and the research conducted by Youn et al[16]. Sequencing and relevant clinical data were obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases for a cohort of 340 patients diagnosed with HCC. The objective was to develop prognostic models. The GEO database was utilized to obtain the datasets that were used for validating the expression of SG genes, specifically, GSE25097 and GSE36376. According to the aforementioned data sources, a comprehensive bioinformatics analysis was conducted. Using the SG genes in HCC, a risk scoring system HCC was developed and subsequently validated using the TCGA dataset. Furthermore, functional analysis and Gene Set Enrichment Analysis (GSEA) of SG genes investigated these genes’ prospective functions and pathways within HCC. These investigation outcomes show that signatures of seven SG genes can serve as independent predictors of overall survival (OS) in patients with HCC, and survival and prognosis analyses have been conducted for each expression of the hub gene.

MATERIALS AND METHODS

Data collection and pre-processing

This analysis was performed by downloading liver cancer data sets GSE25097 and GSE36376 from GEO official website (https://www.ncbi.nlm.nih.gov/geo/) database[17-26], and the datasets utilized in the study were derived from samples of Homo sapiens. The study employed the GPL10687 Rosetta/Merck Human RSTA Affymetrix 1.0 microarray and Custom CDF to analyze the GSE25097 dataset. This dataset consisted of 289 normal samples and 268 cases of HCC. The dataset GSE36376 is based on the platform GPL10558 Illumina HumanHT-12 V4.0 expression beadchip, and 0 normal samples and 433 HCC samples are included in the analysis.

Utilization of the UCSC Xena platform (https://xenabrowser.net/datapages/) allowed for the acquisition of TCGA dataset for this investigation. We downloaded TPM expression profile data, count data matrix, survival data, and clinical data of bulk-seq transcriptome sequencing of cancer tissues (371 cases) and neighboring tissues (50 cases) of HCC patients from UCSC Xena[27]. UCSC Xena has performed log2 normalization on its stored data to eliminate the dimension of gene expression data within samples. We used the annotation file (gencode.v22.annotation.gene.probeMap) provided by UCSC Xena to perform Gene Symbol re-annotation on the expression data to conduct an upcoming analysis. Download the clinical data of respective individuals through transcriptome sequencing, including gender, age, pathological condition, and survival prognosis. To achieve diverse analytical objectives, we decided to exclude patients with incomplete prognostic data in the context of endometrial cancer. We then integrated this data with GSE25097 and GSE36376 to obtain sequencing data and pertinent clinical information for 340 individuals with liver malignancy (Table 1). This was done to construct a prognosis model. The gene set pertaining to SGs utilized in this investigation was obtained from the RNA granule database (rnagranuledb.lunenfeld.ca) as well as from the study conducted by Youn et al[16]. SG genes with relevance score > 1 were extracted, and 241 genes were finally identified.

Table 1 Baseline clinical characteristics.

ID	Dataset	Platforms	Sample	Source types	Ref.
1	TCGA-LIHC	-	50 HC, 374 HCC	Liver tissue	-
2	GSE25097	GPL10687	289 HC, 268 HCC	Liver tissue	[17-22]
3	GSE36376	GPL10588	433 HCC	Liver tissue	[23-25]

TCGA-LIHC: The Cancer Genome Atlas-Liver Hepatocellular Carcinoma; HCC: Hepatocellular carcinoma; HC: Healthy control.

Least absolute shrinkage and selection operator-Cox regression analysis to screen prognostic genes

To assess the prognostic capacity of gene expression for normal and HCC groups, we extracted SG-related gene expression using a limma package based on the combined data set. Subsequently, the study employed univariate lasso, multivariate Cox, and Cox regression analyses to screen the prognostic genes and construct the prognostic model. Initially, the study used a univariate Cox proportional regression analysis to investigate the connection between differential gene expression and OS and retained genes with a P value of less than 0.05 (or 0.1). After that, The Least absolute shrinkage and selection operator (LASSO) algorithm was employed in the univariate Cox regression analysis in order to eliminate multicollinearity and determine which variables were significant. Multivariate Cox regression analysis was used to obtain more accurate independent prognostic factors, specifically hub genes. The final screening was conducted through stepwise regression. Ultimately, the algorithm for calculating risk scores was developed by integrating optimized gene expression data with the corresponding estimated Cox regression coefficients:

Where C_i exhibits the level of gene expression, G_i and is the coefficient related to the gene.

The specimens were categorized into two distinct groups, specifically the high- and low-risk groups, in accordance with the given risk score. The OS of the test set was analyzed through the utilization of a survival package, wherein both of log-rank test and Kaplan-Meier analysis were conducted. Survival and prognosis for each hub gene expression were analyzed. The study conducted an assessment of the reliability of survival prognostication through the application of time-dependent receiver operating characteristic (ROC) curves. The precision of prognosis or prediction was evaluated by computing the area under the curve (AUC) values using the time ROC R package. To explore differentially expressed genes (DEGs) across risk groups, we visualized differences in expression based on the Wilcoxon test and plotted boxplots.

Kyoto Encyclopedia of Genes and Genomes and Gene Ontology functional analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG)[28] database is a prevalent resource for the storage of data related to biological pathways, genomes, drugs, and diseases. KEGG[29] pathway enrichment of necroptosis differential genes was conducted utilizing R package clusterProfiler (version 4.2.0) to identify biological processes with significant enrichment, and the enrichment analysis was conducted and revealed a significance threshold of P value < 0.05. The present study performed a Gene Ontology (GO) enrichment analysis on the DEGs derived from the combined datasets of TCGA-Liver Hepatocellular Carcinoma (TCGA-LIHC), GSE25097, and GSE36376[30], which is a prevalent approach utilized in genes large-scale functional enrichment analyses across various categories and stages, the utilization of Fisher’s exact test for the computation of biological processes, respectively, (biological process), molecular function, and cellular component, the process of quantifying the differential genes contained within each GO entry was performed, and the study employed the hypergeometric distribution algorithm to assess the statistical significance level of the differential gene enrichment within every GO entry. The result of the calculation returns an enrichment significance P value, with lower values indicating more statistically significant. A P value less than 0.05 is reported as statistically significant. This article plays a part in visualizing it through bubble plots.

Gene Set Variation Analysis and GSEA enrichment assays

After screening DEGs among different risk groups, for exploring the pathway enrichment of differential genes, the present study conducted an enrichment analysis utilizing the Gene Set Variation Analysis (GSVA) approach on the set of (DEGs) in the combined datasets of TCGA-LIHC, GSE25097, and GSE36376 utilizing the GSVA package[31]. By downloading the gene set “c2.cp.kegg.v7.5.1.symbols.gmt” and setting the threshold as |log(fold change)| > 0.1 and P value < 0.05, we found pathways with significant enrichment differences among different subtypes. The GSVA technique is an unsupervised, non-parametric analytical approach. It is primarily employed to assess the outcomes of gene set enrichment analyses conducted on microarray and transcriptome data. The primary purpose of this technique is to assess the enrichment of distinct metabolic pathways across various samples. This is achieved by converting the gene expression matrix from one sample to another, resulting in a gene set expression matrix between samples.

Subsequently, GSEA has been conducted on all genes in the combined datasets of TCGA-LIHC, GSE25097, and GSE36376 using the “clusterprofiler” R package. GSEA is a statistical method employed to estimate the distribution pattern of genes belonging to a pre-defined gene set within a gene table that is ranked based on the degree of association with a phenotype to judge their significance to phenotype[32]. The gene set is denoted as “c2.cp.v7.5.1.symbols.gmt”[33] was acquired from the MSigDB database (v7.5.1) and carried out GSEA analysis on all genes respectively, and set P value < 0.05 was reported as statistically significant.

Protein-protein interaction

The STRING database[34] serves as a platform for identifying potential interactions between established and projected proteins. Based on the combined datasets of TCGA-LIHC, GSE25097, and GSE36376, we established the networks of protein-protein interaction (PPI) for DEGs utilizing the STRING database. The STRING database (version 11.5) was utilized to upload the gene list, and the reliability was set to 0.4; PPI table data was subsequently purchased into Cytoscape (version 3.9.1) software for constructing PPI networks for visualization. Finally, the utilization of CytoHubba plug-in was employed to examine the hub genes that are present within the PPI network, and the top ten core genes were identified and chosen according to MCC score (Maximum clique centrality)[35]. Utilizing the Human Protein Atlas, database facilitated the visualization of core gene expression in various tissues.

Building a prediction nomogram

Based on the TCGA-LIHC, GSE25097, and GSE36376 combined datasets, we incorporated risk scores and clinicopathological characteristics into the model, followed by the application of multivariate Cox regression for assessing the risk scores impact and clinicopathological features on prognosis. The RMS package (version:5.1-4, http://cran.r-project.org/web/packages/rms/index.html) was then utilized for the production of nomograms and calibration plots that encompassed clinically relevant characteristics. The evaluation of calibration was conducted by means of a process that involved the mapping of anticipated probability derived from the nomogram to the actual probability of occurrence, considering the 45° line representing the best predictive value. The utilization of decision curve analysis (DCA) allowed for the evaluation of risk score clinical values in addition to the clinical pathological characteristics exhibited by the model with respect to the potential advantages for the patient. The larger the area under the curve, the more significant clinical benefit the patient would have.

Statistical method

The present investigation’s data processing and analysis were conducted using R software (Version 4.1.0). The continuous variables are presented as mean ± SD. The comparison of progression-free survival within the two groups was carried out through the utilization of Kaplan-Meier analysis in combination with the log-rank test. χ² or Fisher’s exact tests were utilized to assess and examine the statistically significant level between two sets of categorical variables. If not specified, the results were calculated and obtained by Spearman correlation analysis between the correlation coefficients of different molecules, and all statistical P values utilized in the analysis were two-sided, with P < 0.05, which was deemed statistically significant.

RESULTS

Project flow chart

In this study, the liver cancer data were downloaded from the GEO official website database and TCGA database for analysis. The project flow chart is shown in Figure 1.

Figure 1 Project flow chart. TCGA-LIHC: The Cancer Genome Atlas-Liver Hepatocellular Carcinoma; LASSO: Least absolute shrinkage and selection operator; PPI: Protein-protein interaction; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GSEA: Gene Set Enrichment Analysis; GSVA: Gene Set Variation Analysis; ROC: Receiver operating characteristic curve; DCA: Decision curve analysis.

Risk model construction

Following the extraction of SG gene expression in HCC, aiming to predict the clinical risk and survival outcome of samples through the difference of SG gene expression and function, based on the combined dataset of TCGA-LIHC, GSE25097, and GSE36376, we used LASSO model method to determine the optimal hub gene number (Figure 2A) and correspondence coefficient (Figure 2B), screened the best hub genes according to the ratio of 5:5, and divided the HCC samples into train group and test group. Finally, the genes for constructing risk scores were determined. The corresponding correlation coefficients were calculated based on multi-factor COX regression (Table 2) for calculating the risk score according to the hub genes expression level (Supplementary Table 1), and training specimens were separated into two different groups, which relied on the HCC samples median risk score (Figure 2C), and test group samples (Figure 2D) have been separated into two distinct groups according to their level of risk; namely, the high-risk group (high risk) and the low-risk group (low risk), and hub genes expression in train group specimens (Figure 2E) and test group samples (Figure 2F) within the high- and low-risk groups was visualized through the use of heat map. The red color indicates the upregulation of gene expression level within the sample, and the blue color represents the under-expression of gene expression level within the sample risk score = DDX1* 1.109698909 + DKC1* 0.624264345 + BICC1* - 0.489537897 + HNRNPUL1* - 2.330604694 + CNOT6* 1.986647774 + DYRK3* 0.584581122 + CCDC124* 0.635540329.

Figure 2 Risk grouping construction based on The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Least absolute shrinkage and selection operator is a gradient descent method for determining hub genes; B: Number of hub genes and their corresponding coefficients; C: Train group samples were classified into high- and low-risk groups based on the risk scores. Red indicates high-risk group samples, and blue indicates low-risk group samples; D: Test group samples have been classified into high-risk and low-risk groups based on the risk scores. Red indicates high-risk group samples, and blue indicates low-risk group samples; E: The hub genes expression within the high- and low-risk groups of the train group. Red indicates up-regulated expression, blue indicates under-expression of the expression, and white indicates zero expression; F: The hub genes expression within high- and low-risk groups of the test group. Red indicates up-regulated expression, blue indicates under-expression of the expression, and white indicates zero expression.

Table 2 Hubgenes screened by Least absolute shrinkage and selection operator with coefficients.

ID	coef	HR	HR (95L)	HR (95H)	P value
DDX1	1.109698909	3.033444913	1.220924802	7.536736107	0.01685573
DKC1	0.624264345	1.866872079	0.966659311	3.605418495	0.063029058
BICC1	-0.489537897	0.612909556	0.431027993	0.87153997	0.006421535
HNRNPUL1	-2.330604694	0.097236931	0.034257958	0.275994874	1.19E-05
CNOT6	1.986647774	7.291051498	2.746308357	19.35668725	6.67E-05
DYRK3	0.584581122	1.794239261	0.975801337	3.299129036	0.059953298
CCDC124	0.635540329	1.88804203	0.99423886	3.58535846	0.052100047

HR: Hazard ratio.

We then performed survival analysis on specimens from high- and low-risk groups within the combined dataset of TCGA-LIHC, GSE25097, and GSE36376. Survival analysis revealed that the samples of the low-risk group within the training group (Figure 3A) (P < 0.001) and the test group (Figure 4A) (P = 0.022) had longer survival times. Based on this, a survival analysis of all hub genes within the train group and test group was performed, respectively, and it was revealed that each hub gene expression within the train group and test group had significant differences in HCC survival. CCDC124 (P < 0.001) (Figure 3B), DDX1 (P < 0.001) (Figure 3C), DKC1 (P < 0.001) (Figure 3D), HNRNPUL1 (P = 0.043) (Figure 3E), CNOT6 (P < 0.001) (Figure 3F), BICC1 (P = 0.001) (Figure 3G), DYRK3 (P = 0.004) (Figure 3H) in train group; in the test group, CCDC124 (P < 0.001) (Figure 4B), DDX1 (P < 0.001) (Figure 4C), DKC1 (P < 0.001) (Figure 4D), HNRNPUL1 (P = 0.002) (Figure 4E), CNOT6 (P < 0.001) (Figure 4F), BICC1 (P = 0.002) (Figure 4G), DYRK3 (P = 0.006) (Figure 4H). Then the differential genes were screened among high and low-risk groups in accordance with the Wilcoxon rank sum test (Supplementary Table 2).

Figure 3 Survival analysis of high- and low-risk groups and hub genes in train group specimens according to The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Survival analysis of risk scores for train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; B: Survival analysis of CCDC124 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; C: Survival analysis of DDX1 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; D: Survival analysis of DKC1 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; E: Survival analysis of HNRNPUL1 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; F: Survival analysis of CNOT6 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; G: Survival analysis of BICC1 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples; H: Survival analysis of DYRK3 in train group samples. Red is indicative of high-risk group samples, and blue is indicative of low-risk group samples.

Figure 4 Hub genes survival analysis within the test group. A: Survival analysis of risk scores for test group samples. Red represents high-risk group samples, and blue indicates low-risk group samples; B: Survival analysis of CCDC124 samples in the test group. Red indicates high-risk group samples, and blue represents low-risk group samples; C: Survival analysis of DDX1 in test group samples. Red indicates high-risk group samples, and blue indicates low-risk group samples; D: Survival analysis of DKC1 in test group samples. Red represents high-risk group samples, and blue indicates low-risk group samples; E: Survival analysis of HNRNPUL1 in the test group. Red represents high-risk group samples, and blue indicates low-risk group samples; F: Survival analysis of CNOT6 in test group samples. Red indicates high-risk group samples, and blue indicates low-risk group samples; G: Survival analysis of BICC1 in test group samples. Red indicates high-risk group samples, and blue indicates low-risk group samples; H: Survival analysis of DYRK3 in test group samples. Red indicates high-risk group samples, and blue indicates low-risk group samples.

Functional analysis

The present study conducted functional typing of DEGs among the high-risk and low-risk groups within the combined TCGA-LIHC, GSE25097, and GSE36376 datasets. In order to determine biological processes that are significantly enriched, the differential gene expression analysis underwent GO and KEGG enrichment analysis via the clusterProfiler R package (Tables 3 and 4). Each term in GO enrichment results is sorted in ascending order according to P value, i.e., the more significant the enrichment, the higher. We selected the first 12 GO terms for circle diagram display (Figure 5A) and then selected the first 4 terms of the three functional categories to visualize with a bar chart (Figure 5B). It can be seen that in GO analysis, significantly different necroptosis genes (P < 0.05) have primarily participated in a collection of biological activities like the regulation and translation of the mRNA metabolic process, RNA catabolic process, and mRNA catabolic process (Figure 5D), and also include cell components such as cytoplasmic SG, cytoplasmic RNP granule, RNP granule and p-body (Figure 5E), and enriched molecular functions such as nucleic acid binding, translation regulator activity, translation factor activity, translation initiation factor activity, and RNA binding (Figure 5F). As shown in the findings of KEGG enrichment analysis, the genes whose expression levels were found to be differential have been enriched in the RNA degradation process, amyotrophic lateral sclerosis, spliceosome and nucleocytoplasmic transport pathways (Figure 5C), and the connection among gene expression and pathway enrichment was shown by string diagram for the first five pathways (Figure 5G).

Figure 5 Gene Ontology analysis and Kyoto Encyclopedia of Genes and Genomes analysis of differentially expressed genes between high- and low-risk groups in The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: The circle diagram visualization Gene Ontology (GO) enrichment analysis outcomes among the high- and low-risk groups, the outer circle is GO terms, red dots indicate genes overexpressed genes, blue dots indicate under-expressed genes, quadrilateral colors indicate the z-score of GO terms, blue indicates the z-score is negative, and the corresponding GO terms would probably exhibit inhibition, and red indicates z-score is positive, and it is more probable for the activation of the term to take place in the corresponding GO terms; B: Dot plot of GO enrichment analysis for differentially expressed genes; C: The abscissa represents the gene ratio. The ordinate represents the pathway name, the node size indicates the genes enriched number included within the pathway, and the node color is indicative of the P value; D: Biological processes function enrichment chord diagram, node color represents gene expression level. Red indicates the gene overexpressed gene, and blue indicates the under-expressed gene; E: Cellular component function enrichment chord diagram, node color represents gene expression level. Red indicates the gene overexpressed gene, and blue indicates the under-expressed gene; F: Molecular function enrichment chord diagram, node color represents gene expression level. Red indicates the gene overexpressed gene, and blue indicates the under-expressed gene; G: Kyoto Encyclopedia of Genes and Genomes function enrichment chord diagram, node color represents gene expression level. Red indicates the gene overexpressed gene, and blue indicates the under-expressed gene. FC: Fold change.

Table 3 Gene Ontology analysis results.

Ontology	ID	Description	Gene ratio	Bg ratio	P value	adj P value
BP	GO:0006417	Regulation of translation	46/93	463/18800	3.06e-49	2.93e-46
BP	GO:1903311	Regulation of mRNA metabolic process	31/93	294/18800	3.99e-33	1.91e-30
BP	GO:0006401	RNA catabolic process	29/93	280/18800	1.06e-30	3.39e-28
BP	GO:0006402	mRNA catabolic process	27/93	233/18800	6.75e-30	1.62e-27
BP	GO:0017148	Negative regulation of translation	27/93	248/18800	3.81e-29	7.3e-27
CC	GO:0035770	Ribonucleoprotein granule	44/95	261/19594	9.7e-58	1.16e-55
CC	GO:0036464	Cytoplasmic ribonucleoprotein granule	40/95	243/19594	1.09e-51	6.51e-50
CC	GO:0010494	Cytoplasmic stress granule	23/95	79/19594	1.3e-35	5.21e-34
CC	GO:0000932	P-body	16/95	93/19594	6.99e-21	2.1e-19
CC	GO:0016281	Eukaryotic translation initiation factor 4F complex	7/95	11/19594	1.63e-14	3.92e-13
MF	GO:0008135	Translation factor activity, RNA binding	16/95	84/18410	3.24e-21	3.23e-19
MF	GO:0003743	Translation initiation factor activity	14/95	51/18410	3.87e-21	3.23e-19
MF	GO:0090079	Translation regulator activity, nucleic acid binding	17/95	108/18410	5.48e-21	3.23e-19
MF	GO:0045182	Translation regulator activity	24/95	406/18410	4.16e-19	1.84e-17
MF	GO:0003730	mRNA 3’-UTR binding	18/95	222/18410	7.11e-17	2.52e-15

GO: Gene Ontology; BP: Biological process; CC: Cellular component; MF: Molecular function; UTR: Untranslated region.

Table 4 Kyoto Encyclopedia of Genes and Genomes analysis results.

Ontology	ID	Description	Gene ratio	Bg ratio	P value	adj P value
KEGG	hsa03018	RNA degradation	10/47	79/8164	1.55e-11	6.69e-10
KEGG	hsa03040	Spliceosome	7/47	147/8164	1.83e-05	0.0004
KEGG	hsa03013	Nucleocytoplasmic transport	6/47	108/8164	3.22e-05	0.0005
KEGG	hsa03015	mRNA surveillance pathway	5/47	97/8164	0.0002	0.0024
KEGG	hsa05014	Amyotrophic lateral sclerosis	7/47	364/8164	0.0045	0.0349
KEGG	hsa05416	Viral myocarditis	3/47	60/8164	0.0049	0.0349

KEGG: Kyoto Encyclopedia of Genes and Genomes.

We then performed GSVA and GSEA analyses using the GSVA package and clusterProfiler R package for all genes in the TCGA-LIHC, GSE25097, and GSE36376 combined datasets (Tables 5 and 6). GSVA analyses showed that the low-risk group exhibited primarily enriched results and up-regulated in vasopressin-regulated water reabsorption, cell cycle, progesterone-mediated oocyte maturation, basal transcription factors, and RNA degradation pathways compared to the high-risk group (Figure 6); It is enriched and down-regulated in maturity-onset diabetes of the Young, calcium signaling pathway, neuroactive ligand-receptor interaction, glycine serine and threonine metabolism and peroxisome proliferator-activated receptor (PPAR) signaling pathway. GSEA analysis suggests low risk versus high risk in fatty acid metabolism (Figure 7A), retinal metabolism (Figure 7B), drug metabolism cytochrome P450 (Figure 7C), PPAR signaling pathway (Figure 7D), steroid hormone biosynthesis (Figure 7E), threonine and glycine serine metabolism (Figure 7F), cell cycle (Figure 7G), butanoate metabolism (Figure 7H), complement and coagulation cascades (Figure 7I) and other pathways have significant enrichment differences (P < 0.05).

Figure 6 Gene Set Variation Analysis and Gene Set Enrichment Analysis of differential genes among high- and low-risk groups within The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. The main biological functions and pathway heat maps obtained from Gene Set Variation Analysis are displayed. Blue indicates low risk, and yellow indicates high risk. Red represents overexpression in this pathway, and blue indicates under-expression in this pathway. KEGG: Kyoto Encyclopedia of Genes and Genomes.

Figure 7 Gene Set Enrichment Analysis of differentially expressed genes between high- and low-risk groups in The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Gene Set Enrichment Analysis (GSEA) analysis illustrated that low risk was significantly up-regulated within the KEGG_FATTY_ACID_METABOLISM pathway compared to high risk; B: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_RETINOL_METABOLISM pathway compared to high risk; C: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_DRUG_METABOLISM_CYTOCHROME_P450 pathway compared to high risk; D: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_PPAR_SIGNALING_PATHWAY pathway compared to high risk; E: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_STEROID_HORMONE_BIOSYNTHESIS pathway compared to high risk; F: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM pathway compared to high risk; G: GSEA analysis illustrated that low risk was significantly down-regulated within the KEGG_CELL_CYCLE pathway compared to high risk; H: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_BUTANOATE_METABOLISM pathway compared to high risk; I: GSEA analysis illustrated that low risk was significantly up-regulated within the KEGG_COMPLEMENT_AND_COAGULATION_CASCADES pathway compared to high risk. KEGG: Kyoto Encyclopedia of Genes and Genomes.

Table 5 Gene Set Enrichment Analysis results.

ID	Description	Set size	Enrichment score	NES	P value	adj P value	Q values
KEGG_FATTY_ACID_METABOLISM	KEGG_FATTY_ACID_METABOLISM	34	0.71018864	3.228627657	1.00E-10	6.65E-09	4.53E-09
KEGG_RETINOL_METABOLISM	KEGG_RETINOL_METABOLISM	27	0.800774501	3.477938283	1.00E-10	6.65E-09	4.53E-09
KEGG_DRUG_METABOLISM_CYTOCHROME_P450	KEGG_DRUG_METABOLISM_CYTOCHROME_P450	35	0.696799368	3.236251677	2.13E-10	9.45E-09	6.43E-09
KEGG_PPAR_SIGNALING_PATHWAY	KEGG_PPAR_SIGNALING_PATHWAY	38	0.63438858	3.008127457	1.34E-09	4.44E-08	3.02E-08
KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM	KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM	20	0.794030134	3.009053132	7.28E-09	1.94E-07	1.32E-07
KEGG_STEROID_HORMONE_BIOSYNTHESIS	KEGG_STEROID_HORMONE_BIOSYNTHESIS	28	0.69286184	2.925063353	1.14E-08	2.52E-07	1.71E-07
KEGG_CELL_CYCLE	KEGG_CELL_CYCLE	63	-0.616156999	-1.947443161	4.85E-08	9.22E-07	6.27E-07
KEGG_BUTANOATE_METABOLISM	KEGG_BUTANOATE_METABOLISM	28	0.643803698	2.717953993	2.54E-07	4.22E-06	2.88E-06
KEGG_COMPLEMENT_AND_COAGULATION_CASCADES	KEGG_COMPLEMENT_AND_COAGULATION_CASCADES	47	0.53140243	2.615817223	3.37E-07	4.97E-06	3.39E-06
KEGG_TRYPTOPHAN_METABOLISM	KEGG_TRYPTOPHAN_METABOLISM	29	0.61199465	2.664359791	1.07E-06	1.43E-05	9.72E-06
KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450	KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450	35	0.551628438	2.562012164	3.01E-06	3.64E-05	2.48E-05
KEGG_BETA_ALANINE_METABOLISM	KEGG_BETA_ALANINE_METABOLISM	17	0.692581145	2.532317703	2.55E-05	0.000282589	0.000192344
KEGG_PEROXISOME	KEGG_PEROXISOME	38	0.464005438	2.200209055	6.76E-05	0.000691885	0.000470931
KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION	KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION	33	0.487548344	2.182626135	0.000121713	0.001156276	0.000787018
KEGG_ARACHIDONIC_ACID_METABOLISM	KEGG_ARACHIDONIC_ACID_METABOLISM	28	0.496563499	2.096348232	0.000261033	0.002314494	0.001575358
KEGG_PROPANOATE_METABOLISM	KEGG_PROPANOATE_METABOLISM	21	0.516414812	1.982973473	0.001656507	0.013769716	0.009372343
KEGG_TYROSINE_METABOLISM	KEGG_TYROSINE_METABOLISM	25	0.468326241	1.954109762	0.001866878	0.014605575	0.00994127
KEGG_ARGININE_AND_PROLINE_METABOLISM	KEGG_ARGININE_AND_PROLINE_METABOLISM	26	0.432080752	1.85908878	0.002731901	0.018388891	0.01251638
KEGG_SMALL_CELL_LUNG_CANCER	KEGG_SMALL_CELL_LUNG_CANCER	48	-0.516209717	-1.594193363	0.002759563	0.018388891	0.01251638
KEGG_OOCYTE_MEIOSIS	KEGG_OOCYTE_MEIOSIS	51	-0.503867568	-1.561948773	0.002765247	0.018388891	0.01251638

KEGG: Kyoto Encyclopedia of Genes and Genomes.

Table 6 Gene Set Variation Analysis results.

	logFC	AveExpr	t	P value	adj P value	B
KEGG_RNA_DEGRADATION	0.26385	-0.0063	9.67948	6.21E-20	9.44E-18	34.57459
KEGG_DRUG_METABOLISM_CYTOCHROME_P450	-0.2879	0.02526	-7.7606	8.06E-14	6.12E-12	20.81127
KEGG_PPAR_SIGNALING_PATHWAY	-0.2706	0.04111	-7.3521	1.23E-12	4.36E-11	18.15604
KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM	-0.3098	0.03838	-7.3347	1.38E-12	4.36E-11	18.04547
KEGG_CELL_CYCLE	0.21445	-0.0119	7.31519	1.56E-12	4.36E-11	17.92124
KEGG_STEROID_HORMONE_BIOSYNTHESIS	-0.258	0.03831	-7.2895	1.85E-12	4.36E-11	17.75812
KEGG_RETINOL_METABOLISM	-0.2937	0.04154	-7.277	2.01E-12	4.36E-11	17.67921
KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG	-0.2291	0.02587	-7.0965	6.41E-12	1.22E-10	16.54879
KEGG_SPLICEOSOME	0.23491	0.00117	6.77412	4.84E-11	8.17E-10	14.58431
KEGG_ARACHIDONIC_ACID_METABOLISM	-0.1958	0.0193	-6.5965	1.43E-10	2.03E-09	13.53201
KEGG_LINOLEIC_ACID_METABOLISM	-0.3017	0.04341	-6.5919	1.47E-10	2.03E-09	13.50512
KEGG_FATTY_ACID_METABOLISM	-0.278	0.02396	-6.5447	1.96E-10	2.48E-09	13.22943
KEGG_GLYCOLYSIS_GLUCONEOGENESIS	-0.1779	-0.007	-6.482	2.84E-10	3.33E-09	12.86608
KEGG_VASOPRESSIN_REGULATED_WATER_REABSORPTION	0.14008	-0.0067	6.40042	4.62E-10	5.01E-09	12.39671
KEGG_PROGESTERONE_MEDIATED_OOCYTE_MATURATION	0.14218	0.00097	6.26602	1.01E-09	9.71E-09	11.63418
KEGG_CALCIUM_SIGNALING_PATHWAY	-0.1189	0.00875	-6.2615	1.04E-09	9.71E-09	11.60868
KEGG_COMPLEMENT_AND_COAGULATION_CASCADES	-0.2481	0.03366	-6.2543	1.09E-09	9.71E-09	11.56852
KEGG_BASAL_TRANSCRIPTION_FACTORS	0.20569	-0.0133	6.20578	1.44E-09	1.18E-08	11.29665
KEGG_TYROSINE_METABOLISM	-0.2007	0.01699	-6.2017	1.47E-09	1.18E-08	11.27395
KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION	-0.1315	0.01542	-6.1916	1.56E-09	1.19E-08	11.21751

KEGG: Kyoto Encyclopedia of Genes and Genomes.

Table 7 Single variate Cox analysis of hub genes.

ID	HR	HR (95L)	HR (95H)	P value
DDX1	2.811787564	1.79184186	4.412303052	6.89E-06
DKC1	2.088261011	1.488787267	2.929118315	2.00E-05
BICC1	1.203898304	1.01529958	1.427530509	0.032791365
HNRNPUL1	1.717777659	1.128350776	2.615108838	0.011631909
CNOT6	2.223105295	1.547745966	3.193157832	1.53E-05
DYRK3	1.970782211	1.324333062	2.932783779	0.000822902
CCDC124	1.951959781	1.340432388	2.84247607	0.000486868

HR: Hazard ratio.

PPI analysis

A SG gene-based PPI network was established, utilizing differential expression data among the high-risk and low-risk groups within the combined dataset of TCGA-LIHC, GSE25097, and GSE36376 through the STRING website, and evaluated the obtained data and visualized the PPI network utilizing Cytoscape software, where red indicates overexpressed genes and blue indicates underexpressed genes (Figure 8A); subsequently, a set of ten genes were chosen as the core genes (EIF3D, EIF4G1, EIF4A2, EIF3B, EIF3A, EIF3E, EIF4B, EIF4E, DDX3X, and EIF3H) using cytoHubba plug-in analysis of MCC scores (Figure 8B).

Figure 8 Protein-protein interaction analysis. A: Protein-protein interaction analysis was conducted utilizing the STRING database for kinetokines-related genes, and the visualization of interaction relationship was carried out, where red indicates up-regulated genes and blue indicates downregulated genes; B: The resulting top 10 hub genes were analyzed utilizing the CytoHubba plug-in.

Hub gene correlation and prognosis analysis

For evaluating the influence of hub genes on HCC prognosis, we applied univariate Cox analysis to hub genes in TCGA-LIHC, GSE25097, and GSE36376 combined datasets and found that hub genes had a significant influence on HCC (P < 0.05) (Figure 9A). Concurrently, the calculation of hazard ratio value for each hub gene was performed to determine whether the gene was a protective or risk factor (Table 7). We then calculated and visualized the association among the hub genes (Figure 9B), and all seven of them showed positive correlations. The ROC curve for the prognostic model utilizing hub genes was plotted through the implementation of time ROC package for evaluating its performance and calculated the model’s AUC values at 1, 3, and 5 years within the training group (Figure 9C) and the test group (Figure 9D). The 1-, 3-, and 5-year survival rate AUC of train group were 0.886, 0.806, and 0.755, respectively (Figure 9C). The 1-, 3-, and 5-year survival AUCs in the Test group were 0.730, 0.694, and 0.683, respectively (Figure 9D).

Figure 9 Hub gene association and prognosis analysis relied on The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Univariate Cox analysis of Hub genes in hepatocellular carcinoma (HCC). The left side of the dotted line represents hazard ratio (HR) < 1, revealing that this gene expression serves as a protective factor for the prognosis of HCC; the right side of the dotted line represents HR > 1, suggesting that this gene expression is a predictive risk factor for HCC; B: Visualization of interactions between Hub-genes. The red line indicates a positive association between the two different genes; the size of dot indicates the P value of the gene’s effect on the prognosis of HCC samples, and the orange dot represents that the gene is a risk factor for HCC; C: Receiver operating characteristic (ROC) curve. Green represents 1-year survival, blue represents 3-year survival, and red represents 5-year survival; D: ROC curve of test group. Green represents 1-year survival, blue represents 3-year survival, and red represents 5-year survival. AUC: Area under the curve.

Prognostic model construction

A prognostic model was constructed by incorporating risk score, age, T stage, and N stage, and based on the combined dataset of TCGA-LIHC, GSE25097, and GSE36376, multivariate Cox regression revealed that risk score (P < 0.001), grade (P < 0.001), age (P = 0.022), and T stage (P < 0.001) were significantly correlated with HCC prognosis (Figure 10A). The risk stratification curve demonstrates that the prognostic model can well anticipate the survival time of HCC samples in different risk groups (Figure 10B), while the calibration curve shows that the overfitting control of the model is reasonable and the preparation is high (Figure 10C). We then visualized the prognostic model via nomograms, based on which 1-, 3-, and 5-year survival rates for HCC could be anticipated (Figure 11A). The DCA analysis showed that nomograms based on prognostic models were associated with a more significant clinical benefit in HCC patients at 1-year (Figure 11B), 3-year (Figure 11C), and 5-year (Figure 11D) survival.

Figure 10 Multivariate Cox analysis and model performance of prognostic factors based on The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Visualization of the correlation between risk scores and clinical characteristics and hepatocellular carcinoma (HCC) prognosis, ^aP ≤ 0.05, ^bP ≤ 0.01, ^cP ≤ 0.001; B: The prognostic model classified the HCC sample into high-risk and low-risk groups, and there was a significant variance in survival among the high- and low-risk groups. Red is indicative of the high-risk group, and blue is indicative of the low-risk group; C: Calibration curve of HCC clinical prognosis model. Green represents 1-year survival, blue represents 3-year survival, and red represents 5-year survival. OS: Overall survival.

Figure 11 Prognostic model nomograms and decision curve analysis curves based on The Cancer Genome Atlas-Liver Hepatocellular Carcinoma, GSE25097, and GSE36376 combined datasets. A: Clinical prognostic model nomogram of hepatocellular carcinoma (HCC). ^aP ≤ 0.05, ^cP ≤ 0.001; B: Decision curve analysis (DCA) curve for 1-year survival. The area under the curve represents the possible prognostic benefit for HCC patients based on this factor; C: DCA curve for 3-year survival. The area under the curve represents the possible prognostic benefit for HCC patients based on this factor; D: DCA curve for 5-year survival. The area under the curve represents the possible prognostic benefit for HCC patients based on this factor.

DISCUSSION

The prevalence of liver cancer has been observed to increase recently. HCC is the most commonly occurring primary liver cancer disorder. Since HCC does not respond to radiotherapy as well as chemotherapy, and there is no effective treatment for unresectable HCC, the mortality of HCC currently ranks second globally, and the 5-year survival rate is lower than 15%[36,37]. Recently, several studies have shown that SGs’ biogenesis and stability are related to cancer progression[14]. Cancer cells frequently encounter hypoxic, nutrient-deprived, and hyperosmotic environments to obtain their heightened metabolic requirements for proliferation[38]. Therefore, the tumor microenvironment is often full of various distresses, including high concentrations of reactive oxygen species and hypoglycemia, almost all of which may closely induce the formation of SGs[11]. SG is involved in tumorigenesis and cancer metastasis through tumor-related signaling pathways and other mechanisms. In addition, some chemotherapeutic agents were shown to promote SGs[39]. Hence, SGs’ role in cancer has garnered significant attention. Upregulation of SG components has been observed in HCC[39]. Nevertheless, additional research is required to discover the prognostic significance of SGs in HCC.

Throughout this investigation, the sequencing and relevant clinical data of individuals diagnosed with HCC were acquired from the combined datasets of TCGA-LIHC, GSE25097, and GSE36376. SG genes with relevance score > 1 were extracted from the RNA granule database (rnagranuledb.lunenfeld.ca), and the study of Youn et al[16], and a total of 241 genes were finally detected. Then the optimal number of hub genes was identified by the LASSO model, and the optimal hub genes were screened according to the ratio of 5:5. HCC samples were separated into train and test groups. Therefore, the median risk score is computed according to the hub genes expression, and the train samples and test samples are separated into high- (high risk) and low-risk groups (low risk). Survival analysis revealed that the low-risk group within the train and test groups had a higher OS rate and longer survival time. Between them, the7 hub genes expression levels (CCDC124, DDX1, DKC1, HNRNPUL1, CNOT6, BICC1, DYRK3) in the train group and test group had significant variations in survival time and OS of HCC (P < 0.05). Functional analysis has been conducted on these seven differential genes. In GO analysis, significant differential genes (P < 0.05) were primarily contributed within various biological processes, including modulation of mRNA metabolic, mRNA catabolic, RNA catabolic processes, and translation regulation. Also, they contained cell components that were rich in RNP, cytoplasmic RNP, cytoplasmic SGs, and p-body. According to the KEGG enrichment analysis, genes that exhibited differential expression were rich in the RNA degradation process, amyotrophic lateral sclerosis, spliceosome, and nucleocytoplasmic transport pathways. GSEA analysis revealed that low risk was significantly enriched in retinol metabolism, fatty acid metabolism, drug metabolism cytochrome P450, steroid hormone biosynthesis, PPAR signaling pathway, serine, glycine and threonine metabolism, cell cycle, butanoate metabolism, complement, and coagulation cascades compared to high risk (P < 0.05). We further analyzed the PPI of kinetogranin-related genes using the STRING database and revealed the top ten up-regulated core genes (EIF3D, EIF4G1, EIF4A2, EIF3B, EIF3A, EIF3E, EIF4B, EIF4E, DDX3X, EIF3H). Univariate cox analysis showed that seven hub genes (DDX1, DKC1, BICC1, HNRNPUL1, CNOT6, DYRK3, CCDC124) were risk factors for HCC prognosis. Through the utilization of risk score, multivariate Cox regression was identified (P < 0.001), grade (P < 0.001), age (P = 0.022), and T stage (P < 0.001) were significantly connected with HCC prognosis. Through risk stratification curve analysis, we confirmed that the prognostic model can well anticipate the HCC samples’ survival time within different risk groups. DCA analysis demonstrated that the prognostic model can make HCC individuals have higher clinical benefits in 1-, 3- and 5-year survival rates. Our findings suggest that the differential expression of SG-related genes is significantly connected to the OS rate, survival time, and prognosis of HCC.

SGs are induced by various stress conditions (heat, oxidative stress, and hypoxia). In cancer cells, the formation of SGs is thought to protect them from apoptosis and induce resistance to anti-cancer drugs or radiotherapy, making SGs a potential target for cancer therapy[40]. Most current understanding of stress particles is still in the molecular and cellular biological stages, and their practical impact on the human clinical environment is still insufficient. This study further confirmed that seven SG genes (DDX1, DKC1, BICC1, HNRNPUL1, CNOT6, DYRK3, and CCDC124) were significantly correlated with survival time of individuals with HCC and the OS rate and were also risk factors affecting the prognosis of HCC, providing a more theoretical basis for SGs as an innovative biomarker and prospective therapeutic target for the management of HCC. DDX1, DYRK3, and CCDC124 have been reported to be overexpressed in HCC[41-43], but DKC1, BICC1, HNRNPUL1, and CNOT6 are rarely reported in HCC.

DDX1 belongs to the DEAD-box RNA helicase family and contributes to various biological processes, including mRNA translation, tRNA splicing, rRNA processing, miRNA maturation, and DNA double-strand break (DSB) repair[44,45]. DDX1 has been found to promote the development of cancers, including neuroblastoma, retinoblastoma, testicular cancer, colorectal malignancy, breast malignancy, and liver cancer[41,46-48]. The study by Yuan et al[41] suggested that DDX1 is a poor prognostic factor of liver cancer, which is consistent with this study. In addition, the DDX1 promoter methylation level within the HCC tissues exhibited reduced levels than in normal tissues. The positive association between DDX1 expression and infiltration of immune cells (including dendritic cells, macrophages, B cells, and T cells) is primarily involved in the cell cycle and RNA editing[41].

Kan et al[49] showed that the DKC1 gene performs a vital function in the progression of cells in colorectal cancer, and a higher DKC1 expression level induces colorectal cancer progression by increasing ribosomal protein expression in a pseudouridine synthase activity-dependent manner. Current studies have shown that BICC1 can regulate biological processes, including apoptosis and proliferation, and has a poor association, prognosis for oral cancer and gastric cancer, and immune infiltration[50-52].

HNRNPUL stimulates the end excision of DNA, induces DNA DSB repair via homologous recombination, and participates in DNA damage response. It has been discovered that higher HNRNPUL1 expression levels are associated with an increased likelihood of future occurrences and are an independent risk factor for disease-free survival in esophageal squamous cell carcinoma individuals administrating platinum-based chemotherapy following the surgery[53].

CNOT6 has been considered to be up-regulated within acute lymphoblastic leukemia, acute myeloid leukemia, androgen-independent prostate cancer, and lung cancer cells, and expression of CNOT6 affects DNA mismatch repair[54]. The expression of DYRK1A and DYRK3 exhibits a correlation with immune-infiltrating cells within the tumor microenvironment and is overexpressed in MSI subtypes, which can serve as biomarkers for immunotherapy[42].

CCDC124 is a cell division-associated RBP significantly up-regulated in esophageal, adrenal, endometrial, liver, ovarian, thyroid, and bladder cancers[43]. In the present investigation, differential expression of the seven SG genes in HCC is found to affect the OS rate, survival time, and prognosis of HCC patients. The SG becomes a biomarker for predicting the stage of progression and the therapeutic effect of tumors, and the targeted SG can serve as a new method for tumor treatment. However, additional research is necessary to understand the specific mechanism fully.

Furthermore, the differential expression of SG-related genes in HCC was found, which was related to RNA degradation, nucleic acid transport, and splicing pathways, as well as fatty acid, retinol, cytochrome P450, PPAR, sterol synthesis, amino acid metabolism, and other pathways. In GO analysis, the significantly different necroptosis genes primarily participated in mRNA metabolic process regulation, translation regulation, RNA catabolism process, mRNA catabolism process, and other different biological processes, as well as cellular components such as RNP particles, cytoplasmic RNP particles, cytoplasmic stress particles, and p-bodies. Moreover, the findings of KEGG enrichment analysis indicate that the genes exhibiting differential expression are significantly enriched in the process of RNA degradation, amyotrophic lateral sclerosis, spliceosome, and nucleoplasmic transport pathways. Further GSVA analysis showed that according to the low-risk group, the high-risk group primarily enriched and up-regulated the vasopressin-regulated water reabsorption, cell cycle, progesterone-mediated oocyte maturation, basal transcription factors, and RNA degradation pathways. It is enriched and downregulated in young-onset adult-onset diabetes, calcium signaling pathway, neuroactive ligand-receptor interaction, glycine, serine and threonine metabolism, and PPAR signaling pathway. GSEA analysis revealed that there were significant enrichment variations between the high- and low-risk groups throughout drug metabolism, retinal metabolism, fatty acid metabolism, cytochrome P450, PPAR signaling pathway, steroid hormone biosynthesis, serine, glycine and threonine metabolism, cell cycle, butyrate metabolism, complement, and coagulation cascade.

Although our present investigation offers novel perspectives on the correlation between stress particles and prognostic significance in patients with HCC, it is crucial to acknowledge the limitations that must be taken into account. At first, this investigation involved only one dataset, which may have led to selection bias. Second, due to the use of online databases as the primary source of data in this study, certain critical clinical data, including specific chemotherapy treatments administered to the patients, could not be obtained. Third, additional research and thorough in vivo and in vitro experimental verification are crucial to clarify the biological roles and fundamental mechanisms of SGs in HCC.

CONCLUSION

In conclusion, from the TCGA and GSE datasets, we found that seven SG genes (DDX1, DKC1, BICC1, HNRNPUL1, CNOT6, DYRK3, and CCDC124) exhibited differential expression within HCC. GO and GSEA analysis revealed that these SGs may be involved in RNA degradation, nucleic acid transport, splicing pathways, as well as fatty acid, retinol, cytochrome P450, PPAR, sterol synthesis, and amino acid metabolism. Nomogram models with good prognostic accuracy were constructed for 1-, 3-, and 5-year OS predictions. These may contribute novel ideas to the study of HCC. Differential expression of the seven SG genes in HCC is found to affect the OS rate, survival time, and prognosis of HCC patients. The SG becomes a biomarker for predicting the stage of progression and the therapeutic effect of tumors, and the targeted SG can serve as a new method for tumor treatment. Although, the data mentioned above were obtained through bioinformatics analysis and required confirmation from a significant number of prospective study samples.

ACKNOWLEDGEMENTS

We acknowledge The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) database for providing their platforms and contributors for uploading their meaningful datasets. TCGA and GEO belong to public databases. The patients involved in the database have obtained ethical approval. Users can download relevant data for free for research and publish relevant articles. Our study is based on an open source data, so there are no ethical issues and other conflicts of interest.