Identification of pyroptosis subtypes and prognosis model of hepatocellular carcinoma based on pyroptosis‐related genes

Abstract

Background

Hepatocellular carcinoma (HCC) is a common malignant tumor with poor prognosis. Pyroptosis, a type of programmed cell death, regulates tumor cell development. However, the role of pyroptosis‐related genes (PRGs) in HCC and their association with prognosis are unclear.

Methods

We conducted bioinformatics analysis to identify PRGs in The Cancer Genome Atlas‐Liver Hepatocellular Carcinoma (TCGA‐LIHC) patients. Consensus clustering classified patients into different subtypes. We used LASSO regression to established a pyroptosis subtype‐related score (PSRS) related to prognosis. OncoPredict identified potential pharmaceuticals based on PSRS.

Results

We found 20 HCC‐related PRGs in 335 TCGA‐LIHC patients. Consensus clustering classified patients into two subtypes. Subtype I had better overall survival and higher response to anti‐PD1 treatment. The prognostic model involving 20 genes predicted poorer prognosis for high‐PSRS group. The model was validated in two external cohorts. OncoPredict identified 65 potential pharmaceuticals based on PSRS.

Conclusion

Our investigation revealed a correlation between pyroptosis and HCC. We established PSRS as independent risk factors for predicting prognosis. The study paves the way for using PRGs as prognostic biomarkers and exploring personalized therapy for HCC.

INTRODUCTION

Hepatocellular carcinoma (HCC) is the primary form of liver cancer, the sixth most common malignancies and the fourth leading cause of cancer-related deaths worldwide.¹ HCC typically develops after prolonged periods of chronic hepatitis, liver fibrosis, and cirrhosis.² HCC has extreme heterogeneity and poor prognosis,³ with a 5-year survival rate of only 18%.^4,5 The current prognosis of HCC relies mainly on clinicopathological staging frameworks, yet several marker characteristics, including TP53 gene mutations and genes associated with cellular proliferation, have been identified as survival predictors.^6–8 Therefore, there is an urgent need for novel prognostic biomarkers to predict survival and design personalized treatment strategies for patients with HCC.

Pyroptosis is a type of cell death mediated by gasdermins and inflammasomes, which occurs in vertebrates as an innate immune response mechanism.^9,10 Dysregulation of pyroptosis may lead to impaired pathogen clearance efficacy, dysfunctional adaptive immune response, and tissue damage. Evidence suggests that pyroptosis is involved in the onset and progression of various human diseases, including malignant tumors. Recent studies have shown that chemically-induced pyroptosis of tumor cells can occur in vitro and in vivo.¹¹ As pyroptosis may potentially promote cancer progression by triggering the release of multiple inflammatory mediators, including IL-1 and IL-18, some researchers believe that it is a protumorigenic mechanism.^12–14 Nevertheless, recent researches have validated that externally induced pyroptosis exhibits robust antitumor activity.^15,16

In this study, we investigated the mRNA expression of pyroptosis-related genes (PRGs) in 355 patients with HCC from The Cancer Genome Atlas (TCGA)-Liver Hepatocellular Carcinoma (LIHC) database and classified them into distinct pyroptosis subtypes. We further elucidated the relationship between pyroptosis subtypes and prognosis, and developed a prognostic model comprising 20 pyroptosis subtype-related genes. The model's accuracy was verified in two external HCC cohorts, demonstrating its potential for clinical decision-making in HCC patient management. Furthermore, we used integrated analysis to scrutinize differences in genomic variations, tumor microenvironment, and immunogenomic patterns between two pyroptosis subtypes. We also identified a novel group of potential drug molecules for patients with HCC based on PSRS, which may provide more therapeutic options. Our study seeks to personalize survival prediction and treatment options for patients with HCC through novel pyroptosis-based molecular classification and prognostic model.

MATERIALS AND METHODS

Acquisition of multiomic data of patients with HCC

The current study utilized multiomic data from three HCC cohorts, namely TCGA-LIHC,¹⁷ LIRI-JP,¹⁸ and CHCC-HBV.¹⁹ All transcriptome data were obtained through RNA sequencing, while genomic variations data were obtained through whole genome sequencing. TCGA-LIHC cohort data were acquired from the Genomic Data Commons (GDC) website (, accessed on 12 December 2022). LIRI-JP cohort data were obtained from the International Cancer Genome Consortium (ICGC) portal (, accessed on 27 March 2019). CHCC-HBV cohort data were acquired from The National Omics Data Encyclopedia (NODE) website (, accessed on 20 May 2021), and clinical data was obtained from the reference.¹⁹ Finally, after excluding patients with incomplete clinical information, 335 patients were obtained from TCGA-LIHC, 233 patients from LIRI-JP, and 159 patients from CHCC-HBV.

PRGs selection

PRGs were selected from various sources. First, a group of PRGs were extracted from relevant literatures.^16,20–23 Second, the Gene Ontology (GO) database () was also applied to obtain PRGs, using the keyword “pyroptosis.” Third, the Molecular Signature Database v7.4 (MSigDB) was utilized to download an additional PRGs (Reactome pyroptosis). All PRGs from the above three sources included for further study are summarized in Table S1.

Identification of PRGs with differential expression

TCGA-LIHC cohort was investigated for genes showing distinct expressions between tumor and non-tumor tissues. The “DESeq2” R package was used, with a threshold of |log2Foldchange (log2FC) | > 0.1 and false discovery rate (FDR) <0.05, and instances with insufficient gene expression data were excluded.

Detection of the PRGs-based molecular categorization of TCGA-LIHC patients

To classify the TCGA-LIHC patients into distinct molecular subtypes based on differentially expressed PRGs (DEPRGs), we employed an unsupervised consensus clustering approach. This approach utilized the k-means algorithm, implemented through the ConsensusClusterPlus package.²⁴ The process of consensus clustering was meticulously designed to ensure robustness and reproducibility of the resulting subtypes.

Specifically, the clustering was conducted over 1000 bootstrap iterations, with each iteration randomly sampling 80% of the patient data. This bootstrapping approach enhances the stability and reliability of the clustering outcomes. To determine the most appropriate number of clusters (i.e., molecular subtypes), we relied on a comprehensive evaluation of several metrics:

• The relative change in area under the CDF Curves: This metric evaluates the rate of change in cluster stability as the number of clusters increases. A smaller rate of change suggests a more stable cluster configuration.
• The proportion of ambiguous clustering (PAC) Algorithm: PAC quantifies the degree of uncertainty or ambiguity in assigning samples to clusters. A lower PAC value indicates clearer distinction between clusters.
• Consensus heatmap analysis²⁵: By visualizing the consensus matrix as a heatmap, we assessed the clustering's consistency. The heatmap provides a qualitative measure of the clustering quality, where distinct blocks of high consensus indicate well-defined clusters.

By integrating these methods, we determined the optimal number of molecular subtypes rooted in pyroptosis-related gene expression patterns among TCGA-LIHC patients. This comprehensive and systematic approach ensures that the molecular categorization is robust, reproducible, and reflective of underlying biological differences associated with pyroptosis in liver cancer.

Gene set variation analysis (GSVA)

To analyze significantly enriched molecular pathways of pyroptosis subtypes, the GSVA package was used.²⁶ Furthermore, differential analysis of enrichment scores for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was conducted with the limma package,²⁷ identifying KEGG pathways with |log2FC| > 0.2 and FDR <0.05 as the most differentially enriched pathways between pyroptosis subtypes.²⁸

The immunogenomic features of TCGA-LIHC patients

Based on the gene expression profiles of TCGA-LIHC samples, we analyzed 29 immune-associated gene sets which represented the comprehensive immune activity of tumors, encompassing the types, functionalities, and molecular pathways of tumor infiltrating immune cells.²⁹ Single-sample gene-set enrichment analysis (ssGSEA) was used for quantification of immune gene sets' enrichment.²⁶

Development and validation of the pyroptosis subtype-based prognostic model

To distinguish the differentially expressed genes (DEGs) between the pyroptosis subtypes, we examined the TCGA-LIHC cohort using the “DESeq2” R package, ensuring FDR <0.05 and |log2FC| > 1. We evaluated the predictive significance of DEGs in the cohort using random forest Cox regression analyses. Univariate Cox regression analyses with a p-value cutoff of 0.2 were then employed to assess the prognostic significance. We employed the LASSO-Cox regression model with the glmnet package to screen the candidate genes and construct a predictive model. This process identified the key pyroptosis subtype-related genes (PSRGs) for predicting prognosis in patients with HCC. The risk score, denoted as “pyroptosis subtype-related score (PSRS),” was then computed in the TCGA-LIHC cohort. The PSRS score was calculated using the following formula:$$$ \mathrm{PSRS}\ \mathrm{score}={\mathbf{\sum}}_{k-1}^n{X}_k\times {Y}_k\ \left(X:\mathrm{coefficients},Y:\mathrm{gene}\ \mathrm{expression}\ \mathrm{level}\right) $$$

Patients with HCC were classified into the high-PSRS and low-PSRS groups using the median PSRS, and the difference in overall survival (OS) between the two groups was analyzed. To validate the findings, the patients with HCC in LIRI-JP and CHCC-HBV datasets was also analyzed, where the formula was utilized to compute the PSRS score in these cohorts. Patients with HCC were then classified into high-PSRS and low-PSRS categories based on the median risk score, and their OS rates were compared. Receiver operating characteristic (ROC) curves were generated using the “timeROC,” “survminer” and “survival” R packages. To determine whether the PSRS serves as an independent prognostic factor, both univariate and multivariate Cox regression analyses were performed, taking into account covariates such as age, gender, grade, stage, and alpha-fetoprotein (AFP).

Construction of nomogram and calibration curves

The nomogram was constructed to predict individual survival probability using the “RMS” package. We also generated calibration curves to predict the 1-, 3-, and 5-year survival rates of patients with HCC.

Drug sensitivity analysis

Using the oncoPredict package of R software, the sensitivity score of each small molecular compound was computed for every patient in the high-PSRS and low-PSRS groups.

Statistical analyses

The study employed the independent Student's t-test and Chi-squared test for continuous and categorical data, respectively, to conduct pairwise comparisons across groups. The Mann–Whitney U-test was used to evaluate gene expression, immune cell infiltration, and immune pathway activation between nontumor and tumor tissues, while the Kruskal–Wallis test was used for comparisons across multiple groups. R 4.1.3 was used for all statistical analyses, and a two-tailed p-value of less than 0.05 was considered statistically significant unless stated otherwise. Odds ratios (ORs), hazard ratios (HRs), and 95% confidence intervals (CIs) were reported when necessary.

RESULTS

Identification of differentially expressed PRGs in TCGA-LIHC patients

The study's workflow is shown in Figure 1. First, 33 PRGs were extracted from previous researches.^16,20–23 Besides, 27 and 27 genes were identified respectively in the GOBP_PYROPTOSIS and REACTOME_PYROPTOSIS gene sets. After the removal of duplicated genes from the three gene sets, 61 genes were defined as PRGs (Figure 2A). The expression of these 61 PRGs were evaluated in non-tumor and tumor samples in the TCGA-LIHC dataset. By the criteria of FDR <0.05 and |log2FC| > 0.5, 20 DEPRGs were identified. Among them, 14 genes were upregulated, while 6 genes were downregulated in TCGA-LIHC samples, as demonstrated by volcano plots and heatmaps (Figure 2B,C).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Functional enrichment analysis

To better understand the functions of the DEPRGs, GO and KEGG analyses were conducted using the clusterProfiler package.^30–32 The GO enrichment analysis indicated that the DEPRGs were primarily associated with the regulatory processes that control interleukin-1 beta production, inflammasome complex, and phosphatidylserine binding, etc (Figure S1A). Furthermore, the KEGG analysis demonstrated that these DEPRGs participated in a variety of biological processes, including NOD-like receptor signaling pathway, necroptosis, apoptosis, etc (Figure S1B). These results suggest that DEPRGs are involved in numerous biological pathways beyond pyroptosis.

Identification of two pyroptosis subtypes characterized by distinct survival outcomes, functional annotations, and clinical features

The expression patterns of 20 DEPRGs in TCGA-LIHC patients were analyzed using unsupervised consensus clustering to explore a new molecular characterization. Based on the relative change in the area under the CDF curve, the PAC algorithm, and the consensus heatmap, the optimal number of clusters was determined to be two (k value = 2) (Figure 3A–C). All patients with HCC were classified into two subgroups: Pyroptosis Subtype (PS) I (228 patients, 68.1%) and PS II (107 patients, 31.9%) (Figure 3D). The results of K–M survival analysis showed that PS I patients had significantly enhanced OS compared to those in the PS II group (HR = 0.473, log-rank p = 3.006 × 10⁻⁴) (Figure 3E). The median OS time of patients in the PS I group was longer than that of the PS II group (6.73 vs. 3.37 years). Univariate and multivariate Cox regression analyses were conducted to assess the prognostic significance of PS and various clinicopathological variables.³³ The results of these analyses showed that the PS was significantly correlated with OS in the TCGA-LIHC cohort (Table 1), demonstrating that the PS is an independent prognostic factor.

[IMAGE OMITTED. SEE PDF]

TABLE 1 Univariate and multivariate cox proportional hazards analysis of clinicopathological variables and PS for overall survival in the TCGA-LIHC cohort.

	TCGA-LIHC(n = 260)
Univariate analysis	Multivariate analysis
HR (95% CI)	p-value	HR (95% CI)	p-value
Age (older vs. young)	1.66 (1.00–2.77)	0.052	1.67 (0.98–2.84)	0.059
Gender (female vs. male)	1.49 (0.89–2.47)	0.127	1.30 (0.77–2.22)	0.329
Grade (G3/G4 vs. G1/G2)	1.51 (0.91–2.5)	0.113	1.28 (0.76–2.17)	0.357
Stage (stage III/IV vs. stage I/II)	1.44 (0.82–2.55)	0.205	–	–
AFP (≥400 vs. <400)	0.89 (0.49–1.61)	0.691	–	–
PS (II vs. I)	2.13 (1.29–3.50)	0.003	2.11 (1.26–3.54)	0.005

Furthermore, GSVA was carried out to gain insight into the molecular pathways and underlying mechanisms associated with the pyroptosis subtypes of TCGA-LIHC patients. Thirty-six molecular pathways that were differentially enriched were identified, comprising of 22 pathways positively correlated with PS I while 14 pathways positively associated with PS II (Figure 3F). PS I tumors were primarily correlated with the biosynthesis of unsaturated fatty acids, peroxisome function, fatty acid metabolism, and the PPAR signaling pathway, while PS II tumors were mostly associated with the cell cycle, DNA replication, mismatch repair, and spliceosome function.

Subsequently, a comparison of the demographic and clinicopathological characteristics of TCGA-LIHC patients between the PS I and II groups was also conducted (Table 2). The distribution of OS status, tumor grade, and AFP level were all markedly different between the two subtypes. However, there was no significant difference in the age, gender, tumor stage, Child-Pugh classification grade, and fibrosis Ishak score distribution of patients with HCC between the two subtypes.

TABLE 2 Demographics and clinicopathological features of HCC patients in the TCGA-LIHC cohort.

Variables		p-value
PS I (n = 228)	PS II (n = 107)
Age (years)	60.3 ± 13.6	59.3 ± 12.1	0.244
Gender	-	-	0.373
Male	147 (64.5%)	75 (70.1%)
Female	81 (35.5%)	32 (29.9%)
OS status	-	-	1.162 × 10−3
Alive	166 (72.8%)	58 (54.2%)
Dead	62 (27.2%)	49 (45.8%)
Child-Pugh classification grade	-	-	0.674
A	148 (90.2%)	66 (93.0%)
B	16 (9.8%)	5 (7.0%)
C	1 (−)	0 (−)
NA	63 (−)	36 (−)
Tumor stage		-	0.745
Stage I	115 (55.3%)	51 (49.5%)
Stage II	50 (24.0%)	27 (26.2%)
Stage III	40 (19.2%)	24 (23.3%)
Stage IV	3 (1.4%)	1 (1.0%)
NA	20 (−)	4 (−)
Grade	-	-	3.384 × 10−3
G1	37 (16.4%)	10 (9.5%)
G2	116 (51.6%)	43 (41.0%)
G3	68 (30.2%)	44 (41.9%)
G4	4 (1.8%)	8 (7.6%)
NA	3 (−)	2 (−)
Fibrosis Ishak Score	-	-	0.101
0—No fibrosis	56 (37.3%)	16 (29.1%)
1,2—Portal fibrosis	20 (13.3%)	11 (20.0%)
3,4—Fibrous speta	17 (11.3%)	9 (16.4%)
5—Nodular formation and incomplete cirrhosis	4 (2.7%)	5 (9.1%)
6—Established cirrhosis	53 (35.3%)	14 (25.5%)
NA	78 (−)	52 (−)
AFP	-	-	9.113 × 10−5
<400	157 (84.0%)	53 (61.6%)
≥400	30 (16.0%)	33 (38.4%)
NA	41 (−)	21 (−)

The pyroptosis subtypes possessed distinct somatic mutation pattern

Previous studies have investigated the role of genomic alterations in influencing tumor immunity and immune infiltration profiles.^34,35 Hence, somatic mutational analysis was conducted to explore the distinct genomic variations between the two pyroptosis subtypes and revealed specific highly mutated genes in each pyroptosis subtype (Figure 4A,B). TP53 was the most frequently mutated gene in PS II (34%), whereas CTNNB1 was the most commonly mutated gene in PS I (25%). The somatic mutation analysis also identified three differentially mutated genes out of 13 genes after merging the top 10 mutated genes from each subtype. TP53 mutations were more frequent in the PS II cohort compared to the PS I cohort (34% vs. 24%; p = 8.662 × 10⁻²), which is similar to the findings of RYR1 and CSMD1 mutations. The frequencies of RYR1 (11% vs. 5%; p = 9.142 × 10⁻²) and CSMD1 (11% vs. 3%; p = 5.943 × 10⁻³) mutations in the PS II group were significantly higher than those in the PS I group (Figure 4C).

[IMAGE OMITTED. SEE PDF]

Construction of a prognostic model using PRGs in the training set

As shown in Figure 5A, we filtered 1869 DEGs in the two pyroptosis subtypes, using the “DESeq2” R package, with a threshold of |log2FC| > 1 and FDR <0.05. We preserved a set of 90 characteristic genes after screening them further through the random forest algorithm. We filtered out 47 of the 90 genes with a P-value of less than 0.2 through univariate Cox regression analysis. Subsequently, we constructed a pyroptosis-related prognostic model involving 20 genes through LASSO regression (Figure 5B,C).

[IMAGE OMITTED. SEE PDF]

This method allowed us to derive coefficients for each of the PRGs, thus constructing the PSRS. The formula for calculating the PSRS for each patient is as follows: PSRS = (−0.026 × IGFALS expression level (exp.)) + (0.096 × HAVCR1 exp.) + (0.151 × STAC exp.) + (0.167 × SMOX exp.) + (0.114 × GAL exp.) + (0.002 × GAGE1 exp.) + (0.158 × TFDP3 exp.) + (0.234 × RNF186 exp.) + (0.079 × P2RY6 exp.) + (0.122 × GCG exp.) + (0.555 × DAW1 exp.) + (−0.332 × SIRPG exp.) + (0.072 × PRAC2 exp.) + (0.069 × NELL1 exp.) + (0.038 × MSC exp.) + (0.001 × MACIR exp.) + (1.437 × LOXHD1 exp.) + (0.092 × EPO exp.) + (0.122 × PYDC1 exp.) + (−0.047 × PPARGC1A exp.). This formula provides a PSRS for each patient in the TCGA-LIHC dataset. We subsequently categorized the entire cohort into two groups based on the median value of the PSRS. This resulted in a high-PSRS group (n = 167) and a low-PSRS group (n = 168), enabling us to compare prognostic outcomes between these two cohorts more effectively (Figure 5D). The high-PSRS group showed higher mortality incidence and shorter OS time compared to the low-PSRS group (Figure 5E,F). The ROC analysis indicated that the accuracy of OS prognostics was 0.777 at 1 year (95% CI 0.731–0.824), 0.824 at 3 years (95% CI 0.790–0.858), and 0.747 at 5 years (95% CI 0.695–0.800) (Figure 5G). These findings suggest that our model shows high accuracy in predicting prognosis in the training set.

Validation of the prognostic model in the test sets

To assess the accuracy of the prognostic model we developed, we applied the identical formula to calculate PSRS in 233 patients with HCC from the LIRI-JP cohort. Based on the median PSRS, 117 patients were categorized into the low-PSRS group, and 116 patients into the high-PSRS group (Figure S2A). The high-PSRS group demonstrated higher mortality incidence and shorter OS time than the low-PSRS group (Figure S2B,C), which was consistent with the result of the training set. The ROC analysis revealed that the prognostic accuracy of OS was 0.727 at 1 year (95% CI 0.665–0.788), 0.663 at 2 years (95% CI 0.607–0.719), and 0.716 at 3 years (95% CI 0.662–0.771) (Figure S2D). Similarly, we computed PSRS for 159 patients with HCC from the CHCC-HBV cohort using the same formula as in the training set. Based on the median PSRS, 79 patients were categorized into the low-PSRS group, while the remaining 80 patients were placed in the high-PSRS group (Figure S3A). The patients in the high-PSRS group had higher mortality incidence and shorter OS time than those in the high-PSRS group (Figure S3B,C), which was consistent with the observations made in the training set. The ROC analysis demonstrated that the prognostic accuracy of OS was 0.716 at 1 year (95% CI 0.654–0.778), 0.680 at 2 years (95% CI 0.633–0.726), and 0.606 at 3 years (95% CI 0.554–0.658) (Figure S3D). These results further confirm the accuracy of our established prognostic model.

Independent prognosis analysis of PSRS and clinical characteristics

We conducted a comparative analysis of the PSRS among various clinicopathological characteristic groups in TCGA-LIHC patients. Figure S4A shows a positive correlation between the PSRS and grade (p = 1.256 × 10⁻⁴) and stage (p = 8.174 × 10⁻⁵), while no significant difference was observed in the other four groups (age, gender, Child-Pugh grade, and AFP). To evaluate the potential of PSRS and clinical features as independent prognostic factors, we performed univariate and multivariate analyses. The results of the univariate independent prognosis analysis showed a significant correlation between age, gender, grade, and PSRS with the OS (Figure S4B). Moreover, the multivariate independent prognosis analysis revealed that the PSRS could independently predict the prognosis of TCGA-LIHC patients (p < 0.001) (Figure S4C).

Construction of nomogram and calibration curves

To offer physicians with a better quantitative method to predict the prognosis of the patients with HCC, we established a nomogram model that combined demographic and clinical variables including age, gender, tumor grade, tumor stage, and PSRS based on TCGA-LIHC patients. According to Figure 6A, the PSRS emerged as a significant factor among various demographic and clinical variables in the prognostic model. Calibration curves demonstrated that the prognostic model had a significant correlation with the actual survival data of TCGA-LIHC patients (Figure 6B–D). Our nomogram model demonstrated significantly superior performance in comparison to traditional prognostic scoring systems, with consistently high areas under the ROC curve (AUC) values for 1, 3, and 5 years (Figure 6E–G). Thus, our study underscores that the nomogram model using PSRS provides an accurate prediction of the OS of patients with HCC.

[IMAGE OMITTED. SEE PDF]

Functional annotations in the high-PSRS and low-PSRS groups

GSVA enrichment analysis was conducted separately in the high- and low-PSRS groups in both the training and validation sets. Our findings revealed that in the training set, the low-PSRS group had enriched pathways such as beta-alanine metabolism, valine, leucine and isoleucine degradation, propanoate metabolism, fatty acid metabolism, limonene and pinene degradation, glycine, serine, and threonine metabolism, and primary bile acid biosynthesis (Figure S5A). Likewise, similar enriched pathways were also observed in the low-PSRS group in the validation set (LIRI-JP and CHCC-HBV). Furthermore, the high-PSRS group showed enriched pathways in DNA replication, cell cycle, and mismatch repair (Figure S5B,C).

Analysis of the tumor microenvironment in the high-PSRS and low-PSRS groups

We conducted ssGSEA on 29 immune signatures related to PSRS in TCGA-LIHC patients to investigate the underlying mechanisms. Among these, 14 signatures exhibited differential enrichment. Particularly, six anti-tumor signatures, including effector cells, effector cell traffic, NK cells, T cells, Th1 signature, and antitumor cytokines, were significantly enriched in the low-PSRS group. Conversely, eight pro-tumor signatures, including granulocyte traffic, myeloid cells traffic, protumor cytokines, cancer-associated fibroblasts, matrix, matrix remodeling, angiogenesis, and tumor proliferation rate, were highly enriched in the high-PSRS group (Figure S6). The findings suggest that the low-PSRS group mainly displayed a higher abundance in anti-tumor signatures, while the high-PSRS group exhibited increased enrichment in pro-tumor signatures. These results may offer insights into the differential prognosis associated with PSRS in TCGA-LIHC patients.

Evaluating the therapeutic response in the high-PSRS and low-PSRS groups

Using the oncoPredict algorithm on Genomics of Drug Sensitivity in Cancer (GDSC) database, we estimated chemotherapeutic response through half-maximal inhibitory concentration (IC50) data for TCGA-LIHC patients with high and low PSRS. Our findings identified 65 small molecular compounds with significantly different responses between the two groups (Table S2). Among these, the top four with lowest p values between high-PSRS and low-PSRS were Dasatinib_1079, IAP_5620_1428, SB505124_1194, and TAF1_5496_1732 (Figure S7A–D). Our findings suggest these small molecular compounds as potential agents for treating patients with HCC, with IAP_5620_1428, SB505124_1194, and TAF1_5496_1732 being more effective in low-PSRS tumors, and Dasatinib_1079 better suitable for high-PSRS tumors. Further analysis is necessary to determine the efficacy of these compounds as HCC treatment agents. In a word, our results provide potential molecular chemotherapy compounds for patients with HCC.

DISCUSSION

Our study aimed to explore pyroptosis subtypes based on PRGs and establish a prognostic model to predict the prognosis of patients with HCC. In recent years, the search for prognostic markers and prognostic models of HCC has garnered significant research attention.^36–40 These models present prognostic implications that would largely benefit patients with HCC for the prediction of the prognosis. Consistent with previous studies, the prognostic model we have formulated has shown commendable efficacy in predicting the prognosis of patients with HCC.

In our study, we analyzed the expression of 61 genes, reported in the literature and databases as PRGs, in HCC and non-tumor samples of the TCGA-LIHC dataset. The results indicated that 20 of these PRGs exhibited differential expression. Then, we applied consensus clustering, a widely used unsupervised clustering method,⁴¹ to stratify TCGA-LIHC patients based on the expression of 20 DEPRGs. Furthermore, we compared the clinicopathological parameters to establish the relationship between pyroptosis subtypes and clinical characteristics. Our results suggest that PS II subtype patients exhibited a significantly higher level of PRGs expression, poorer prognosis.

Our study of the signaling pathways related to the two subtypes of pyroptosis in HCC reveals potential mechanisms that contribute to their differing prognoses. Specifically, in the PS I subtype, the activation of pathways such as the PPAR signaling pathway highlights its significant involvement in the regulation of liver metabolism. Dysregulation of PPARα and PPARγ has been strongly associated with the development of HCC. PPARγ has been identified as a tumor-suppressor gene in hepatocarcinogenesis, exerting its inhibitory effects on tumor cell growth through mechanisms such as suppression of cell proliferation, induction of G2/M arrest, and promotion of apoptosis.⁴² Activation of the peroxisome pathway has been associated with a more favorable prognosis for patients with HCC, with studies indicating that individuals with low peroxisomal function, particularly those with HCV-related HCC, are at increased risk of adverse clinical outcomes.⁴³ On the other hand, the PS II subtype exhibits activation of pathways related to cell cycle regulation, DNA replication, mismatch repair, and spliceosome function. The dysregulation of cell cycle regulation is frequently associated with the uncontrolled growth of HCC cells. Studies have shown a correlation between the autophagic degradation machinery and the cell-cycle regulator cyclin D1 in the pathogenesis of HCC tumors.⁴⁴ Additionally, the upregulation of MCM7, a crucial component in DNA replication initiation, is predominantly detected in HCC tumors and is strongly associated with unfavorable prognostic outcomes in patients with HCC.⁴⁵ These results are consistent with our own findings, indicating that distinct activation and inhibition of pathways may contribute to the prognostic differences between the two pyroptosis subtypes, warranting further exploration of the underlying mechanisms of gene expression within these pathways.

To facilitate the distinction between the two pyroptosis subtypes in clinical practice, we developed a novel prognostic model for patients with HCC using the 20 most important PSRGs. Furthermore, the model was further validated in two external cohorts. Additionally, we performed comprehensive analyses on the model scores (PSRS) and elucidated drug sensitivity analyses between the high-PSRS and low-PSRS groups. These results suggest that our analyses on the molecular subtypes are thorough, and the validated prognostic model is reliable and practical for clinical use and future research on HCC. Notably, the high-PSRS group exhibited a considerably worse prognosis compared to the low-PSRS group.

Regarding the mechanism, the impact of pyroptosis on the prognosis and progression of patients with HCC cannot be separated from tumor immunity, as our study suggests. Several immune-related biological functions, such as viral protein interaction with cytokine and cytokine receptor, were significantly enriched. Hepatitis B and C viruses have been established as the two main risk factors for HCC,⁴⁶ and they can affect HCC development in various ways, including inducing inflammation.⁴⁷ Cytokines play a crucial role in promoting HCC carcinogenesis and progression, released in response to infection, inflammation, and carcinogen-induced injury.⁴⁸ Cytokine activity has also been identified as a key indicator of the severity and development of hepatitis B or C virus infections.⁴⁹ According to the ssGSEA results in this study, the antitumor cytokine scores of the low-PSRS group were found to be significantly higher than those of the high-PSRS groups. Conversely, the protumor cytokine scores were observed to be higher in the high-PSRS groups as compared to the low-PSRS group. These findings suggest a potential correlation between the PSRS and the immune response, wherein the high-PSRS group is associated with a predominant protumor cytokine milieu, while the low-PSRS group exhibits a more pronounced antitumor cytokine expression profile. Further investigation is required to confirm and validate these observations in a larger cohort of subjects.

During the drug sensitivity analysis, a total of 65 small molecular compounds were identified to have significant differences in sensitivity between the two PSRS groups. Among these compounds, Dasatinib_1079, IAP_5620_1428, SB505124_1194, and TAF1_5496_1732 exhibited the most pronounced differences. Notably, each of these compounds has been linked to potential therapeutic effects for HCC. Dasatinib is a tyrosine kinase inhibitor used to treat certain types of leukemia. A previous study found that the axitinib, erlotinib, and dasatinib mixture inhibited tumor growth in both SW620 and HT29 CAM tumors.⁵⁰ IAP-5620 targets components such as TRAF2 and RIPK1, which have been associated with poor outcomes in HCC, suggesting its potential utility in this context.⁵¹ SB505124, impacting the TGF-beta signaling pathway, presents a novel approach to targeting HCC, especially considering its association with ferroptosis and cancer metabolism.⁵² TAF1_5496_1732, by targeting TAF1, hints at the critical role of transcription regulation in HCC progression and offers a new angle for therapeutic intervention.⁵³ These results highlight the potential of utilizing PSRS score-based drug sensitivity analysis to discover novel therapeutic strategies for HCC. Although the identification of significant compounds with potential relevance to HCC treatment is promising, it is crucial to acknowledge the challenges associated with translating these findings into clinical applications. Additional research, encompassing experimental validation and clinical trials, is imperative to comprehensively evaluate the therapeutic efficacy and safety profiles of these compounds in the context of HCC.

In addition, our study has some limitations that should be considered in future researches. First, we only applied commonly used mRNA expression data for tumor classification in HCC subtyping. This leaves room for potential underestimation or overestimation of the subtypes, calling for more comprehensive data sources, such as gene mutations, to be considered in future investigations. Second, the lack of experimental validation of the model-related genes leads to ambiguity about their functionality in HCC, which should be addressed in future studies. Last, due to the unavailability of a similar pyroptosis-related prognostic model based on pyroptosis subtype, we could not undertake a comparative analysis to evaluate the superiority of our model. Therefore, future researches should aim to validate our model's effectiveness and explore its potential for further optimization.

CONCLUSION

Here, our investigation has demonstrated the close correlation between pyroptosis and HCC. Our study has established the PSRS as an independent risk factor for the prediction of HCC prognosis in the TCGA-LIHC, LIRI-JP, and CHCC-HBV datasets based on the 20 pyroptosis subtype-related genes. Our findings have also revealed the relationship between the PSRS and tumor immunity, as well as the identification of potential molecular compounds for HCC treatment. This study offers a promising direction for the identification of novel predictive markers for HCC prognosis and serves as a critical foundation for future investigations into the intricate relationship between pyroptosis and HCC immunity.

AUTHOR CONTRIBUTIONS

Haisheng Xie: Conceptualization (lead); data curation (lead); formal analysis (equal); methodology (equal); software (lead); supervision (equal); writing – original draft (lead); writing – review and editing (equal). Guanlin Huang: Conceptualization (equal); data curation (equal); formal analysis (equal); methodology (equal); software (equal); validation (lead); writing – review and editing (equal). Haoming Mai: Formal analysis (equal); software (equal); validation (equal); writing – review and editing (equal). Jiaxuan Chen: Funding acquisition (equal); methodology (equal); validation (equal); writing – review and editing (equal). Rong Na: Validation (equal); writing – review and editing (equal). De-Ke Jiang: Formal analysis (lead); funding acquisition (lead); methodology (lead); project administration (lead); resources (lead); supervision (lead); writing – review and editing (lead).

ACKNOWLEDGMENTS

We acknowledge TCGA database, CHCC-HBV dataset, and ICGC database for generously providing their platforms and to the diligent contributors for uploading their invaluable datasets.

FUNDING INFORMATION

This research was funded by the National Key Research and Development Program of China (No. 2022YFC2303600), the General Programs from the National Natural Science Foundation of China (No. 82272765, 81472618, 81670535, and 82203253), the Local Innovative and Research Teams Project of Guangdong Pearl River Talents Program (No. 2017BT01S131), the General Program from the Natural Science Foundation of Guangdong Province (No. 2022A1515220090, and 2019A1515011423), the Innovative Research Team Project of Guangxi Province (No. 2017GXNSFGA198002), the Dean Fund of Nanfang Hospital, Southern Medical University (No. 2018Z005), the Grant for Recruited Talents to Start Scientific Research from Nanfang Hospital, China, and the Outstanding Youths Development Scheme of Nanfang Hospital at Southern Medical University, China (No. 2017J001), and the Postdoctoral Science Foundation of China (No. 2022M711508).

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The study relies on publicly available data from the TCGA database, CHCC-HBV dataset and ICGC database, which were used and analyzed. Further data are available from authors upon reasonable request.

ETHICS STATEMENT

This article does not contain any studies with human participants or animals performed by any of the authors.

References

Villanueva A. Hepatocellular carcinoma. N Engl J Med. 2019;380(15):1450‐1462.

Sia D, Villanueva A, Friedman SL, Llovet JM. Liver cancer cell of origin, molecular class, and effects on patient prognosis. Gastroenterology. 2017;152(4):745‐761.

Nault JC, Villanueva A. Intratumor molecular and phenotypic diversity in hepatocellular carcinoma. Clin Cancer Res. 2015;21(8):1786‐1788.

Jemal A, Ward EM, Johnson CJ, et al. Annual report to the nation on the status of cancer, 1975‐2014, featuring survival. J Natl Cancer Inst. 2017;109(9):djx030.

Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. CA Cancer J Clin. 2023;73(1):17‐48.

Li X, Xu W, Kang W, et al. Genomic analysis of liver cancer unveils novel driver genes and distinct prognostic features. Theranostics. 2018;8(6):1740‐1751.

Schulze K, Imbeaud S, Letouzé E, et al. Exome sequencing of hepatocellular carcinomas identifies new mutational signatures and potential therapeutic targets. Nat Genet. 2015;47(5):505‐511.

Totoki Y, Tatsuno K, Covington KR, et al. Trans‐ancestry mutational landscape of hepatocellular carcinoma genomes. Nat Genet. 2014;46(12):1267‐1273.

Wallach D, Kang TB, Dillon CP, Green DR. Programmed necrosis in inflammation: toward identification of the effector molecules. Science. 2016;352(6281): [eLocator: aaf2154].

Shi J, Gao W, Shao F. Pyroptosis: gasdermin‐mediated programmed necrotic cell death. Trends Biochem Sci. 2017;42(4):245‐254.

Rébé C, Derangère V, Ghiringhelli F. Induction of pyroptosis in colon cancer cells by LXRβ. Mol Cell Oncol. 2015;2(1): [eLocator: e970094].

Dunn JH, Ellis LZ, Fujita M. Inflammasomes as molecular mediators of inflammation and cancer: potential role in melanoma. Cancer Lett. 2012;314(1):24‐33.

Hu B, Elinav E, Huber S, et al. Inflammation‐induced tumorigenesis in the colon is regulated by caspase‐1 and NLRC4. Proc Natl Acad Sci USA. 2010;107(50):21635‐21640.

Tang R, Xu J, Zhang B, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol. 2020;13(1):110.

Wang YY, Liu XL, Zhao R. Induction of pyroptosis and its implications in cancer management. Front Oncol. 2019;9:971.

Xia X, Wang X, Cheng Z, et al. The role of pyroptosis in cancer: pro‐cancer or pro‐“host”? Cell Death Dis. 2019;10(9):650.

Network. TCGAR. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017;169(7):1327‐1341.

Fujimoto A, Furuta M, Totoki Y, et al. Whole‐genome mutational landscape and characterization of noncoding and structural mutations in liver cancer. Nat Genet. 2016;48(5):500‐509.

Gao Q, Zhu H, Dong L, et al. Integrated proteogenomic characterization of HBV‐related hepatocellular carcinoma. Cell. 2019;179(2):561‐577.

Karki R, Kanneganti TD. Diverging inflammasome signals in tumorigenesis and potential targeting. Nat Rev Cancer. 2019;19(4):197‐214.

Wang B, Yin Q. AIM2 inflammasome activation and regulation: a structural perspective. J Struct Biol. 2017;200(3):279‐282.

Man SM, Kanneganti TD. Regulation of inflammasome activation. Immunol Rev. 2015;265(1):6‐21.

Zheng M, Karki R, Vogel P, Kanneganti TD. Caspase‐6 is a key regulator of innate immunity, Inflammasome activation, and host defense. Cell. 2020;181(3):674‐687.

Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572‐1573.

Șenbabaoğlu Y, Michailidis G, Li JZ. Critical limitations of consensus clustering in class discovery. Sci Rep. 2014;4:6207.

Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA‐seq data. BMC Bioinformatics. 2013;14:7.

Smyth GK, Michaud J, Scott HS. Use of within‐array replicate spots for assessing differential expression in microarray experiments. Bioinformatics. 2005;21(9):2067‐2075.

Yang C, Huang X, Liu Z, Qin W, Wang C. Metabolism‐associated molecular classification of hepatocellular carcinoma. Mol Oncol. 2020;14(4):896‐913.

He Y, Jiang Z, Chen C, Wang X. Classification of triple‐negative breast cancers based on immunogenomic profiling. J Exp Clin Cancer Res. 2018;37(1):327.

Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25‐29.

Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45(D1):D353‐D361.

Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovations. 2021;2(3): [eLocator: 100141].

Wang Z, Guo X, Gao L, Wang Y, Ma W, Xing B. Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients. Aging. 2020;12(18):18297‐18321.

Thorsson V, Gibbs DL, Brown SD, Wolf D. The immune landscape of cancer. Immunity. 2018;48(4):812‐830.

Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48‐61.

Fu XW, Song CQ. Identification and validation of Pyroptosis‐related gene signature to predict prognosis and reveal immune infiltration in hepatocellular carcinoma. Front Cell Dev Biol. 2021;9: [eLocator: 748039].

Wu Q, Jiang S, Cheng T, Xu M, Lu B. A novel Pyroptosis‐related prognostic model for hepatocellular carcinoma. Front Cell Dev Biol. 2021;9: [eLocator: 770301].

Zhang S, Li X, Zhang X, Zhang S, Tang C, Kuang W. The Pyroptosis‐related gene signature predicts the prognosis of hepatocellular carcinoma. Front Mol Biosci. 2021;8: [eLocator: 781427].

Deng M, Sun S, Zhao R, et al. The pyroptosis‐related gene signature predicts prognosis and indicates immune activity in hepatocellular carcinoma. Mol Med. 2022;28(1):16.

Wang J, Wang Y, Steffani M, et al. Novel risk classification based on pyroptosis‐related genes defines immune microenvironment and pharmaceutical landscape for hepatocellular carcinoma. Cancer. 2022;14(2):1‐20.

Brannon AR, Reddy A, Seiler M, et al. Molecular stratification of clear cell renal cell carcinoma by consensus clustering reveals distinct subtypes and survival patterns. Genes Cancer. 2010;1(2):152‐163.

Yu J, Shen B, Chu ES, et al. Inhibitory role of peroxisome proliferator‐activated receptor gamma in hepatocarcinogenesis in mice and in vitro. Hepatology. 2010;51(6):2008‐2019.

Lupberger J, Croonenborghs T, Roca Suarez AA, et al. Combined analysis of metabolomes, proteomes, and transcriptomes of hepatitis C virus‐infected cells and liver to identify pathways associated with disease development. Gastroenterology. 2019;157(2):537‐551.

Wu SY, Lan SH, Wu SR, et al. Hepatocellular carcinoma‐related cyclin D1 is selectively regulated by autophagy degradation system. Hepatology. 2018;68(1):141‐154.

Qu K, Wang Z, Fan H, et al. MCM7 promotes cancer progression through cyclin D1‐dependent signaling and serves as a prognostic marker for patients with hepatocellular carcinoma. Cell Death Dis. 2017;8(2): [eLocator: e2603].

Castello G, Scala S, Palmieri G, Curley SA, Izzo F. HCV‐related hepatocellular carcinoma: from chronic inflammation to cancer. Clin Immunol. 2010;134(3):237‐250.

Michielsen PP, Francque SM, van Dongen JL. Viral hepatitis and hepatocellular carcinoma. World J Surg Oncol. 2005;3:27.

Budhu A, Wang XW. The role of cytokines in hepatocellular carcinoma. J Leukoc Biol. 2006;80(6):1197‐1213.

Akcam FZ, Tigli A, Kaya O, Ciris M, Vural H. Cytokine levels and histopathology in chronic hepatitis B and chronic hepatitis C. J Interf Cytokine Res. 2012;32(12):570‐574.

Berndsen RH, Swier N, van Beijnum JR, Nowak‐Sliwinska P. Colorectal cancer growth retardation through induction of apoptosis, using an optimized synergistic cocktail of axitinib, erlotinib, and dasatinib. Cancer. 2019;11(12): [eLocator: 1878].

Schneider AT, Gautheron J, Feoktistova M, et al. RIPK1 suppresses a TRAF2‐dependent pathway to liver cancer. Cancer Cell. 2017;31(1):94‐109.

Tao Y, Wei L, You H. Ferroptosis‐related gene signature predicts the clinical outcome in pediatric acute myeloid leukemia patients and refines the 2017 ELN classification system. Front Mol Biosci. 2022;9: [eLocator: 954524].

Yan H, Wang Q, Shen Q, et al. Identification of potential transcription factors, long noncoding RNAs, and microRNAs associated with hepatocellular carcinoma. J Cancer Res Ther. 2018;14(Supplement):S622‐s627.