A pyroptosis-related gene signature predicts prognosis and immune microenvironment in hepatocellular carcinoma
World Journal of Surgical Oncology volume 20, Article number: 179 (2022)
Hepatocellular carcinoma (HCC) is a highly malignant tumor with a very poor prognosis. Pyroptosis is an inflammatory form of cell death and plays an important role in cancer development. The prognostic value of pyroptosis-related genes (PRGs) in HCC has not been studied extensively.
Unsupervised consensus clustering analysis was performed to identify two subtypes based on the expression profiles of prognostic PRGs in the The Cancer Genome Atlas (TCGA) database, and the differences between the two subtypes were compared. A prognostic model based on four PRGs was established by further least absolute shrinkage and selection operator (LASSO) Cox regression analysis and multivariate Cox regression analysis.
Two subtypes (clusters 1 and 2) were identified by consensus clustering based on prognostic PRGs in HCC. Survival outcomes, biological function, genomic alterations, immune cell infiltration, and immune checkpoint genes were compared between the subtypes. Cluster 2 had a worse survival outcome than cluster 1. Cluster 2 was enriched for hallmarks of cancer progression, TP53 mutation, tumor-promoting immune cells, and immune checkpoint genes, which may contribute to the poor prognosis. A prognostic risk signature that predicted the overall survival (OS) of patients was constructed and validated. Consequently, a risk score was calculated for each patient. Combined with the clinical characteristics, the risk score was found to be an independent prognostic factor for survival of HCC patients. Further analysis revealed that the risk score was closely associated with the levels of immune cell infiltration and the expression profiles of immune checkpoint genes.
Collectively, our study established a prognostic risk signature for HCC and revealed a significant correlation between pyroptosis and the HCC immune microenvironment.
Primary liver cancer is the sixth most frequently occurring malignancy and the third most common cause of cancer mortality worldwide . Hepatocellular carcinoma (HCC) is the most dominant form of liver cancer and accounts for approximately 90% of cases . Liver resection is the main curative treatment for HCC. However, most patients with HCC are diagnosed with intermediate or advanced stage disease, thus losing the opportunity for surgery. Even for patients who undergo surgical resection, postoperative recurrence or distant metastasis is common. Despite the significant progress made in diagnosis and treatment, the combination of the difficulty of early diagnosis and the ease of tumor recurrence and metastasis contribute to the poor prognosis of HCC. Therefore, prognostic biomarkers are urgently needed to help predict the outcomes for HCC patients and to guide clinical therapy.
Pyroptosis is an inflammatory form of programmed cell death, characterized by cellular swelling, many bubble-like protrusions, rupture, and release of inflammatory cytokines such as interleukin-1β (IL-1β) and and interleukin-18 (IL-18) . Pyroptosis is mainly regulated by two main pathways: the caspase-1-mediated inflammasome pathway, involving the NLRP1, NLRP3, NLRC4, AIM2, and pyrin inflammasome (the canonical inflammasome pathway), and the caspase-4/5/11-mediated inflammasome pathway (the noncanonical inflammasome pathway) . In both the canonical and noncanonical pathways, activated caspase cleaves pro-pyroptotic factor gasdermin D (GSDMD), leading to the generation of an N-terminal fragment (GSDMD-NT) . The resultant N-terminal fragment binds to phosphatidylinositol phosphates and phosphatidylserine on the inner surface of cell membranes and oligomerizes to generate pores in the lipid bilayer . These pores disrupt the osmotic balance within cells, resulting in rupture of the plasma membrane and the release of cell contents and pro-inflammatory cytokines .
Research has increasingly indicated that pyroptosis plays an important role in malignant cell transformation, growth, invasion, metastasis, and therapeutic responses. GSDMD and gasdermin E (GSDME) are two key effectors of pyroptosis. Wang et al.  reported that GSDMD expression was decreased in gastric cancers, and the downregulation of GSDMD could markedly promote the proliferation of tumors through cell cycle-related proteins accelerating S/G2 phase cell transition. GSDME can enhance immune responses to tumors such as triple-negative breast cancer by activating pyroptosis. Reduced GSDME expression is associated with decreased breast cancer survival, suggesting that GSDME might be a tumor suppressor in breast cancer . One of the main characteristics of pyroptosis is the release of immunogenic cellular content and inflammatory cytokines, both of which are closely associated with the tumor microenvironment. Wang et al.  suggested that pyroptosis-induced inflammation triggers robust antitumor immunity and can synergize with the checkpoint blockade. All findings combined indicate that the role of pyroptosis in cancers varies among different tumor cells. In HCC, a recent study revealed a mechanism of action for sorafenib, eliciting macrophage pyroptosis that enhances NK-cell effector function and ultimately effective HCC cell killing . Prognostic signatures based on PRGs were established for lung adenocarcinoma, ovarian cancer, and thyroid cancer [11,12,13]. However, the prognostic value of PRGs in HCC has not been studied extensively. In this study, we investigate the expression profiles of PRGs between normal liver tissue and HCC. Subsequently, clustering subtypes and risk models based on PRGs are established to improve prognostic risk stratification. The relationships between clustering subgroups, risk models, immune cell infiltration, and immune checkpoint genes are also analyzed to further explore the effect of PRGs on the tumor immune microenvironment (TIME).
Materials and methods
The RNA-sequencing (RNA-seq) data of 371 HCC samples and 50 adjacent normal tissues, clinical characteristics, and survival information of liver hepatocellular carcinoma (LIHC) patients were obtained from the UCSC Xena database (https://xenabrowser.net/datapages/). Three-hundred sixty-four histologically diagnosed LIHC patients with both expression profiles and survival information were included for further analysis. The RNA-seq data and survival information of 233 HCC samples in ICGC-LIRI-JP (liver cancer, RIKEN, Japan) cohort as the external validation cohort were downloaded from the International Cancer Genome Consortium (ICGC) data portal (https://icgc.org/). All datasets used in this study are publicly available.
Extraction of pyroptosis-related genes
A total of 30 pyroptosis-related genes were obtained from prior reviews [8, 12, 14,15,16,17], which are shown in Table S1. The expression differences of 30 PRGs between HCC cancer tissue and adjacent normal tissue were compared by the Wilcoxon rank-sum test. Univariate cox regression analysis was performed to evaluate the prognostic significance of the 30 PRGs.
On the basis of the expression profiles of 11 prognostic PRGs, we clustered the patients into two subgroups by the R package “ConsensusClusterPlus” . The clustering was performed using the following settings: 100 iterations, resample rate of 80%, K-means method, and Euclidean distance.
Differential gene expression and functional enrichment analysis
The differentially expressed genes (DEGs) between the two clusters were identified using the R package “limma” . Genes with a |log2 fold change (log2FC)| ≥ 1 and an adjusted p-values < 0.01 were considered as the selection criteria of DEGs. Gene ontology (GO) analysis and gene set enrichment analysis (GSEA) were carried out to functionally annotate genes that are differentially expressed in different clusters by using the R package “clusterProfiler” . The hallmark gene sets retrieved from the molecular signatures database (https://www.gsea-msigdb.org/gsea/index.jsp) were used for the GSEA.
Somatic mutation data (MuTect2 pipeline) of LIHC patients were downloaded using the “TCGAbiolinks” R package . Afterward, the somatic mutation data of different clusters were analyzed, visualized, and compared by using the “maftools” R package . TMB was defined as the number of somatic mutations in the coding region per megabase. TMB was compared using Wilcoxon rank-sum test between the two clusters.
Establishment and validation of the PRG risk signature
LASSO Cox regression analysis was performed to eliminate 11 prognostic PRGs positively correlated with each other to avoid overfitting. Subsequently, stepwise multivariate regression analysis was conducted to further screen prognostic genes by using the lowest Akaike information criterions (AIC) value. Ultimately, a PRG risk signature involved four prognostic PRGs was established. The risk score was calculated as follows: risk score = coef (gene 1) × expr (gene 1) + coef (gene 2) × expr (gene 2) + coef (gene 3) × expr (gene 3) + coef (gene 4) × expr (gene 4), where “coef” indicates the coefficient of gene derived from the multivariate Cox analysis and “expr” is the gene expression level. The LIHC patients were classified into high-risk score group and low-risk score group based on median risk score as the cutoff value in TCGA group. We used Kaplan-Meier survival analysis to compare the survival outcomes between the high-risk and low-risk groups. The receiver operating characteristic curves (ROC) were utilized to assess the specificity and sensitivity of the model as well as the accuracy at 1, 2, and 3 years by “timeROC” package . The accuracy of the PRG risk signature was then evaluated in the ICGC-LIRI-JP cohort.
Correlation of PRG risk score with tumor immune infiltration levels
The immune score was calculated using the ESTIMATE algorithm through the R “estimate package” . The immune cell infiltration was estimated by CIBERSORT algorithm . The algorithm was based on leukocyte gene signature matrix (LM22) gene signature and 1000 permutations. Only samples with a p-value < 0.05 were included to perform the subsequent analysis. Correlation analysis between PRG risk score and the immune infiltration levels in HCC was conducted with Pearson correlation coefficients.
Independence of the PRGs model
Multivariate and univariate Cox regression analysis were conducted to assess whether the risk score was an independent variable considering other clinical characteristics (age, gender, TNM stage, and grade) in the patients with HCC. A nomogram-integrated PRG risk score, age, TNM stage, and grade were established to predict the probable 1-, 2-, and 3-year survival of the patients. The concordance index (C-index) and calibration curves were applied to evaluate the concordance between predicted survival outcomes and practical outcomes.
All statistical analysis was performed using the R software (R version 4.0.2). P < 0.05 represents statistical significance.
Expression profile and prognostic value of PRGs in HCC
We first compared the expression of the 30 PRGs in hepatocellular carcinoma tissues and normal liver tissues using The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) database. In total, 22 PRGs were differentially expressed between tumor tissues and adjacent nontumorous tissues. Among these genes, 10 genes (AIM2, IL1B, IL6, NLRC4, NLRP3, NLRP6, NLRP7, TNF, GZMB, and MEFV) were significantly downregulated, and 12 genes (CASP3, CASP8, GPX4, GSDMB, GSDMC, GSDMD, DFNA5, NLRP1, NOD1, NOD2, PLCG1, and PYCARD) were significantly upregulated in HCC tissues (Fig. 1A). The correlation matrix of the expression of the 30 PRGs is shown in Fig. 1B. Most of the genes were positively correlated with each other. We then explored the prognostic value of the 30 PRGs using univariate Cox regression analysis, and 11 of them (CASP3, CASP6, CASP8, GPX4, GSDMC, DFNA5, NLRC4, NLRP6, NOD1, NOD2, and PLCG1) were identified as significantly associated with the survival outcomes of LIHC patients (Fig. 1C). NLRP6 served as a protective factor for prognosis with a hazard ratios less than 1, whereas the rest of the prognostic PRGs was all risk factors with hazard ratios greater than 1. These results demonstrated that the expression profile of PRGs in HCC significantly differed from normal liver tissue, and that PRGs might play an important role in the prognosis of LIHC patients.
LIHC classification pattern based on the prognostic PRGs
To explore the correlation between the expression of the 11 prognostic PRGs and LIHC subtypes, we performed an unsupervised consensus clustering analysis with all 364 LIHC patients. The clustering variable (k) = 2 was demonstrated to be the most appropriated selection form k = 2 to 10, and we divided the 364 LIHC patients into two clusters, namely cluster 1 (n = 185) and cluster 2 (n = 179) (Fig. 2A, Fig. SA–S1C). Principal component analysis (PCA) was conducted to test the difference between cluster 1 and cluster 2. Figure 2B shows that the two clusters were well separated from each other. The Kaplan-Meier curve indicated that patients in cluster 2 had a significantly worse OS (p < 0.001) than patients in cluster 1 (Fig. 2C). The clinicopathological features between the two clusters were then compared using the chi-squared test. Cluster 1 was associated with a low WHO grade (p < 0.05) and had a longer overall survival time (p < 0.01) compared with cluster 2 (Fig. 2D). The expression pattern of the 11 prognotic PRGs also varied significantly between the two clusters. Nine of the ten poor-prognostic genes (CASP3, CASP6, CASP8, GSDMC, DFNA5, NLRC4, NOD1, NOD2, and PLCG1) were enriched in cluster 2, while the only favorable-prognostic gene (NLRP6) was upregulated in cluster 1 (Fig. 2 D and E).
Distinct immune cell infiltration and immune checkpoint
The immune microenvironment plays crucial roles in the pathogenesis of HCC . Pyroptosis is a pro-inflammatory programmed cell death characterized by release of the immunogenic contents. We further explored the association of pyroptosis with the HCC immune microenvironment. The immune scores of the two clusters were estimated using the ESTIMATE algorithm. The immune score of cluster 2 was significantly higher than that of cluster 1 (Fig. 3B). Subsequently, the enrichment levels of 22 types of immune cells between cluster 1 and cluster 2 were compared. Cluster 1 showed higher infiltration levels of CD4 naive T cells, gamma delta T cells, activated NK cells, monocytes, and resting mast cells, whereas cluster 2 was more correlated with regulatory T cells (Tregs), M0 macrophages, M2 macrophages, resting dendritic cells, and neutrophils (Fig. 3A). The expression profiles of immune checkpoint genes, which play a key role in immune modulation, were also examined. The expression levels of PD-1, PD-L1, and CTLA-4 were all significantly upregulated in cluster 2 in comparison with cluster 1 (Fig. 3 C–E). Taken together, these results demonstrated that PRGs could indeed divided LIHC patients into two molecular subtypes with significantly different clinicopathological features, survival outcomes, and immune microenvironments.
Functional analysis of DEGs
To better illuminate the underlying molecular differences between the two clusters, we identified the DEGs between cluster 1 and cluster 2 and annotated their functions using GO and GSEA analysis. With an adjusted cutoff value of p < 0.01 and |log2(fold change)| > 1, 3935 DEGs were identified. Among them, 3885 genes were upregulated, and 150 genes were downregulated in cluster 2 compared with cluster 1 (Fig. S2A). GO analysis indicated that the DEGs were related to terms including extracellular matrix organization, extracellular structure organization, positive regulation of cell activation, phagocytosis, and immune-related binding (Fig. 4A). The GSEA analysis showed that hallmarks of cancer progression—such as epithelial mesenchymal transition (normalized enrichment score [NES] = 1.704, adjusted p-values < 0.001), G2M checkpoint (normalized enrichment score [NES] = 1.478, adjusted p-values < 0.001), angiogenesis (normalized enrichment score [NES] = 1.513, adjusted p-values < 0.001), and E2F targets (normalized enrichment score [NES] = 1.457, adjusted i-values < 0.001)—were significantly enriched in cluster 2 (Fig. 4B).
Genetic mutation landscapes of the two clusters
With rapid technological advances in genetics, the importance of the mutational profile of tumors, for both prognostic and predictive values, is being increasingly understood. We compared the tumor mutation burden (TMB) of the two clusters and found that the TMBs were not significantly different between the two clusters (Fig. S2B). Subsequently, the somatic mutation characterization of the two clusters was investigated. As shown in Fig. 4C, the top five genes with the highest mutation frequencies in cluster 1 were CTNNB1 (31%), TTN (23%), TP53 (17%), MUC16(15%), and ALB (12%), while those in cluster 2 (Fig. 4D) were TP53 (42%), TTN (26%), CTNNB1 (18%), MUC16 (17%), and LRP1B (10%). CTNNB1, DYNC2H1, TEP1, LAMA3, ATP8A1, RASA1, DUSP27, C3, and HECW2 were found to be highly mutated in cluster 1 compared with cluster 2, while TP53, TSC2, ZNF208, PTPN14, ZFHX4, and DNAH10 were found to be highly mutated in cluster 2 (Fig. 4E). These above functional enrichment analysis and mutational profile analysis indicated that the two clusters based on PRGs were also distinctively different in molecular basis and mutation characterization, which in turn affected their clinical characteristics and outcomes.
Construction and validation of a risk model based on prognostic PRGs
To construct a prognostic gene model, 11 prognostic PRGs identified by univariate Cox analysis were further filtered by LASSO Cox regression analysis and stepwise multivariate regression analysis (Fig. S3 A–C). Four prognostic PRGs—GPX4, NLRP6, NOD2, and CASP8—were identified. A prognostic model was established to evaluate the survival risk for each patient as follows: PRG risk score = 0.3570335 × expression of GPX4 + (−0.1181711) × expression of NLRP6 + 0.1128816 × expression of NOD2 + 0.3197166 × expression of CSAP8. Afterwards, 364 LIHC patients were divided into high- and low-risk groups based on the median PRG risk score. We then investigated the correlation between the cluster groups and risk groups. As shown in Fig. S3D, the high-risk group was made up of 145 (79.67%) cluster 2 patients and 37 (20.33%) cluster 1 patients, whereas the low-risk group included 137 (75.27%) patients in cluster 1 and 45 (24.73%) patients in cluster 2. This result indicated that the risk groups dichotomized based on PRG risk score were largely consistent with clustering subgroups based on the expression of the 11 prognostic PRGs. The distribution of the risk scores, OS, OS status, and expression profiles of the four PRGs-based signatures is displayed in Fig. 5A. As shown in Fig. 5 A and B, patients in the high-risk group exhibited shorter OS times (p < 0.001) and higher mortality. The heatmap (Fig. 5A) and boxplot (Fig. S4A) indicated that risky PRGs, including GPX4, NOD2, and SCAF11, were highly enriched in the high-risk group, whereas the expression levels of the protective gene NLRP6 were significantly enriched in the low-risk group. We also performed ROC curve analysis to evaluate the predictive accuracy of our risk model. In the TCGC database, the AUCs for the PRG signatures at the 1-, 2-, and 3-year OS were 0.749, 0.669, and 0.648, respectively (Fig. 5C). We then verified our prognostic model using an independent verification dataset from the ICGC database. The risk score for each of the 232 LIHC patients from the ICGC was calculated using the same formula in the TCGA dataset, and patients were divided into high- and low-risk groups based on the median risk score. A similar result was observed between the ICGC and TCGA cohorts, where the high-risk group exhibited shorter OS times (p = 0.0018) and higher mortality (Fig. 5 D and E). The three risky PRGs and one protective PRG were also significantly enriched in the high-risk group and low-risk group, respectively (Fig. 5D and Fig. S4B). In the ICGC cohort, the AUCs of the PRG signature corresponding to 1, 2, and 3 years of survival were 0.631, 0.714, and 0.736, respectively (Fig. 5F). Taken together, our prognostic model showed a relatively high accuracy in predicting the prognosis of LIHC patients.
The risk signature correlated with the clinical characteristics and the immune landscape of HCC
We evaluated the relationship between the PRG risk scores and clinical characteristics. As shown in the heatmap, the difference in T stage (p < 0.05), TNM stage (p < 0.05), grade (p < 0.001), and overall survival (p < 0.001) between the high- and low-risk groups was significant (Fig. 6A). Tumor-infiltrating immune cells (TIICs) play essential roles in cancer development and are closely associated with clinical outcomes. We assessed the proportions of different TIICs and explored the correlation between the PRG risk score and TIICs. Figure 6 B–D demonstrates that the PRG risk score was positively correlated with the infiltration levels of Tregs and M2 macrophages but negatively correlated with the infiltration levels of resting mast cells. Moreover, we investigated the relationship between the PRG risk score and the expression of immune checkpoints. Figure 7 A–F shows that the PRG risk score had a significant positive correlation with the expression of common immune checkpoints and their ligands, namely anti-programmed cell death protein 1 (PD-1) and its ligand programmed cell death ligand 1/2 (PD-L1/2) and anti-cytotoxic T-lymphocyte-associated antigen-4 (CTLA-4) and its ligand CD86/CD80. The expression of novel immune checkpoints, including lymphocyte activation gene-3 (LAG-3), T-cell immunoglobulin and mucin-domain containing-3 (TIM-3), T-cell immunoglobulin and ITIM domain (TIGIT), B7H3, Galectin-9, and V-domain Ig suppressor of T cell activation (VISTA), was also significantly associated with the PRG risk score (Fig. 8G–L). Taken together, these results suggested that the PRG risk score is significantly associated with the clinical characteristics, the landscape of immune cell infiltrations, and the expression levels of immune checkpoint genes in LIHC patients.
The risk signature is an independent prognostic factor for patients with HCC
Univariate and multivariate Cox regression analyses were carried out to assess whether the risk signature was an independent prognostic predictor for OS in the TCGA cohort. As described in Fig. 8A, the results of the univariate Cox regression analysis showed that T stage (HR: 2.462, 95% CI :1.695–3.577, p < 0.001), M stage (HR: 3.740, 95% CI: 1.183–11.819, p = 0.025), TNM stage (HR: 2.440, 95% CI: 1.682–3.540, p < 0.001), and the PRG risk score (HR: 2.712, 95% CI:1.863–3.948, p < 0.001) were significantly associated with the prognosis of LIHC patients. After multivariate Cox regression analysis, only the PRG risk score (HR: 2.637, 95% CI: 1.771–3.925, p < 0.001) remained an independent predictor for the prognosis of LIHC patients. To facilitate the use of our PRGs prognostic model, we established a nomogram comprising the PRG risk score and clinical characteristics for predicting LIHC prognosis. By comparison with clinical characteristics, the PRG risk score showed predominant predictive ability in the nomogram (Fig. 8B). A calibration plot demonstrated ideal consistency compared with the ideal model, indicating that the nomogram has stability for predicting LIHC patient prognosis in clinical practice (Fig. 8C). The C-index of the established nomogram for OS prediction was 0.695 (Fig. 8D). These results collectively verified that our PRG risk signature could reliably serve as an independent prognostic factor for LIHC patients.
Pyroptosis is a novel form of programmed cell death triggered by certain inflammasomes. Studies have found that pyroptosis plays a dual role in the proliferation, invasion, and metastasis of tumors. On the one hand, pyroptosis can promote inflammatory cell death of cancer and inhibit the proliferation and migration of cancer cells. On the other hand, inflammatory mediators IL-1 and IL-18, released by the activation of pyroptosis, can promote the progression of various cancers . The expression level of GSDME was significantly increased in esophageal squamous cell carcinoma (ESCC) tissues than in normal esophageal tissues and was significantly correlated with a better prognosis in patients with ESCC. In GSDME high-expressed ESCC cell lines, the serine/threonine protein kinase PLK1 inhibitor BI2536 combined with cisplatin can increase chemosensitivity by inducing GSDMD-mediated pyroptosis . The expression of GSDMD was significantly upregulated in non-small cell lung cancer (NSCLC), and higher GSDMD expression is associated with invasive features, including more advanced tumor-node-metastasis stages and larger tumor sizes . Apoptosis, another form of caspase-mediated cell death, has been found to be extensively involved in the tumorigenesis and development of HCC. Inhibition of HIF1A-AS1, a long noncoding RNA (lncRNA) overexpressing in HCC tissues, could enhanced starvation-induced apoptosis in HCC cell . A prognostic model involved apoptosis-related genes with impressive prognostic predictive power was also established in HCC . Furthermore, numerous prognostic genes and prognostic signature for HCC have been reported by many studies [32,33,34,35] and further evaluated by some meta-analysis [36, 37]. However, the relationship between pyroptosis and HCC is not fully understood. Most studies regarding pyroptosis in human HCC have found pyroptosis plays an antitumor role in HCC. For example, Wei et al.  reported that the expression of NLRP3 was either completely lost or significantly downregulated in human HCC, and that the deficiency correlated significantly with poor pathological differentiation and advanced stages, indicating that the NLRP3 inflammasome was involved in the progression of HCC. One distinct characteristic which separates pyroptosis from other kinds of cell death is the activation of caspase-1. Chu et al.  observed a significantly decreased expression of caspase-1 in HCC tissues from patients and found that activation of caspase-1-dependent pyroptosis shows therapeutic potential against HCC. In this study, we first studied the mRNA levels of 30 PRGs in HCC samples and adjacent normal tissues and found that 22 PRGs were differentially expressed. In detail, the expression levels of AIM2, IL1B, IL6, NLRC4, NLRP3, NLRP6, NLRP7, TNF, GZMB, and MEFV were decreased, while the expression levels of CASP3, CASP8, GPX4, GSDMB, GSDMC, GSDMD, DFNA5, NLRP1, NOD1, NOD2, PLCG1, and PYCARD were increased in HCC samples compared with normal tissues. These results indicated that pyroptosis might be extensively involved in HCC tumorigenesis. Univariate Cox regression analysis was performed to assess the prognostic value of the 30 PRGs and found that 11 PRGs were significantly associated with overall survival. Based on the expression profiles of the 11 prognostic genes, LIHC patients were classified into two clusters with distinct survival outcomes and clinicopathological features. Patients in clusters 2 were more likely to have a shorter survival time, higher WHO grade, and higher immune score. Further analysis revealed that cluster 2 had increased levels of immune checkpoint genes. The CIBERSORT algorithm was used to calculate the proportion of different types of tumor-infiltrating immune cells. The result showed that compared with cluster 1, immune cells that promote tumor proliferation, therapeutic resistance, and metastasis, such as M2 macrophages, regulatory T cells, and neutrophils, were enriched in cluster 2 [40,41,42]. These results suggest a significant correlation between PRGs and tumor-immune infiltration.
Certain gene mutations in HCC are closely related to the survival outcomes and vary significantly between different subgroups of HCC. TERT promoter mutations, the tumor protein p53 (TP53), and WNT pathway oncogene catenin beta 1 (CTNNB1) are the most frequent somatic genetic alterations in HCC . Several studies have shown that HCCs with mutations in CTNNB1 display a particular phenotype with well-differentiated tumors and better prognostic outcomes. HCCs with TP53 mutation and an absence of CTNNB1 mutation display aggressive tumor characteristics and worse prognosis . We found that our subgroups based on PRGs had significantly different mutational profiles. In cluster 1, which had a longer overall survival time, 31% of patients had mutations in CTNNB1, and 17% had mutations in TP53. In cluster 2, which had worse prognosis, 42% of patients had mutations in TP53, and 18% of patients had mutations in CTNNB1.
Our study constructed a prognostic gene model based on four prognostic PRGs (GPX4, NLRP6, NOD2, and CASP8) and generated a PRG risk score for each LIHC patient. Glutathione peroxidase 4 (GPX4) is an essential regulator of ferroptosis, and an antitumor effect of GPX4 inhibition-induced ferroptosis has been found in a variety of cancers, such as breast, kidney, and colorectal cancers [45,46,47]. Recent studies have suggested that GPX4 also functions as an important gateway to pyroptosis. Kang et al.  found that GPX4 is upregulated in innate immune cells to counter-regulate GSDMD-N-mediated pyroptotic cell death, thereby preventing lethal systemic inflammation. Guerriero E. reported a significantly overexpression of GPX4 in HCC compared with nontumor tissues, and this overexpression was associated with an increased malignancy grade . However, another study found that GPX4 suppressed the formation and progression of HCC by inhibition of angiogenesis and tumor cell proliferation as well as by immune-mediated mechanisms . In our study, multivariate Cox regression analysis showed that GPX4 was a risky gene with an HR larger than 1 (HR = 1.45, 95% CI = 1.14–1.85), and GPX4 expression levels were significantly upregulated in the high-risk group compared with those in the low-risk group. NLR family pyrin domain containing 6 (NLRP6) is a member of the NLR (nucleotide-oligomerization domain-like receptor) family that patrols the cytosolic compartment of cells to detect pathogen- and damage-associated molecular patterns . Wang et al.  demonstrated that NLRP6 was downregulated in approximately 75% of primary gastric cancer cases and functioned as a negative regulator of gastric cancer. In another study, NLRP6 was found to perform essential functions in the regulation of tissue repair necessary for protection against chemically induced injury . Several studies have suggested a protective role of NLRP6 in a variety of liver diseases, such as liver cirrhosis, liver injury, acute liver failure, and alcoholic hepatitis [54,55,56,57]. However, the relationships between NLRP6 and HCC remain largely unknown. In our study, we found that NLRP6 was one of the pyroptosis-related prognostic biomarkers in HCC and was likely to serve as a protective factor against HCC. Further in vivo and in vitro studies should be performed to elucidate the role of NLRP6 in HCC. Nucleotide-binding oligomerization domain 2 (NOD2) is a member of the family of pattern recognition receptors (PRRs) that can initiate potent immune response against pathogens . Ma et al.  reported that NOD2 deficiency promoted hepatocarcinogenesis, while overexpression of NOD2 in HCC cells inhibited tumorigenesis and reversed resistance to chemotherapy. However, a recent study by Zhou et al.  suggested an opposite role for NOD2 in the development of HCC. Zhou et al. showed that NOD2 was overexpressed in HCC samples and closely correlated with poor prognosis of LIHC patients. Loss of hepatic NOD2 attenuated the tumorigenesis of DEN/CCl4-induced HCC in hepatocyte-specific Nod2-knockout mice. We found that NOD2 was a cancer-promoting gene, as it was enriched in the high-risk group and associated with poor prognosis among LIHC patients. Given these contradictory results, further studies will be needed to understand the relationship between NOD2 and the development of HCC. Caspase-8 (CASP8) has long been considered as an initiator of the extrinsic apoptotic pathway . Recent studies have suggested that CASP8 also cleaves and activates GSDMD, inducing pyroptosis in response to activation of the extrinsic apoptosis signaling pathway . An increasing number of studies have confirmed that CASP8 is associated with tumor growth and invasion, angiogenesis, metastasis, therapeutic resistance, and poor clinical outcomes [63, 64]. In our research, we found that CASP8 was upregulated in the high-risk group and was prognostic for poor survival.
The prognostic value of the PRG risk model was evaluated in patients with LIHC and validated in the external ICGC cohort. The risk score obtained from the risk signatures effectively classified the patients with LIHC into high- and low-risk groups. The OS of the patients in the high-risk group was shorter than that of the patients in the low-risk group in the TCGA cohort. Consistent results were also obtained in an independent ICGC validation cohort. The risk groups classified by the risk score were largely consistent with the classification based on the expression profile of prognostic PRGs. These results indicated that the risk signature based on PRGs may help to improve the clinical assessment of patients, optimize medical intervention, and identify potential therapeutic targets for different subgroups.
Accumulating evidence has revealed that immune infiltrates in tumor microenvironment play a significant role in the prognosis of HCC. SNRPC, PKB, and DCK have been discovered to significantly correlate with patient outcomes and immune cell infiltration in hepatocellular carcinoma [65,66,67]. We further explored the relationship between the PRG risk score and tumor-infiltrating immune cells. The PRG risk score was positively correlated with Tregs and M2 macrophages but negatively correlated with the infiltration levels of resting mast cells. Notably, the expression levels of several immune checkpoint-related genes, including PD-1, PD-L1, PD-L2, CTLA-4, CD86, CD80, LAG-3, TIM-3, TIGIT, B7H3, Galectin-9, and VISTA, were all significantly enriched in the high-risk group and correlated with the risk score. Taken together, we speculated that overexpression of immune checkpoint and tumor-promoting immune cells impaired the efficacy of the antitumor immune response, thus leading to the unfavorable prognosis in the high-risk group. Immunotherapy is a promising treatment option for HCC, and its efficacy is significantly influenced by the HCC immune microenvironment . Given its close association with the HCC immune microenvironment, our PRG risk signature may serve as a useful tool to stratify LIHC patients into different immune subtypes and predict the sensitivity to immunotherapy.
Univariate and multivariate Cox regression analysis showed that the risk score was an independent prognostic factor for LIHC patients. We integrated risk score, TNM stage, grade, age, and gender to construct a nomogram, and we found a good agreement between the actual OS and predicted OS at 1, 2, and 3 years. These results indicated that the PRG signature combined with commonly available clinical characteristics might be a promising prognostic tool for LIHC patients.
Our study has a few limitations. First, our study is based on the TCGA public database and only verified by the ICGC public database. Second, given the high heterogeneity of HCC, relatively few tumor samples are included in our research. Furthermore, only protein-coding genes were explored in the present study. However, it has been established that noncoding RNAs, such as lncRNA and microRNAs, also play essential role in tumorigenesis, progression, and immunity of HCC. For example, Xue et al. reported that lncRNA ZEB1-AS1 inhibited HCC progression through miR-23c . In another study, a prognostic signature and a nomogram integrating three lncRNAs were constructed for HCC patients . Further studies that integrated multicenter datasets and different types of RNA genes are needed to test the results. Thirdly, our research should be further explored by both in vitro and in vivo experiments to elucidate the precise roles of pyroptosis in HCC.
In conclusion, our study showed that pyroptosis is heavily involved in HCC. Most of the PRGs were differently expressed between normal and HCC tissues. Further comprehensive bioinformatics analysis identified HCC subtypes based on PRG signatures. The survival outcomes, clinical characteristics, genomic alterations, and immune microenvironment differed significantly between the subgroups. The risk score generated by the PRG signature is closely correlated with the infiltration levels of immune cells and the expression profile of immune checkpoint genes.
Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49.
Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2021;7(1):6.
Yu P, Zhang X, Liu N, Tang L, Peng C, Chen X. Pyroptosis: mechanisms and diseases. Signal Transduct Target Ther. 2021;6(1):128.
Kayagaki N, Stowe IB, Lee BL, et al. Caspase-11 cleaves gasdermin D for non-canonical inflammasome signalling. Nature. 2015;526(7575):666–71.
Shi J, Zhao Y, Wang K, et al. Cleavage of GSDMD by inflammatory caspases determines pyroptotic cell death. Nature. 2015;526(7575):660–5.
Liu X, Zhang Z, Ruan J, et al. Inflammasome-activated gasdermin D causes pyroptosis by forming membrane pores. Nature. 2016;535(7610):153–8.
Wang WJ, Chen D, Jiang MZ, et al. Downregulation of gasdermin D promotes gastric cancer proliferation by regulating cell cycle-related proteins. J Dig Dis. 2018;19(2):74–83.
Zhang Z, Zhang Y, Xia S, et al. Gasdermin E suppresses tumour growth by activating anti-tumour immunity. Nature. 2020;579(7799):415–20.
Wang Q, Wang Y, Ding J, et al. A bioorthogonal system reveals antitumour immune function of pyroptosis. Nature. 2020;579(7799):421–6.
Hage C, Hoves S, Strauss L, et al. Sorafenib induces pyroptosis in macrophages and triggers natural killer cell-mediated cytotoxicity against hepatocellular carcinoma. Hepatology. 2019;70(4):1280–97.
Lin W, Chen Y, Wu B, Chen Y, Li Z. Identification of the pyroptosisrelated prognostic gene signature and the associated regulation axis in lung adenocarcinoma. Cell Death Discov. 2021;7(1):161.
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(1):71.
Wu P, Shi J, Sun W, Zhang H. Identification and validation of a pyroptosis-related prognostic signature for thyroid cancer. Cancer Cell Int. 2021;21(1):523.
Kovacs SB, Miao EA. Gasdermins: effectors of pyroptosis. Trends Cell Biol. 2017;27(9):673–84.
Man SM, Karki R, Kanneganti TD. Molecular mechanisms and functions of pyroptosis, inflammatory caspases and inflammasomes in infectious diseases. Immunol Rev. 2017;277(1):61–75.
Fan R, Sui J, Dong X, Jing B, Gao Z. Wedelolactone alleviates acute pancreatitis and associated lung injury via GPX4 mediated suppression of pyroptosis and ferroptosis. Free Radic Biol Med. 2021;173:29–40.
Zhou Z, He H, Wang K, et al. Granzyme A from cytotoxic lymphocytes cleaves GSDMB to trigger pyroptosis in target cells. Science. 2020;368(6494):eaaz7548.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
Ritchie ME, Phipson B, Wu D, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44(8):e71.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.
Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32(30):5381–97.
Yoshihara K, Shahmoradgoli M, Martinez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Llovet JM, Castet F, Heikenwalder M, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol. 2022;19(3):151–72.
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.
Wu M, Wang Y, Yang D, et al. A PLK1 kinase inhibitor enhances the chemosensitivity of cisplatin by inducing pyroptosis in oesophageal squamous cell carcinoma. EBioMedicine. 2019;41:244–55.
Gao J, Qiu X, Xi G, et al. Downregulation of GSDMD attenuates tumor proliferation via the intrinsic mitochondrial apoptotic pathway and inhibition of EGFR/Akt signaling and predicts a good prognosis in nonsmall cell lung cancer. Oncol Rep. 2018;40(4):1971–84.
Hong F, Gao Y, Li Y, Zheng L, Xu F, Li X. Inhibition of HIF1A-AS1 promoted starvation-induced hepatocellular carcinoma cell apoptosis by reducing HIF-1alpha/mTOR-mediated autophagy. World J Surg Oncol. 2020;18(1):113.
Liu R, Wang G, Zhang C, Bai D. A prognostic model for hepatocellular carcinoma based on apoptosis-related genes. World J Surg Oncol. 2021;19(1):70.
Yang H, Huo J, Li X. Identification and validation of a five-gene prognostic signature for hepatocellular carcinoma. World J Surg Oncol. 2021;19(1):90.
Yu M, Xu W, Jie Y, et al. Identification and validation of three core genes in p53 signaling pathway in hepatitis B virus-related hepatocellular carcinoma. World J Surg Oncol. 2021;19(1):66.
Sun G, Sun K, Shen C. Human nuclear receptors (NRs) genes have prognostic significance in hepatocellular carcinoma patients. World J Surg Oncol. 2021;19(1):137.
Wu M, Liu Z, Li X, Zhang A, Lin D, Li N. Analysis of potential key genes in very early hepatocellular carcinoma. World J Surg Oncol. 2019;17(1):77.
Chen J, Pan W, Chen Y, Wen L, Tu J, Liu K. Relationship of ALDH2 rs671 and CYP2E1 rs2031920 with hepatocellular carcinoma susceptibility in East Asians: a meta-analysis. World J Surg Oncol. 2020;18(1):21.
Quan Y, Yang J, Qin T, Hu Y. Associations between twelve common gene polymorphisms and susceptibility to hepatocellular carcinoma: evidence from a meta-analysis. World J Surg Oncol. 2019;17(1):216.
Wei Q, Zhu R, Zhu J, Zhao R, Li M. E2-induced activation of the NLRP3 inflammasome triggers pyroptosis and inhibits autophagy in HCC cells. Oncol Res. 2019;27(7):827–34.
Chu Q, Jiang Y, Zhang W, et al. Pyroptosis is involved in the pathogenesis of human hepatocellular carcinoma. Oncotarget. 2016;7(51):84658–65.
Yeung OW, Lo CM, Ling CC, et al. Alternatively activated (M2) macrophages promote tumour growth and invasiveness in hepatocellular carcinoma. J Hepatol. 2015;62(3):607–16.
Gao Q, Qiu SJ, Fan J, et al. Intratumoral balance of regulatory and cytotoxic T cells is associated with prognosis of hepatocellular carcinoma after resection. J Clin Oncol. 2007;25(18):2586–93.
Zhou SL, Zhou ZJ, Hu ZQ, et al. Tumor-associated neutrophils recruit macrophages and T-regulatory cells to promote progression of hepatocellular carcinoma and resistance to sorafenib. Gastroenterology. 2016;150(7):1646–1658 e1617.
Zucman-Rossi J, Villanueva A, Nault JC, Llovet JM. Genetic landscape and biomarkers of hepatocellular carcinoma. Gastroenterology. 2015;149(5):1226–1239 e1224.
Nishida N, Nishimura T, Kaido T, et al. Molecular scoring of hepatocellular carcinoma for predicting metastatic recurrence and requirements of systemic chemotherapy. Cancers (Basel). 2018;10(10):367.
Ding Y, Chen X, Liu C, et al. Identification of a small molecule as inducer of ferroptosis and apoptosis through ubiquitination of GPX4 in triple negative breast cancer cells. J Hematol Oncol. 2021;14(1):19.
Zou Y, Palte MJ, Deik AA, et al. A GPX4-dependent cancer cell state underlies the clear-cell morphology and confers sensitivity to ferroptosis. Nat Commun. 2019;10(1):1617.
Sui X, Zhang R, Liu S, et al. RSL3 drives ferroptosis through GPX4 inactivation and ROS production in colorectal cancer. Front Pharmacol. 2018;9:1371.
Kang R, Zeng L, Zhu S, et al. Lipid peroxidation drives gasdermin D-mediated pyroptosis in lethal polymicrobial sepsis. Cell Host Microbe. 2018;24(1):97–108 e104.
Guerriero E, Capone F, Accardo M, et al. GPX4 and GPX7 over-expression in human hepatocellular carcinoma tissues. Eur J Histochem. 2015;59(4):2540.
Rohr-Udilova N, Bauer E, Timelthaler G, et al. Impact of glutathione peroxidase 4 on cell proliferation, angiogenesis and cytokine production in hepatocellular carcinoma. Oncotarget. 2018;9(11):10054–68.
Li R, Zhu S. NLRP6 inflammasome. Mol Aspects Med. 2020;76:100859.
Wang X, Wu X, Wang Q, Zhang Y, Wang C, Chen J. NLRP6 suppresses gastric cancer growth via GRP78 ubiquitination. Exp Cell Res. 2020;395(1):112177.
Normand S, Delanoye-Crespin A, Bressenot A, et al. Nod-like receptor pyrin domain-containing protein 6 (NLRP6) controls epithelial self-renewal and colorectal carcinogenesis upon injury. Proc Natl Acad Sci U S A. 2011;108(23):9601–6.
Zhu Y, Ni T, Deng W, et al. Effects of NLRP6 on the proliferation and activation of human hepatic stellate cells. Exp Cell Res. 2018;370(2):383–8.
Li M, Chen Y, Shi J, et al. NLRP6 deficiency aggravates liver injury after allogeneic hematopoietic stem cell transplantation. Int Immunopharmacol. 2019;74:105740.
Schneider KM, Elfers C, Ghallab A, et al. Intestinal dysbiosis amplifies acetaminophen-induced acute liver injury. Cell Mol Gastroenterol Hepatol. 2021;11(4):909–33.
Ji X, Li L, Lu P, Li X, Tian D, Liu M. NLRP6 exerts a protective role via NF-kB with involvement of CCL20 in a mouse model of alcoholic hepatitis. Biochem Biophys Res Commun. 2020;528(3):485–92.
Caruso R, Warner N, Inohara N, Nunez G. NOD1 and NOD2: signaling, host defense, and inflammatory disease. Immunity. 2014;41(6):898–908.
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(3):174.
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(1):9.
Mandal R, Barron JC, Kostova I, Becker S, Strebhardt K. Caspase-8: the double-edged sword. Biochim Biophys Acta Rev Cancer. 2020;1873(2):188357.
Schwarzer R, Laurien L, Pasparakis M. New insights into the regulation of apoptosis, necroptosis, and pyroptosis by receptor interacting protein kinase 1 and caspase-8. Curr Opin Cell Biol. 2020;63:186–93.
Muller I, Strozyk E, Schindler S, et al. Cancer cells employ nuclear Caspase-8 to overcome the p53-dependent G2/M checkpoint through cleavage of USP28. Mol Cell. 2020;77(5):970–984 e977.
Fianco G, Mongiardi MP, Levi A, et al. Caspase-8 contributes to angiogenesis and chemotherapy resistance in glioblastoma. Elife. 2017;6:e22593.
Cai J, Zhou M, Xu J. N6-methyladenosine (m6A) RNA methylation regulator SNRPC is a prognostic biomarker and is correlated with immunotherapy in hepatocellular carcinoma. World J Surg Oncol. 2021;19(1):241.
Mu W, Xie Y, Li J, et al. High expression of PDZ-binding kinase is correlated with poor prognosis and immune infiltrates in hepatocellular carcinoma. World J Surg Oncol. 2022;20(1):22.
Song D, Wang Y, Zhu K, et al. DCK is a promising prognostic biomarker and correlated with immune infiltrates in hepatocellular carcinoma. World J Surg Oncol. 2020;18(1):176.
Sangro B, Sarobe P, Hervas-Stubbs S, Melero I. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2021;18(8):525–43.
Xue S, Lu F, Sun C, Zhao J, Zhen H, Li X. LncRNA ZEB1-AS1 regulates hepatocellular carcinoma progression by targeting miR-23c. World J Surg Oncol. 2021;19(1):121.
Zhang C, Bao T, Ke Y, et al. Integrated analysis of ceRNA network reveals potential prognostic Hint1-related lncRNAs involved in hepatocellular carcinoma progression. World J Surg Oncol. 2022;20(1):67.
We would like to thank the other members of our research team for their help and TCGA and ICGC databases for the availability of the data. We thank Catherine Perfect, MA (Cantab), from Liwen Bianji (Edanz) (www.liwenbianji.cn), for editing the English text of a draft of this manuscript.
This work was supported by the National Natural Science Foundation of China (81971479), Major Research Project of Science Technology Department of Zhejiang Province (2019C03048), and Foundation Project for Medical Science and Technology of Zhejiang province (2019ZH022).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Consensus clusters by PRGs in TCGA cohort. Figure S2. DEGs and TMB scores of the clusters 1 and 2. Figure S3. Screening of four PRGs signature genes. Figure S4. The expressions of four prognostic PRGs in high- and low-risk groups. Table S1. Names of 30 pyroptosis-related genes.
About this article
Cite this article
Jin, Y., Pu, X., Ping, D. et al. A pyroptosis-related gene signature predicts prognosis and immune microenvironment in hepatocellular carcinoma. World J Surg Onc 20, 179 (2022). https://doi.org/10.1186/s12957-022-02617-y
- Immune checkpoint gene
- Immune infiltrates