The RNA sequence data and related clinical information of 374 liver cancer patients (TCGA-LIHC) were acquired from the TCGA website (https://portal.gdc.cancer.gov/repository). The gene expression data were normalized by scale method using the “limma” package (Yuan C. et al., 2021). After excluding the missing clinical information of patients, 370 HCC patients were randomly separated into the training and the test groups by the “caret” package. Besides, transcriptomics data with clinical features of 231 HCC patients (LIRI-JP) were downloaded from the ICGA database (https://dcc.icgc.org/projects/LIRI-JP). These HCC patients were HBV or HCV carriers from Japan. The gene read count values of these patients were also normalized, and both the TCGA and ICGC data were public.

A total of thirty-three PRGs were retrieved from the previous literature and are listed in Supplementary Table S1 (Man and Kanneganti, 2015; Tang et al., 2020; Ye J. et al., 2021). The differentially expressed genes (DEGs) were analyzed by the “limma” package in the R software and visualized by the heatmap and volcano plot with adjusted p-value < 0.05 in the TCGA cohort. In addition, a protein-protein interaction network for the PRGs was generated using the software GENEMANIA (http://genemania.org/).

The Univariate Cox analysis was employed to screen the PRGs with the prognostic value, and the cut-off p-value was set at 0.05, and 17 survival-related genes were used for further study. The prognostic model was established to minimize overfitting using the LASSO-penalized Cox regression analysis via the “glmnet” R package (Du et al., 2021). Eventually, the three genes and their coefficients were retained, and the minimum criteria determined the penalty parameter (λ). The risk score was obtained using the formula: risk score = (βA × Gene A expression) + (βB × Gene B expression) ⋯ + (βN × Gene N expression), in which β represents regression coefficient (Yang et al., 2021). The HCC patients were separated into the high- and low-risk groups based on the median risk score. Then, the principal component analysis (PCA) for the two risk groups was constructed using the “limma” and “scatterplot3d” R packages in terms of gene expressions in the prognostic model (Yuan M. et al., 2021). The survival analysis between the two risk groups was carried out using the “survminer” R package, and the ROC curve analysis was performed via the “survival” and “timeROC” R packages (Fang et al., 2021). Besides, the univariate and multivariate Cox regression was utilized to determine the independent prognostic value of the 3-gene signature.

To validate the efficiency of our model, the patients in the TCGA internal test cohort or ICGC external validation cohort (LIRI-JP) were applied. The mRNA levels of PRGs were normalized according to the “scale” function, and the risk score was calculated using the same formula applied in the TCGA training cohort. The TCGA test or ICGC cohort patients were separated into the low- and high-risk groups using the median risk score obtained from the TCGA training cohort.

Based on the risk score and different clinical features (gender, age, histologic grade, and pathological stage), a nomogram model was established to predict the survival years for the HCC patients using the “rms” and “survival” packages (Dai et al., 2021).

The GSEA analysis was carried out between the two risk groups using the GSEA software 4.0.1 to identify the differential KEGG pathways. The normalized enrichment scores and nominal p-values were determined for analyzing the enrichment levels and statistical significance. In addition, the infiltrating scores of the 16 immune cells and activities of 13 immune-related pathways were analyzed using ssGSEA provided in the “GSVA” R package (Zhao et al., 2021).

The NCI-60 dataset covering nine cancer categories was available on the CellMiner homepage (https://discover.nci.nih.gov/cellminer) (Shankavaram et al., 2007; Shankavaram et al., 2009). In addition, Pearson correlation analysis was performed to indicate the relevance between the independent prognostic PRGs and drug sensitivity. Drugs used in this sensitivity analysis are those approved by the FDA or those in clinical tests.

The FVB mice (8 weeks) were purchased from Westlake University (Zhejiang, China). Plasmids for injection were harvested using the Endotoxin Free Maxiprep kit (Qiagen) and delivered to the FVB mice as a mixture using the hydrodynamic tail vein injection for HCC formation as described before (Song et al., 2017). Dosages of these plasmids were: px330-U6-sgP53 (mouse) 20 μg, pT3-EF1α-c-Myc (human) 5 μg, pCMV-sleeping beauty transposase 2 μg. The daily abdominal palpation of mice was observed, and the mice were then sacrificed when they suffered from high burdens of liver tumors. The mice were raised according to the protocols approved by the Institutional Animal Care and Use Committee of Westlake University. The researcher followed the standard biosecurity and institutional safety procedures in this study.

The human normal liver cell line HL-7702 and HCC cell lines (SK-Hep1 and Huh7) were cultured in a DMEM medium containing 10% FBS in a 5% CO₂ incubator at 37°C. In addition, huh7 cells were treated by Lenvatinib (HY-10981, MCE) with different concentrations (0, 0.5, 1, 5, 10 μM) for 48 h based on the previous study (Jin et al., 2021).

RNA extraction and quantitative real-time reverse transcriptase-polymerase chain reaction (qRT-PCR).

The total RNA in the liver tissues or liver cell lines was extracted using the TRIZOL (Invitrogen) reagent. Then, the cDNA was obtained by reverse transcription using the cDNA Synthesis Mix (Novoprotein, E047-01B) and analyzed using quantitative PCR (Novoprotein, E096-01B). Finally, the mRNA expression levels were normalized with the ACTB gene. All the primers used in this study are listed in Supplementary Table S2.

The total proteins in Huh7 cells treated with Lenvatinib were extracted by RIPA buffer containing 1% PMSF (Beyotime) and quantified by the BCA method (Beyotime). Twenty microgram proteins were used for western blotting against GSDME (1:2000, ab215191, Abcam) and Hsp90 (1:2000, 610419, BD Bioscience).

The one-way ANOVA was adopted to compare the DEGs between the HCC and normal liver tissues. The Mann—Whitney test using the BH method adjusted p-value was adopted to measure the ssGSEA scores. The entire statistical analyses were achieved using the R software 4.0.1. The statistical significance was considered as a P -value less than 0.05 if not otherwise specified.

We performed the differential expression analysis of 33 PRGs involved between the normal liver and HCC samples in the TCGA database (Figure 1A). The 33 PRGs were chosen based on their roles in pyroptosis according to previous studies (Karki and Kanneganti, 2019; Ye Y. et al., 2021; Ju et al., 2021). The results were presented in the heatmaps and volcano plots (Figures 1B,C). We identified 26 differential genes (p < 0.05), among which the three genes were significantly downregulated in the HCC samples, including IL6, IL-1β, and NLRP3. Conversely, another 23 genes were significantly upregulated, including GSDMC, PLCG1, PYCARD, GSDME, NLRP1, GSDMB, GSDMD, CASP8, NOD1, CASP8, NOD1, CASP3, GSDMA, PJVK, TIRAP, NOD2, AIM2, GPX4, NLRP7, CASP9, CASP6, PRKACA, SCAF11, CASP4, and NLRP6. The interactions within these PRGs were identified using GeneMANIA and constituted three interconnected subnetworks: the NLR superfamily regulatory network, the caspase family regulatory network, and the GSDM family regulatory network (Figure 1D). Finally, the correlation for these PRGs was analyzed and indicated in Figure 1E.

We first analyzed the prognostic values of PRGs in the TCGA cohort with the Univariate Cox regression analysis. The high expression of CASP1, CASP3, CASP5, CASP6, CASP8, GPX4, GSDMA, GSDME, NLRC4, NLRP3, NLRP7, NOD1, NOD2, PLCG1, and SCAF11 correlated with the poor survival of the HCC patients, as indicated by the Hazard ratio (HR) > 1 (Figure 2A). Then, we performed the consensus clustering analysis to investigate the relationship between these prognostic genes and HCC subtypes. According to the CDF value, we classified the 370 HCC patients into two clusters (k = 2, Figures 2B–D), and we found that the patients from cluster 1 tended to survive longer than the patients from cluster 2 (Figure 2E), implying a significant prognostic value of these PRGs. Furthermore, the two clusters did not differ in the clinical parameters such as stage, grade, gender, age, and TMN (Figure 2F).

We split the samples in the TCGA-LIHC dataset into two equal cohorts: training cohort and test cohort at random. We found no significant difference between the training and test cohorts from TCGA-LIHC in the leading clinical indicators (Supplementary Table S3). We performed lasso regression analysis using 17 prognostic genes in the training cohort to build the prognostic model. According to the minimum criteria, a risk model consisting of GSDME, GPX4, and SCAF11 was built (Figures 3A,B). The risk score was obtained by the formula: risk score = (0.0182*GSDME exp.) + (0.0005*GPX4 exp.) + (0.0188*SCAF11 exp.)

We separated the patients in the training cohort into the high- and low-risk groups based on the median risk score (Figure 3C). The high-risk group had more deaths and shorter survival years (Figure 3D). The heatmap analysis indicated that high-risk patients had increased expression levels of three risk genes (Figure 3E). In addition, the Kaplan-Meier curve indicated that the high-risk patients had worse overall survival (OS) than low-risk patients (p < 0.001). The subsequent ROC analysis demonstrated that this 3-gene risk model could robustly evaluate and predict the survival of the HCC patients (AUC = 0.763, Figure 3G). Furthermore, the univariate and multivariate Cox regression analyses determined whether the risk score derived from the prognostic risk model could act as an independent prognostic indicator. In the univariate Cox regression analysis, the risk score (p < 0.001, HR = 3.996, 95% CI: 2.253–7.088, Figure 3H) was a potential hazard factor. The multivariate analysis also indicated that the risk score could serve as an independent prognostic factor (p < 0.001, HR = 3.282, 95% CI: 1.787–6.026, Figure 3I). Finally, the PCA plot indicated that the high- and low-risk groups could be well-separated with the three risk genes, but not all the genes or all the pyroptosis genes (Figures 3J–L).

The efficiency of the risk model was validated in a TCGA internal test cohort which included 184 patients and an ICGC external validation cohort. In the TCGA internal test cohort, we divided patients into high or low-risk groups according to the median risk score calculated by the formula in the training cohort. The risk score is an independent prognostic factor in the test cohort as indicated by univariate and multivariate Cox regression analysis (p = 0.001, HR = 3.225, 95% CI = 1.566–6.643 for univariate, Figure 4A; p = 0.01, HR = 2.843, 95% CI = 1.288–6.277 for multivariate, Figure 4B). The high-risk group tended to have more deaths and higher expression of risk genes (Figures 4C–E). The KM curve indicated that overall survival was lower in the high-risk group (p = 0.015, Figure 4F).

Moreover, the ROC curve suggested that the model exhibited an excellent predictive capability (AUC = 0.677, Figure 4G). PCA analysis indicated that the expression levels of the 3-risk gene could separate high-risk patients from low-risk in the TCGA test cohort (Figure 4H). We split the ICGC external validation cohort into the high- and low-risk groups based on the risk score (Figure 5A). The low-risk group had fewer deaths and lower expression of the risk genes (Figures 5B,C). The KM curve and ROC analysis suggested that the overall survival of the low-risk group was higher (p = 0.02, Figure 5D), and the model was reliable (AUC = 0.638, Figure 5E). Finally, the PCA plot indicated that the risk genes were well able to separate the two risk groups (Figure 5F).

We separated the patients in the TCGA-LIHC dataset into several subgroups according to the different clinical parameters. First, we investigated whether the high- and low-risk patients determined differences in survival. The clinical stratifications for the study included age (>65 vs ≤ 65), gender (female vs male), tumor grade (G3/4 vs G1/2), and AJCC stage (I/II vs III/IV). The KM curve showed that the high-risk patients had a poorer survival probability than the low-risk patients under the condition of age >65, female, male, G1–G2, G3–G4, or Stage I–II, except for age ≤65 or stage III-IV (Figures 6A–H).

We developed a novel prognostic nomogram to offer a reliable and quantifiable method for predicting the survival of the HCC patients based on the risk scores and clinical features, such as age, gender, grade, and stage (Figure 7). The nomogram could effectively predict the probability of the 1, 3, and 5-years overall survival in the HCC patients.

Given clinical information, the high- and low-risk groups were significantly correlated with the immune scores (Figure 8A). Specifically, the patients with lower immune scores favored high-risk scores compared to those with higher immune scores. In addition, the patients with higher stages picked higher risk scores (Figures 8B–E). Interestingly, the high-risk group had a high expression level of PD-L1 (Figure 8F). Furthermore, the GSEA analysis significantly enriched the immune-related signaling pathways, such as the B cell and T cell receptor signaling pathways (Figure 8G).

Further, the ssGSEA analyzed the differences in 16 kinds of immune cell infiltrations and 13 types of immune signal pathways. The high-risk group had decreased infiltrations of the B cells, mast cells, NK cells, and TIL and increased infiltrations of the DCS and macrophage compared to the low-risk group (p < 0.05, Figure 8H). Besides, the high-risk group exhibited suppression in the immune pathways, including the cytolytic activity and type II IFN response (p < 0.05, Figure 8I).

We analyzed the differential mRNA levels of the three independent prognostic genes in the normal liver tissues (from TCGA and GTEx) and HCC tissues (from TCGA) by GEPIA. The HCC tissues favored an increased mRNA level of GSDME, GPX4, and SCAF11 (Figures 9A–C). In addition, we examined the protein levels of these three genes using the HPA database. The results showed that the HCC tissues had higher protein levels of GSDME, GPX4, and SCAF11 (Figures 9D–I).

We further validated the expression levels of GSDME, GPX4, and SCAF11 in a mouse HCC model, which was constructed by knocking out the p53 and overexpressing the myc in the mouse liver (Figure 10A). The HCC model formed multiple tumors as shown in Figure 10B. As expected, the mouse HCC livers had remarkably increased GSDME, GPX4, and SCAF11 mRNA levels compared to the normal livers (Figures 10C–E).

In addition, we confirmed the mRNA levels of these independent prognostic genes in human culture cells. As shown in Figures 10F–H, the HCC cell line SK-Hep1 expressed significantly higher mRNA levels of GSDME, GPX4, and SCAF11 compared to the normal liver cell line HL-7702.

By analyzing the CellMiner database, the potential drugs were found to be correlated to these independent prognostic genes (Supplementary Table S4). Among the top 16 gene-drug correlations, 15 correlations pointed to the GSDME; the other correlation was SCAF11 (Figure 11). Surprisingly, the HCC drug lenvatinib positively correlated with the expression of GSDME (Cor = 0.453, p < 0.001, Figure 11K). To confirm this correlation, we treated HCC cell line Huh7 with Lenvatinib. As was expected, Lenvatinib treatment upregulated both the mRNA and protein levels of GSDME. Besides, Lenvatinib was able to induce the active form of GSDME (GSDME N-terminal).