Identification of a pyroptosis-related gene prognostic signature in patients with hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC), which accounts for the highest proportion of liver cancer cases, is one of the most common malignancies in humans and the fourth leading cause of cancer-related death worldwide (1,2). A study has shown that the pathogenic factors for HCC include hepatitis B and C viruses, alcohol abuse, dietary toxins such as aflatoxin and aristolochic acid, metabolic liver disease (mainly nonalcoholic fatty liver), and exposure to carcinogenic substances such as hexachlorobenzene (3). Over the past few decades, due to the dissemination of health knowledge and advanced screening procedures, including alpha-fetoprotein (AFP) detection and noninvasive imaging technology, significant progress has been made in the treatment of HCC (3-7). However, the survival rate of advanced HCC remains poor (8). Tumor-lymph node–metastasis (TNM) staging is widely applied in clinical neoplasm management to inform predictive risk assessment and treatment decisions. Nevertheless, due to the high degree of molecular heterogeneity, the risk of recurrence and death may vary greatly even among patients with the same clinicopathological characteristics (9,10). Therefore, there is an urgent need to develop novel prognostic signatures for more accurate risk stratification among patients with HCC.

Pyroptosis, also referred to as inflammatory cell death, is a type of gasdermin-mediated programmed necrosis (11-13) and has been increasingly recognized for its dual roles in cancer progression and immune response. As an innate immune mechanism, pyroptosis can suppress tumor growth by inducing anti-tumor immunity while simultaneously promoting tumor progression in certain inflammatory microenvironments (14). In cell pyroptosis, gasdermin is activated to form cell membrane pores, which differs from the mechanism of cell apoptosis. Subsequently, the intact cell membrane is destroyed under the dual action of proinflammatory cytokines and alarm factors, eventually leading to cell death (11). The initial research on pyroptosis was mainly focused on its function in fighting infections. More recently, a large number of studies have shown that pyroptosis plays a critical role in the evolution of human malignancies. For instance: the interactions between pyroptosis and cancers are highly complex (15). On the one hand, the growth of tumors can be inhibited by mediating the pyroptosis of cancer cells (16). On the other hand, pyroptosis provides nutritional and environmental support for the occurrence and development of neoplasms (17,18). Nonetheless, evidence indicating the role of pyroptosis in promoting cancer is mounting (17,19).

Research on pyroptosis-related genes (PRGs) in ovarian cancer, lung adenocarcinoma, and gastric cancer has been published recently (20-22). The molecular characteristics, clinical significance, and immune interaction of PRGs in the above-mentioned cancers have been elucidated by studies employing bioinformatics methods. These studies have demonstrated that PRGs can be used as a highly valuable gene signature for predicting the survival outcomes of patients with cancer.

However, the expression levels, clinical prognostic characteristics, and immune infiltration status of the PRGs in HCC have not yet been clarified. Therefore, in this study, we used gene expression data from The Cancer Genome Atlas (TCGA) to identify molecular the subtypes of HCC based on genes related to pyroptosis. To develop a prognostic risk model, we screened 30 genes associated with clinical outcomes from 35 PRGs using the least absolute shrinkage and selection operator (LASSO) Cox regression. The risk model was able to assess the survival status of patients with HCC and was externally validated in cohorts from the International Cancer Genome Consortium (ICGC). This risk score model can be used as an independent prognostic evaluation marker in patients with HCC. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-954/rc).

Methods

Datasets and preprocessing

The RNA-sequencing (RNA-seq) gene expression data of HCC and adjacent normal tissues and corresponding clinical prognostic traits were obtained from the liver hepatocellular carcinoma (LIHC) data set in TCGA database (https://portal.gdc.cancer.gov/). The validation group, including the gene expression data of HCC specimens with clinical characteristics, was obtained from the ICGC data portal (https://icgc.org/). After removal of data with no clinical information and patients with an overall survival (OS) <30 days, 343 HCC and 50 adjacent normal samples from TCGA cohort and 228 HCC samples from the ICGC cohort were included. The gene expression RNA-seq data of 110 normal human liver tissue samples were downloaded from the Genotype-Tissue Expression (GTEx) data set (https://xenabrowser.net/datapages). The somatic mutation profile was obtained from TCGA data portal Subsequently, we used the “maftools” package in R (The R Foundation of Statistical Computing) to process the somatic mutation profiles sorted in the form of mutation annotation format (MAF). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of differentially expressed PRGs

We identified 35 PRGs from the published literature (17,20-23), including GSDMD, CASP5, GSDME, GSDMB, PRKACA, SCAF11, NOD2, NLRP2, CASP3, CASP9, IL18, CASP4, ELANE, GPX4, GSDMC, PLCG1, PYCARD, AIM2, CASP8, TIRAP, NLRP6, NLRP7, NLRP3, IL6, CASP6, TNF, GSDMA, NOD1, NLRP1, PJVK, IL1B, NLRC4, CASP1, GZMA, and GZMB. In the differential gene expression analysis, there are differences between HCC samples and the normal tumor-adjacent tissue samples in TCGA database, which could have led to statistical errors in the results. Therefore, we also included normal liver tissue samples from the GTEx database. The already-processed data were collected from the UCSC Xena browser (https://xenabrowser.net/datapages/). Then, these data were corrected to eliminate the batch effect for subsequent analysis. Afterwards, the messenger RNA (mRNA) expression of PRGs was compared between HCC and normal samples via the Wilcoxon test.

Protein-protein interaction (PPI) networks and hub genes

The online database, Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (24) (v. 11.0; http://www.string-db.org), was used to visualize the PPI networks between the 35 PRGs. Cytoscape (25,26) (https://cytoscape.org/) was used to visualize the PPI network of PRGs and calculate hub genes.

Identification of pyroptosis-related subtypes

The RNA-seq expression profiles of 35 PRGs were used to perform k-means unsupervised clustering through the “ConsensusClusterPlus” R package, which was repeated 1,000 times (27,28).

Estimation of immune infiltration and immune-related pathway activity

In order to clarify the immune infiltration status of each specimen in TCGA dataset corresponding to the pyroptosis subtype, we used single-sample gene set enrichment analysis (ssGSEA) through the “GSVA” package based on immune cell gene sets and immune-related pathways obtained from the published literature (29,30) (supplementary tables available at https://cdn.amegroups.cn/static/public/jgo-2024-954-1.xlsx and https://cdn.amegroups.cn/static/public/jgo-2024-954-2.xlsx). In addition, we applied the ESTIMATE algorithm tool to assess each tumor specimen to calculate the immune score, stromal score, ESTIMATE score, and tumor purity (31).

Construction of the risk model based on PRGs

We used the LASSO penalty method to further screen the prognostic-related pyroptosis genes to construct a prognostic risk model. We then identified the best parameter λ through 10-fold cross-validation. According to the optimal cutoff value calculated by the “survminer” package, we divided TCGA cohort into the high-risk and low-risk groups based on the prognostic risk score obtained from the screened genes. In addition, we plotted time-dependent receiver operating characteristic (ROC) curves via the “timeROC” package to assess the sensitivity and specificity of the risk model in TCGA and ICGC cohorts. Kaplan-Meier survival plots were used to show the prognosis for predicting the OS in HCC.

Correlations of the risk score with immune infiltration

We applied the abovementioned ssGSEA method to estimate the difference in infiltration levels of 28 immune cells between the high-risk score group and the low-risk score group in TCGA cohort. However, using a single algorithm and marker gene set alone to evaluate tumor-infiltrating immune cells (TIICs) could have resulted in calculation errors. To avoid these problems, we additionally used the Tumor Immune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer/) to quantify the immune infiltration per specimen and generate more reliable results (32). The correlation analysis was performed between immune status and PRGs by using Pearson’s correlation coefficients.

Gene function and pathway enrichment analysis

WebGestalt (https://www.webgestalt.org/option.php) was used to characterize the function and pathway enrichment of the target genes (33). First, we used the overexpression enrichment analysis method on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases and input the genes of interest into the WebGestalt server. Subsequently, we applied the “ggplot2” package to plot the results.

Construction and validation of the prognostic nomogram

Univariate and multivariate Cox regression analyses were used to evaluate the survival prognosis associated with pyroptosis risk score and clinical characteristics. Subsequently, we integrated the risk score and clinical characteristics that demonstrated prognostic value to construct a prognostic nomogram. Calibration curves for 1-year, 3-year, and 5-year survival rates were applied to assess the deviation of the nomogram from the ideal model. The ROC curve was constructed, and the area under the curve (AUC) was calculated to evaluate the ability of the nomogram that contained the pyroptosis gene risk score and TNM stage to predict the survival outcomes of patients with HCC.

Statistical analysis

Data were analyzed with the R version 4.0.5 and R Bioconductor packages. Comparisons between two groups were evaluated with the t-test or Wilcoxon test. ROC curves were generated with the “timeROC” package, and the nomogram was developed with the “rms” package. OS was compared between the pyroptosis molecular subgroups and risk score groups with Kaplan-Meier curves and log-rank tests. Unless otherwise stated, all statistical tests were two-way, and a P value <0.05 was considered significant.

Results

Differential expression of PRGs between normal and tumor tissues

The gene expressions of normal tissues included normal human liver tissue samples in the GTEx and tumor-adjacent normal tissues TCGA, while tumor tissues were derived from TCGA. We identified 34 differentially expressed genes (DEGs) related to pyroptosis between 160 normal tissues and 343 tumor tissues (all P values <0.01) (Figure 1A). Among them, 28 genes (PLCG1, CASP1, CASP4, CASP5, CASP6, CASP8, CASP9, PJVK, ELANE, PRKACA, GPX4, GSDMB, TIRAP, GSDMD, GSDME, SCAF11, GZMB, IL18, IL1B, IL6, NOD2, NLRC4, NLRP1, TNF, NLRP2, NLRP3, NLRP6, and NOD1) were downregulated in the tumor group, while 6 genes (CASP3, GSDMA, GSDMC, GZMA, PYCARD, and AIM2) were upregulated in the tumor group.

Figure 1 Expressions and interactions of PRGs and identification of molecular subtypes of HCC. (A) The expression of 35 PRGs in tumor tissues and normal tissues. **, P<0.01; ***, P<0.001; ****, P<0.0001. (B) PPI network showing the interactions of the pyroptosis-related genes. Hub genes: red. (C) Consensus matrix heatmap to cluster patients into 2 to 5 clusters showing the clustering stability after application of 1,000× k-means clustering. (D) Heatmap of the expression of PRGs in two molecular subtypes. FPKM, fragments per kilobase of transcript per million mapped reads; HCC, hepatocellular carcinoma; PRG, pyroptosis-related gene; PPI, protein-protein interaction; TNF, tumor necrosis factor.

In order to further characterize the interaction of the 34 PRGs at the protein level, we constructed a PPI network. The minimum interaction score required for PPI analysis was set to 0.4 (the default recommended value). Through the maximum neighborhood component (MNC), density of maximum neighborhood component (DMNC), maximal clique centrality (MCC) algorithms in Cytoscape, we determined that NLRP1, PYCARD, AIM2, IL18, CASP1, NLRC4, CASP4, CASP8, NLRP3, and TNF were hub genes (Figure 1B). These hub genes were also DEGs between the normal specimens and HCC.

Molecular subtypes of HCCs based on PRGs

To investigate the relationship between the expression profiles of 35 PRGs collected from other studies and the HCC subtypes, we conducted a consensus cluster analysis on all 343 HCC samples in TCGA training cohort. By adjusting the clustering variable (k), we observed that when k=2, the intergroup correlation was lowest and the intragroup correlation was highest compared to those derived from other k values, suggesting that based on 35 PRGs, 343 patients with HCC could be divided into two robust clusters (cluster A: n=59; cluster B: n=284) (Figure 1C,1D).

Prognostic differences and immune landscape of the HCC molecular subtypes

The prognosis signature among molecular subtypes was further analyzed. Analysis of OS showed that those in cluster A had a worse prognosis than did those in cluster B (Figure 2A). We applied ssGSEA to characterize the immune cell infiltration of the two molecular subtypes of the pyroptosis genes and observed that compared with cluster B, cluster A not only had a high level of infiltration of innate immune cells (macrophages, mast cells, and immature dendrites cells, and myeloid-derived suppressor cells [MDSCs]) but also of adaptive immune cells, such as activated B cells and activated CD4 T cells, effector memory CD4 T cells, and effector memory CD8 T cells (Figure 2B). Interestingly, cluster A exhibited a greater abundance of 17 immune-related pathways than that in cluster B (Figure 2C). In addition, through the ESTIMATE algorithm, we found that compared to cluster B, cluster A had higher stromal scores, immune scores, and ESTIMATE scores but a lower tumor purity score (Figure 2D), which was consistent with the above findings that cluster A had higher relative infiltration levels of immune cells compared with cluster B (Figure 2B,2C). Additionally, NOD2 had positive correlation with stromal score, immune score and ESTIMATE score; NLRP6 had negative correlation with immune score and GPX4 had negative correlation with stromal score (Figure S1A-S1E). Additionally, the expression of NLRP6, CASP8, NOD2, and GSDME had significant difference between high-risk and low-risk groups, not GPX4 (Figure S1F-S1J).

Figure 2 The different OS and immune profiles between the two molecular subtypes. (A) Kaplan-Meier curves of OS between the two clusters. (B) Discrepancies in immune cell infiltration between cluster A and cluster B. (C) Difference in expression of the immune-related pathways between cluster A and cluster B. (D) Distributions of the immune score, stromal score, ESTIMATE score, and tumor purity between cluster A and cluster B. The distance of both ends of the boxes represents the interquartile range of values, and the thick line represents the median value. *, P<0.05; ***, P<0.001; ****, P<0.0001. HR, hazard ratio; OS, overall survival.

We noticed that patients with poor prognoS1ses among the pyroptosis molecular subtypes had a greater degree of immune infiltration, revealing the complex relationship between the pyroptosis pathway and the immune microenvironment in HCC.

Development of a prognostic risk model based on training dataset

The prognostic risk prediction model of PRGs was generated based on TCGA cohort. First, we used LASSO Cox regression analysis to reduce the amount of PRGs. The minimum lambda (0.06544434) was obtained via 10-fold cross-validation to construct our five-gene prognostic model (Figure 3A,3B). The formula of the pyroptosis risk model was as follows: Risk score = 0.004182175 × GPX4 + 0.091401425 × CASP8 + 0.018589062 × NOD2 + 0.091845249 × GSDME − 0.031080493 × NLRP6. According to this formula, the upregulation of GPX4, CASP8, NOD2, and GSDME was associated with poor OS, while high expression of NLRP6 was associated with good prognosis (Figure S1K-S1O). The “survminer” package was used to obtain the optimal cutoff value of the risk score (cutoff = –0.01471409), with which patients in the TCGA cohort were divided into a high-risk group and low-risk group. Figure 3C shows the marked difference in the Kaplan-Meier survival curves between the two subgroups (P<0.001). Furthermore, as seen in Figure 3D, we visualized the survival status, population distribution trends, and the expression heatmap of the five genes that composed the risk score of the high-risk and low-risk groups. The AUC in the 1-year, 3-year, and 5-year OS in ROC curves were 0.63, 0.65, and 0.60, respectively (Figure 3E), indicating that the risk score had moderate sensitivity and accuracy in predicting the prognosis of those with HCC.

Figure 3 Construction of a prognostic PRG model in HCC. (A) LASSO and Cox regression coefficient selection and variable screening. (B) Cross-validation in the LASSO Cox regression model to select the tuning parameter. Lambda minimum = 0.06544434. (C) Kaplan-Meier curves of OS of patients in the high-risk and low-risk groups in the TCGA cohort. (D) Distribution of risk score, survival status, and the expression of five prognostic PRGs in TCGA cohort. (E) ROC curves for predicting the 1-, 3-, and 5-year OS in TCGA cohort. HCC, hepatocellular carcinoma; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; OS, overall survival; PRG, pyroptosis-related gene; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

External validation of the prognostic risk model

The ICGC cohort was used as an external validation dataset to verify the robustness of the risk model established via TCGA training cohort. Using the same method as described above, we divided the ICGC cohort into high-risk and low-risk groups based on the optimal cutoff value (cutoff =–0.00717451). The Kaplan-Meier survival curve indicated that the prognosis of the high-risk group was markedly worse than that of the low-risk group (Figure 4A). Figure 4B presents the overall trend of the ICGC cohort, including survival and death status, population distribution, and the expression levels of five prognostic PRGs, which are generally consistent with the data of TCGA cohort. In addition, the ROC curve generated based on the risk score was applied to evaluate the 1-year, 3-year, and 4-year OS in the validation cohort, as shown in Figure 4C.

Figure 4 Validation of the model in ICGC and associations of risk-score with immune characteristics in TCGA. (A) Kaplan-Meier curves of the OS of patients in the high-risk and low-risk groups in the ICGC cohort. (B) Distribution of risk score, survival status, and the expression of five prognostic PRGs in the ICGC cohort. (C) ROC curves for predicting 1-, 3-, and 4-year OS in the ICGC cohort. (D) Immune cell infiltration discrepancy between the high-risk and low-risk groups according to ssGSEA. (E) Immune cell infiltration discrepancy between the high-risk and low-risk groups in the TIMER database. (F) Relationships between molecular subtypes of pyroptosis and risk score. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. HR, hazard ratio; ICGC, International Cancer Genome Consortium; OS, overall survival; PRG, pyroptosis-related gene; ROC, receiver operating characteristic; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TIMER, Tumor Immune Estimation Resource.

Comparison of the immune activity between the risk score subgroups

We further investigated the difference in the abundance of immune cell infiltration between the high-risk group and the low-risk group. Through the “GSVA” package, we observed that CD56 dim natural killer cell and eosinophil infiltration in the low-risk group was relatively higher than that in the high-risk group. In contrast, immune cells with higher infiltration in the high-risk group included activated CD4 T cells, central memory CD4 T cells, effector memory CD4 T cells, γδ T cells, immature dendritic cells, MDSCs, memory B cells, natural killer T cells, plasmacytoid dendritic cells, follicular helper T cells, and type 2 helper T cells (Figure 4D). Subsequently, through use of the TIMER database, we found that compared with the low-risk group, the high-risk group had a higher degree of infiltration of B cells, CD4 T cells, CD8 T cells, neutrophils, macrophages, and dendritic cells (Figure 4E).

We next examined the relationships between the molecular subtypes of pyroptosis and the risk score. We found that cluster A had a significantly higher proportion of patients with high-risk scores, who were previously confirmed to have a poor prognosis (Figure 4F).

Analysis of enrichment pathways and genome changes between risk groups

We employed the Wilcoxon rank-sum test via R software to evaluate the DEGs between the high-risk and low-risk groups. The threshold was a false-discovery rate (FDR) <0.05 and fold change [expressed in the log2 (average fragments per kilobase per million mapped fragments ratio) ≥1. We consequently found that in the high-risk group in TCGA cohort, there were 2,629 significantly upregulated genes and 76 significantly downregulated genes (Figure S2). The details of the genes with high and low expression are provided in the supplementary tables available at https://cdn.amegroups.cn/static/public/jgo-2024-954-3.xlsx and https://cdn.amegroups.cn/static/public/jgo-2024-954-4.xlsx.

GO analysis showed that the genes with high expression in the high-risk group were mainly enriched in various synaptic signaling pathways, embryonic organ morphogenesis, and G protein−coupled receptor signaling pathway (Figure 5A). Meanwhile, genes with low expression in the high-risk group were mainly enriched in various biological metabolic processes, such as the xenobiotic metabolic process, lipid hydroxylation, and cellular response to the xenobiotic stimulus (Figure 5B). KEGG analysis revealed that upregulated genes were mainly enriched in neuroactive ligand−receptor interaction, nicotine addiction, and rheumatoid arthritis (Figure 5C), while downregulated genes were enriched primarily in pathways such as chemical carcinogenesis, retinol metabolism, and drug metabolism (Figure 5D).

Figure 5 Comprehensive analyses of enriched pathways and genomic alterations between the different risk groups. (A,B) Gene Ontology enrichment analysis was performed with significantly upregulated and downregulated genes. (C,D) Kyoto Encyclopedia of Genes and Genomes analysis was performed with significantly upregulated and downregulated genes. (E,F) Top 20 most frequently mutated genes in the high-risk cohort and low-risk cohort. (G) TP53 was the most differentially mutated gene between the high-risk cohort and low-risk cohort. (H) Lollipop plot showing the different mutation sites of TP53 between the two cohorts. (I) Kaplan-Meier analysis showing the independent association between OS and the TP53 mutation. (J,K) The heatmap illustrating the co-occurrence and mutually exclusive mutations of the top 25 frequently mutated genes in each cohort. **, P<0.01; ***, P<0.001; NS, not significant. FDR, false discovery rate; HR, hazard ratio; OS, overall survival; WT, wild type.

As shown in Figure 5E,5F, the 20 most frequently mutated genes in the high-risk and low-risk cohorts were not wholly consistent, which may suggest differences in the genomic alteration of the risk score groups. Based on the frequency of mutation, the mutated genes in the high-risk cohort included TP53, TTN, CTNNB1, MUC16, PCLO, AXIN1, XIRP2, APOB, RYR1, ABCA13, CSMD3, FAT3, LRP1B, ARID1A, CACNA1E, CCDC168, RYR2, ABCB5, COL6A6, and FRAS1 (Figure 5E). The mutated genes in the low-risk cohort were CTNNB1, TTN, TP53, MUC16, ALB, MUC4, RYR2, ABCA13, APOB, OBSCN, PCLO, BAP1, USH2A, ARID2, FLG, HMCN1, LRP1B, SPTA1, CUBN, and PRKDC (Figure 5F). We used the Fisher exact test (with a threshold of P<0.05) to examine the differences between gene mutations and wild-type in the high- or low-risk groups (Figure 5G) and found TP53 mutations to be significantly correlated with pyroptosis in patients with HCC. The lollipop chart in Figure 5H demonstrates the different mutation sites of TP53 between the two cohorts, while Figure 5I shows a difference in the survival prognosis of patients with HCC associated with the TP53 mutation. In addition, as shown in Figure 5J,5K, in the co-occurrence and mutually exclusive mutations of the high and low-risk groups, we observed a particular TP53-CSDM1 co-occurrence mutation, suggesting that their respective mutations in HCC might exert a combined effect.

Establishment and assessment of a nomogram based on the risk score

Univariate and multivariate Cox analyses confirmed the pyroptosis risk score as an independent prognostic marker for predicting the survival status of patients with HCC (Table 1). We further integrated the risk score and prognostic clinical characteristics to develop a nomogram for predicting the 1-year, 3-year, and 5-year OS rates of patients with HCC (Figure 6A). The AUC of the time-dependent ROC curve indicated that the performance of the nomogram prediction was better than that of the risk score alone (Figure 6B). The Kaplan-Meier curve for the OS of the nomogram was significant in nomo-low group compared with nomo-high group (P<0.001). The nomo-high group showed poor OS (Figure 6C). Furthermore, the calibration chart showed that compared with the ideal model, that the nomogram with the risk score and TNM stage performed well (Figure 6D). Additionally, to further investigate the core genes in risk model, we performed the competitive endogenous RNA (ceRNA) analysis based on these novel five genes in risk model. Our results showed that there were 7 microRNAs (miRNAs) obtained according to three database including starbase, miRDB, miRanda (Figure S3A). Then, we screened potential long non-coding RNAs (lncRNAs) associated with 7 miRNAs using ENCORI database and the ceRNA network was visualized using cytoscape (3.8.2) software (Figure S3B).

Figure 6 Development of the nomogram in TCGA cohort. (A) Nomogram for predicting the 1-year, 3-year, and 5-year OS of patients with hepatocellular carcinoma. (B) ROC curves and AUC for the 1-year, 3-year, and 5-year survival of the nomogram. (C) Kaplan-Meier curves of the OS of the nomogram. (D) Calibration plot of the nomogram for predicting the 1-year, 3-year, and 5-year OS. AUC, area under the curve; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Discussion

HCC is a highly malignant liver cancer with rapid progression, often diagnosed at advanced stages, leading to poor prognosis and limited therapeutic options (34,35). The progression from early HCC is insidious, and thus most patients with HCC are at the advanced stage when they are diagnosed, which means that most have already lost the opportunity to undergo surgical liver resection. Although radiotherapy, chemotherapy, liver transplantation, and other potential treatment methods can yield a variety of therapeutic effects and prolong the survival time of patients, the prognosis of those with HCC remains unsatisfactory owing to the high risk of recurrence and intrahepatic spread (36). Pyroptosis is a unique type of programmed cell death characterized by the swelling and lysis of cells and the release of many proinflammatory factors (37). The relationship between malignancies and pyroptosis is complex, and the effects of pyroptosis vary according to the different human tissues, genetic landscapes, and cancers (38). With the continued progress in the research on tumor molecular biology, prognostic gene signatures reflecting malignancy progression have been increasingly applied in predicting survival in HCC, which may help achieve more accurate individualized treatment management (39,40). Our study contributes to this field by providing insights into the relationship between pyroptosis genes and HCC.

We determined the mRNA levels of 35 PRGs derived from the HCC literature and of normal tissues and discovered that the majority of them are differentially expressed. Compared with that in healthy liver tissue, the expression of 28 PRGs and 6 PRGs was lower in HCC tissues and higher in HCC, respectively. Subsequently, we identified PYCARD, AIM2, IL18, CASP1, CASP4, CASP8, NLRC4, NLRP1, NLRP3, and TNF as the hub genes in the pyroptosis pathway.

We additionally applied the unsupervised clustering algorithm to distinguish pyroptosis-related subtypes (cluster A and cluster B). Cluster A was overexpressed in immune cells and immune-related pathways, and its associated prognosis was worse than that of cluster B. Given the prognostic capacity demonstrated by these abovementioned PRGs, we employed LASSO Cox regression analysis to develop a predictive risk model based on five PRGs (GPX4, CASP8, NOD2, GSDME, and NLRP6). This model had high accuracy and practicability in predicting the 1-year, 3-year, and 5-year OS of patients with HCC. Subsequently, validation of the gene signature’s prognostic potential in the ICGC cohort was done using the same formula and statistical methods. The results indicated that for patients with HCC, GPX4, CASP8, NOD2, and GSDME were risk factors while NLRP6 was a protective factor.

In recent years, glutathione peroxidase 4 (GPX4) has been considered to be a vital regulatory factor of ferroptosis tumor cell death. Unlike other GPXs, GPX4 has the function of catalyzing the reduction of lipid peroxides in a complex cell membrane environment, indicating that GPX4 has a distinct role in physiology (41). Inhibition of GPX4 can induce the death of cancer cells that are resistant to conventional chemotherapy or radiotherapy (42-44). In our study, we discovered that the upregulation of GPX4 is associated with poor prognosis and a high risk of HCC in patients. Thus, inhibiting the high expression of GPX4 may represent a novel therapeutic strategy for reducing the mortality rate in patients with HCC.

Caspase-8 (CASP8) is linked to an essential enzyme in the pyroptosis pathway (45,46). The significance of CASP8 in initiating necrosis receptor–induced pyroptosis and maintaining immune homeostasis and surveillance has been well established. It has long been held that the expression levels of CASP8 in most tumors are usually downregulated, leading to pyroptosis escape and enhanced resistance to anticancer therapies (47). Multiple studies (48-50) have shown that compared with those in healthy liver tissues, the expression levels of CASP8 are downregulated in HCC tissues, which is consistent with our results. The gasdermin protein family, particularly GSDMD and GSDME, plays a pivotal role in pyroptotic cell death through pore formation in the cell membrane. These proteins are activated by inflammatory caspases, such as caspase-1, -4, and -11, and their cleavage results in the release of IL-1β and IL-18, which could amplify local immune responses (51,52). However, the relationship between the levels of CASP8 expression and the prognosis of patients with HCC has not yet been clarified. Our analysis revealed that the upregulation of CASP8 could increase the mortality risk in patients HCC; this may be related to the RIPK1/FADD/CASP8 signal complex causing liver parenchymal cell damage, thereby impairing the normal functioning of the liver (53).

Nucleotide-binding oligomerization domain 2 (NOD2), a recognized innate immune sensor, has a significant effect on carcinogenesis (54). It has been reported that NOD2 dysregulation is involved in the pathogenesis of Crohn’s disease (CD) and colitis-related colon cancer. The NOD2 gene polymorphism is associated with colorectal cancer, lymphoma, lung cancer, gastric cancer, ovarian cancer, laryngeal cancer, and breast cancer (55). In their research on HCC, Ma et al. showed that NOD2 restrains tumorigenesis and enhances tumor chemotherapeutic sensitivity through targeted regulation of the AMPK pathway, which may reveal a novel strategy for the treatment of HCC based on NOD2 adjustment (56). A recent study on HCC demonstrated that NOD2 initiates the nuclear autophagy pathway, leading to DNA damage and increased genomic instability; moreover, NOD2 was found to be a poor prognostic factor closely related to the increased risk of death in patients with HCC (57). Our finding that the overexpression of NOD2 is associated with a poor prognosis and high-risk score is consistent with the above research.

GSDME, also known as DFNA5, was first identified as a mutated gene inherited in an autosomal-dominant manner, causing the loss of cochlear hair cells, which eventually leads to progressive, sensorineural, and nonsyndromic hearing impairment. In recent years, multiple studies have shown that GSDME, a member of the gasdermin family, is closely related to cancer (11,17,58). Due to epigenetic suppression of methylation, the expression of GSDME is downregulated in gastric cancer, colorectal cancer, breast cancer, and most human cancer cell lines (11) compared with that in healthy human tissues. This is consistent with the findings of our study. GSDME is a molecular drug target for treating human malignancies, and clinicians can select appropriate chemotherapeutic drugs according to its expression levels to improve the sensitivity of chemotherapeutic drugs and reduce drug resistance (17).

NLRP6 (originally known as PYPAF5) is a novel node-like receptor (NLR) family member that produces inflammasomes. NLRP6 consists of an N-terminal pyrin domain, a central NACHT domain, and a terminal leucine-rich repeat (LRR) domain (59). The majority of research on NLRP6 has concentrated on the intestine in mouse models, and only a few studies have involved human patients (59,60). In the study by Domblides et al., the downregulation of NLRP6 was related to the poor prognosis of patients with colorectal cancer (61). Our research indicated that NLRP6 is a protective factor for the prognosis of patients with HCC, and this offers a direction of further research into NLRP6 in HCC.

The tumor microenvironment (TME) plays a vital role in antitumor molecular therapy. Therefore, we also sought to characterize the abundance of immune and stromal cells in the high-risk and low-risk groups. Interestingly, our findings revealed a higher degree of immune and stromal cell infiltration in the high-risk group compared with the low-risk one, unlike the observations in ovarian cancer and gastric cancer (21,22). Our findings support the well-known evidences in literature, suggesting that, through this mechanism, pyroptosis contributes to the modulation of the TME and the persistent release of damage-associated molecular patterns (DAMPs) could paradoxically result in immune suppression, promoting HCC progression (62). This may suggest a more complex mechanism underlying the influence of pyroptosis genes on the microenvironment of hepatocellular tumors. The low-risk group was associated with various biological metabolic pathways, while the high-risk group was mainly associated with signal transduction pathways. It is worth noting that the high-risk group was also significantly associated with nicotine addiction and autoimmune diseases (rheumatoid arthritis). By comparing the genomic changes in different risk groups, we observed that the high-risk group was significantly associated with more aggressive molecular alterations (such as TP53 mutations). More importantly, the univariate and multivariate Cox analyses confirmed that the nomogram incorporating the pyroptosis risk score and TNM stage had good prediction ability and may serve as a meaningful indicator in personalized medicine. Lastly, the pyroptosis-related mechanisms may, in the future, offer valuable insights into therapeutic applications. Pyroptosis-inducing agents, including certain chemotherapy and immune checkpoint inhibitors, have demonstrated the potential to enhance tumor sensitivity to treatment by overcoming resistance to apoptosis, providing a promising avenue for HCC therapy (63). Additionally, a preclinical study using animal models could validate the therapeutic efficacy of targeting pyroptosis-related pathways in combination with immune checkpoint blockade (64).

Our study has several limitations which should be mentioned. First, the prognostic model for the pyroptosis risk score was constructed and validated using retrospective cohorts, which may limit its applicability in prospective clinical settings. Second, the heterogeneity of the TCGA and ICGC datasets, including variations in patient demographics and treatment histories, may introduce biases. Furthermore, we are validating the expression of these genes in liver cancer tissues and the biological mechanisms linking PRGs and immune infiltration require validation through in vivo and in vitro experiments in the future.

Future in vivo and in vitro studies should also assess the proposed model performance in independent cohorts and afterwards explore its potential role for guiding therapeutic strategies. We will prospectively collect clinical samples to further validate the prognostic ability of risk. Additionally, we will analyze the predictive ability of liver cancer efficacy before and after treatment, such as efficacy evaluation neoadjuvant therapy or conversion therapy. Finally, the molecular functions of key genes in risk models will be explored to uncover the occurrence and development of liver cancer through molecular experiments.

Conclusions

We delineated two pyroptosis molecules subtypes in patients with HCC patients that were associated with distinct prognoses. A risk model comprising five genes (GPX4, CASP8, NOD2, GSDME, and NLRP6) was constructed and demonstrated good prognostic performance. This risk prognosis model could be applied as a prognostic predictor for determining the survival status of patients with HCC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-954/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-954/prf

Funding: This study was supported by a grant from the Natural Science Foundation of Hunan Province, China (No. 2019JJ40422).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-954/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers 2021;7:6. [Crossref] [PubMed]
2. Yang X, Yang C, Zhang S, et al. Precision treatment in advanced hepatocellular carcinoma. Cancer Cell 2024;42:180-97. [Crossref] [PubMed]
3. Wang W, Lou W, Ding B, et al. A novel mRNA-miRNA-lncRNA competing endogenous RNA triple sub-network associated with prognosis of pancreatic cancer. Aging (Albany NY) 2019;11:2610-27. [Crossref] [PubMed]
4. Llovet JM, Montal R, Sia D, et al. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol 2018;15:599-616. [Crossref] [PubMed]
5. European Association for the Study of the Liver. Corrigendum to "EASL Clinical Practice Guidelines: Management of hepatocellular carcinoma" J Hepatol 2019;70:817. [J Hepatol 69 (2018) 182-236]. [Crossref] [PubMed]
6. Marrero JA, Kulik LM, Sirlin CB, et al. Diagnosis, Staging, and Management of Hepatocellular Carcinoma: 2018 Practice Guidance by the American Association for the Study of Liver Diseases. Hepatology 2018;68:723-50. [Crossref] [PubMed]
7. Llovet JM, De Baere T, Kulik L, et al. Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2021;18:293-313. [Crossref] [PubMed]
8. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
9. Craig AJ, von Felden J, Garcia-Lezana T, et al. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2020;17:139-52. [Crossref] [PubMed]
10. Zheng H, Pomyen Y, Hernandez MO, et al. Single-cell analysis reveals cancer stem cell heterogeneity in hepatocellular carcinoma. Hepatology 2018;68:127-40. [Crossref] [PubMed]
11. Liu X, Xia S, Zhang Z, et al. Channelling inflammation: gasdermins in physiology and disease. Nat Rev Drug Discov 2021;20:384-405. [Crossref] [PubMed]
12. Shi J, Gao W, Shao F. Pyroptosis: Gasdermin-Mediated Programmed Necrotic Cell Death. Trends Biochem Sci 2017;42:245-54. [Crossref] [PubMed]
13. Cui L, Xu X, Fan H, et al. Reuterin promotes pyroptosis in hepatocellular cancer cells through mtDNA-mediated STING activation and caspase 8 expression. Cancer Lett 2024;601:217183. [Crossref] [PubMed]
14. Tan Y, Chen Q, Li X, et al. Pyroptosis: a new paradigm of cell death for fighting against cancer. J Exp Clin Cancer Res 2021;40:153. [Crossref] [PubMed]
15. Xue Y, Enosi Tuipulotu D, Tan WH, et al. Emerging Activators and Regulators of Inflammasomes and Pyroptosis. Trends Immunol 2019;40:1035-52. [Crossref] [PubMed]
16. Zhou Y, Zhang W, Wang B, et al. Mitochondria-targeted photodynamic therapy triggers GSDME-mediated pyroptosis and sensitizes anti-PD-1 therapy in colorectal cancer. J Immunother Cancer 2024;12:e008054. [Crossref] [PubMed]
17. Xia X, Wang X, Cheng Z, et al. The role of pyroptosis in cancer: pro-cancer or pro-"host"? Cell Death Dis 2019;10:650. [Crossref] [PubMed]
18. Zhang H, Li L, Zhang Z, et al. Pyroptotic macrophages promote proliferation and chemotherapy resistance of peripheral T-cell lymphoma via TLR4 signaling pathway. Cancer Sci 2024;115:2444-60. [Crossref] [PubMed]
19. Broz P, Pelegrín P, Shao F. The gasdermins, a protein family executing cell death and inflammation. Nat Rev Immunol 2020;20:143-57. [Crossref] [PubMed]
20. Lin W, Chen Y, Wu B, et al. Identification of the pyroptosis-related prognostic gene signature and the associated regulation axis in lung adenocarcinoma. Cell Death Discov 2021;7:161. [Crossref] [PubMed]
21. Shao W, Yang Z, Fu Y, et al. The Pyroptosis-Related Signature Predicts Prognosis and Indicates Immune Microenvironment Infiltration in Gastric Cancer. Front Cell Dev Biol 2021;9:676485. [Crossref] [PubMed]
22. Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov 2021;7:71. [Crossref] [PubMed]
23. Karki R, Kanneganti TD. Diverging inflammasome signals in tumorigenesis and potential targeting. Nat Rev Cancer 2019;19:197-214. [Crossref] [PubMed]
24. Szklarczyk D, Franceschini A, Wyder S, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 2015;43:D447-52. [Crossref] [PubMed]
25. Doncheva NT, Morris JH, Gorodkin J, et al. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res 2019;18:623-32. [Crossref] [PubMed]
26. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
27. Bello-Orgaz G, Menéndez HD, Camacho D. Adaptive k-means algorithm for overlapped graph clustering. Int J Neural Syst 2012;22:1250018. [Crossref] [PubMed]
28. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
29. Jia Q, Wu W, Wang Y, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun 2018;9:5361. [Crossref] [PubMed]
30. Li Y, Jiang T, Zhou W, et al. Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nat Commun 2020;11:1000. [Crossref] [PubMed]
31. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
32. Li T, Fu J, Zeng Z, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res 2020;48:W509-14. [Crossref] [PubMed]
33. Liao Y, Wang J, Jaehnig EJ, et al. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res 2019;47:W199-205. [Crossref] [PubMed]
34. Rich NE. Changing Epidemiology of Hepatocellular Carcinoma Within the United States and Worldwide. Surg Oncol Clin N Am 2024;33:1-12. [Crossref] [PubMed]
35. Lehrich BM, Zhang J, Monga SP, et al. Battle of the biopsies: Role of tissue and liquid biopsy in hepatocellular carcinoma. J Hepatol 2024;80:515-30. [Crossref] [PubMed]
36. Wang X, Lu J. Immunotherapy for hepatocellular carcinoma. Chin Med J (Engl) 2024;137:1765-76. [Crossref] [PubMed]
37. Panganiban RA, Nadeau KC, Lu Q. Pyroptosis, gasdermins and allergic diseases. Allergy 2024;79:2380-95. [Crossref] [PubMed]
38. Liu P, Zhang Z, Chen H, et al. Pyroptosis: Mechanisms and links with diabetic cardiomyopathy. Ageing Res Rev 2024;94:102182. [Crossref] [PubMed]
39. Liu S, Meng Y, Zhang Y, et al. Integrative analysis of senescence-related genes identifies robust prognostic clusters with distinct features in hepatocellular carcinoma. J Adv Res 2024;S2090-1232(24)00150-4.
40. Xu K, Dai C, Yang J, et al. Disulfidptosis-related lncRNA signatures assess immune microenvironment and drug sensitivity in hepatocellular carcinoma. Comput Biol Med 2024;169:107930. [Crossref] [PubMed]
41. Brigelius-Flohé R, Maiorino M. Glutathione peroxidases. Biochim Biophys Acta 2013;1830:3289-303. [Crossref] [PubMed]
42. Kim DH, Kim WD, Kim SK, et al. TGF-β1-mediated repression of SLC7A11 drives vulnerability to GPX4 inhibition in hepatocellular carcinoma cells. Cell Death Dis 2020;11:406. [Crossref] [PubMed]
43. Shin D, Kim EH, Lee J, et al. Nrf2 inhibition reverses resistance to GPX4 inhibitor-induced ferroptosis in head and neck cancer. Free Radic Biol Med 2018;129:454-62. [Crossref] [PubMed]
44. Li J, Li Y, Wang D, et al. PLAG1 interacts with GPX4 to conquer vulnerability to sorafenib induced ferroptosis through a PVT1/miR-195-5p axis-dependent manner in hepatocellular carcinoma. J Exp Clin Cancer Res 2024;43:143. [Crossref] [PubMed]
45. Ni C, Ye Y, Wang M, et al. A six-nucleotide insertion-deletion polymorphism in the CASP8 promoter is associated with risk of coal workers' pneumoconiosis. J Toxicol Environ Health A 2009;72:712-6. [Crossref] [PubMed]
46. Ai Y, Meng Y, Yan B, et al. The biochemical pathways of apoptotic, necroptotic, pyroptotic, and ferroptotic cell death. Mol Cell 2024;84:170-9. [Crossref] [PubMed]
47. Mandal R, Barrón JC, Kostova I, et al. Caspase-8: The double-edged sword. Biochim Biophys Acta Rev Cancer 2020;1873:188357. [Crossref] [PubMed]
48. Cho S, Lee JH, Cho SB, et al. Epigenetic methylation and expression of caspase 8 and survivin in hepatocellular carcinoma. Pathol Int 2010;60:203-11. [Crossref] [PubMed]
49. Yamaguchi Y, Shiraki K, Fuke H, et al. Adenovirus-mediated transfection of caspase-8 sensitizes hepatocellular carcinoma to TRAIL- and chemotherapeutic agent-induced cell death. Biochim Biophys Acta 2006;1763:844-53. [Crossref] [PubMed]
50. Boege Y, Malehmir M, Healy ME, et al. A Dual Role of Caspase-8 in Triggering and Sensing Proliferation-Associated DNA Damage, a Key Determinant of Liver Cancer Development. Cancer Cell 2017;32:342-359.e10. [Crossref] [PubMed]
51. Yang F, Bettadapura SN, Smeltzer MS, et al. Pyroptosis and pyroptosis-inducing cancer drugs. Acta Pharmacol Sin 2022;43:2462-73. [Crossref] [PubMed]
52. Rao Z, Zhu Y, Yang P, et al. Pyroptosis in inflammatory diseases and cancer. Theranostics 2022;12:4310-29. [Crossref] [PubMed]
53. Kondylis V, Polykratis A, Ehlken H, et al. NEMO Prevents Steatohepatitis and Hepatocellular Carcinoma by Inhibiting RIPK1 Kinase Activity-Mediated Hepatocyte Apoptosis. Cancer Cell 2015;28:582-98. [Crossref] [PubMed]
54. Caruso R, Warner N, Inohara N, et al. NOD1 and NOD2: signaling, host defense, and inflammatory disease. Immunity 2014;41:898-908. [Crossref] [PubMed]
55. Kutikhin AG. Role of NOD1/CARD4 and NOD2/CARD15 gene polymorphisms in cancer etiology. Hum Immunol 2011;72:955-68. [Crossref] [PubMed]
56. Ma X, Qiu Y, Sun Y, et al. NOD2 inhibits tumorigenesis and increases chemosensitivity of hepatocellular carcinoma by targeting AMPK pathway. Cell Death Dis 2020;11:174. [Crossref] [PubMed]
57. Zhou Y, Hu L, Tang W, et al. Hepatic NOD2 promotes hepatocarcinogenesis via a RIP2-mediated proinflammatory response and a novel nuclear autophagy-mediated DNA damage mechanism. J Hematol Oncol 2021;14:9. [Crossref] [PubMed]
58. De Schutter E, Croes L, Ibrahim J, et al. GSDME and its role in cancer: From behind the scenes to the front of the stage. Int J Cancer 2021;148:2872-83. [Crossref] [PubMed]
59. Li R, Zhu S. NLRP6 inflammasome. Mol Aspects Med 2020;76:100859. [Crossref] [PubMed]
60. Ghimire L, Paudel S, Jin L, et al. The NLRP6 inflammasome in health and disease. Mucosal Immunol 2020;13:388-98. [Crossref] [PubMed]
61. Domblides C, Soubeyran I, Lartigue L, et al. Prognostic Role of Inflammasome Components in Human Colorectal Cancer. Cancers (Basel) 2020;12:3500. [Crossref] [PubMed]
62. Wei X, Xie F, Zhou X, et al. Role of pyroptosis in inflammation and cancer. Cell Mol Immunol 2022;19:971-92. [Crossref] [PubMed]
63. Huang Y, Wang JW, Huang J, et al. Pyroptosis, a target for cancer treatment? Apoptosis 2022;27:1-13. [Crossref] [PubMed]
64. Fang Y, Tian S, Pan Y, et al. Pyroptosis: A new frontier in cancer. Biomed Pharmacother 2020;121:109595. [Crossref] [PubMed]

(English Language Editor: J. Gray)