Pyroptosis-related genes (PRGs) were obtained from msigdb (https://www.gsea-msigdb.org/) ( Supplementary Table 1 ). The clinical and mRNA expression data in the training cohort are downloaded from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) (https://xenabrowser.net/datapages/; https://dcc.icgc.org/), respectively. The validation cohorts including GSE10142 and GSE14520 datasets were downloaded from https://www.ncbi.nlm.nih.gov/geo/. Batch effects were removed with the R package (sva), and batch-removed results were shown in Supplementary Figures 1A, B . The CPTAC cohort was download from https://www.biosino.org/node/project/detail/OEP000321, and transcriptome and proteome were used in this study and the clinical information of patients is presented in Supplementary Table 2 . The treatment cohort GSE109211 for sorafenib was also downloaded from the GEO database.

Unsupervised clustering analysis based on 40 PRGs was done in the R package ConsensusClusterPlus (17). Differentially expressed genes (DEGs) analysis was identified with fold change>2 and p<0.05 using R package edgeR (18).

DEGs significantly associated with survival were analyzed using Cox regression of the R packages survminer and survival (19, 20). The R package glmnet was used to construct the risk model (21).

In the LASSO regression, a 5-cross-validation model with binomial deviance minimization criteria was implemented in the training set. The lambda with min was used for feature selection. Based on the selected gene markers, the predicted probability of each patient could be calculated using the following formula:

where n is the number of genes, Exp_i is the value of each gene, and Coef_i is the estimated regression coefficient of the selected genes.

The median value was used to divide patients into the high-risk and low-risk groups. The 1-, 2- and 3-year ROC curves and the ROC curves compared to predict the survival status with the R package survivalROC (22).

The prognostic performance of the PRG-related risk model was further validated in the external validation cohort (GSE10142 and GSE14520 datasets) and the CPTAC dataset. Patients are stratified into low- and high-risk groups using the median value as the cutoff. The independent prognostic value of the PRGs-related risk score was assessed using univariable and multivariable analysis, respectively. The clinical information of GSE10142 other than overall survival (OS) could not be downloaded from the public database, so the independence test of the risk score was not performed in GSE10142. The prediction efficiency of the nomogram was evaluated by calibration curves.

Continuous variables were compared by Mann-Whitney U test or Kruskal-Wallis H-test and categorical variables were compared by Chi-Squared test or Fisher exact test. Principal Components Analysis (PCA) was performed by R package FactoMineR. Survival was estimated by Kaplan-Meier curves, with the p determined by a log-rank test. CIBERSORT (https://cibersort.stanford.edu/) was used to estimate the infiltration of immune cells based on mRNA data. KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis of genes with DEGs and Gene Set Enrichment Analysis (GSEA) were performed using R package clusterProfiler. The non-supervised clustering of mRNA data was performed by the k-means method (R package, ConsensusClusterPlus). The area under curve (AUC) and 95% confidence interval (CI) were generated to evaluate the model performance. Immunohistochemical data of the proteins were downloaded from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). Drug susceptibility prediction was performed by R package oncoPredict. The correlation between risk score and immune cells was evaluated using person correlation. P<0.05 was considered to be statistically significant and all p were two-sided. All the statistical tests were performed using R software, version 4.0.1 (R Foundation for Statistical Computing Vienna, Austria).

The 40 core genes that regulate or mediate pyroptosis signaling have been identified previously (23) and were available in the Molecular Signatures Database (MSigDB) ( Supplemental Table 1 ). Gene mutations of pyroptosis-related genes (PRGs) occurred in the majority of hepatocellular carcinoma samples (54.4%) ( Figure 1A ), implying a potential role of PRGs in tumorigenesis and development of HCC. The TP53 gene was the most commonly mutated (30%), whereas mutations in other genes (DHX9, GSDMC, etc.) were less frequent. Frameshift deletions were the most common type of mutation, followed by frameshift insertions and missense mutations. Moreover, most PRGs had CNV gain or loss ( Figure 1B ). GSDMC was the most amplified gene with a frequency of 14.8% in HCC patients, in contrast, HMGB1 had the highest CNV loss with a frequency of 11.6%. In addition, the mRNA expression of 77.5% (31/40) PRGs was significantly different between tumor and normal tissues ( Figure 1C ). Most of PRGs displayed higher expression levels in tumor samples compared with normal controls except for AIM2, IL1B, and NLRC4, of which the mRNA expression was decreased in tumor samples. Meanwhile, PCA based on the mRNA expression of PRGs could well differentiate HCC samples from normal samples ( Figure 1D ), suggesting the abnormal expression of PRGs in hepatocellular carcinoma.

To investigate the role of pyroptosis signaling in HCC, we further explored the role of PRGs in the prognosis of patients with HCC depending on different transcriptomic patterns. After eliminating the batch effect, two datasets derived from TCGA and ICGC were grouped as the training cohort. The patient characteristics were depicted in Supplementary Table 2 . Unsupervised hierarchical clustering of transcriptome profiles of 40 PRGs identified three clusters ( Figure 2A ) in the training cohort. As delineated in the heatmap ( Figure 2B ), ranging from the cluster 1 to 3, the mRNA expression of PRGs was gradually increased. The principal components analysis also depicted the difference between the three clusters ( Figure 2C ). Gene Set Enrichment Analysis (GSEA) analysis on the signatures of cancer hallmarks confirmed that adipogenesis and oxidative phosphorylation were enriched in the cluster 1 but TNF-α signaling via NF-KB and WNT beta-catenin signaling were enriched in the cluster 2; conversely, the cluster 3 was characterized by genes in interferon-alpha response signaling ( Figure 2D ). Survival analysis found significant survival differences among the three clusters (log-rank p for trend = 0.004 Figure 2E ). The cluster 2 exhibited shorter OS compared to the cluster 1 (p=0.003; HR, 1.74; 95% CI 1.26-2.39, Figure 2E ). A two-by-two comparison of differentially expressed gene (DEG) analysis was performed and a total of 318 overlapped genes were found among the three pairs ( Figure 2F ). KEGG functional enrichment analysis was performed on the 318 overlapped DEGs, and significant pathways were partially presented in Figure 2G , such as the HIF-1 signaling pathway, natural killer cell-mediated cytotoxicity, and PPAR signaling pathway. Altogether, these results suggested that PRGs might play an important role in the prognosis of HCC.

To specifically identify the prognosis-related core genes, a prognostic risk model for HCC was constructed based on DEGs using the LASSO regression method. We first selected 211 candidate genes which both existed in the training (TCGA and ICGC) and the validation cohorts (GSE14520 and GSE10143, Figure 3A ). Among these genes, totally 68 genes were significant associated with OS by Cox proportional regression, and then a 24-gene risk model was constructed by LASSO regression ( Figures 3B, C ). Risk score = CYP2A6*(0.054) + PPT1*(0.4319) + G6PD*(0.6332) + CD97*(1.3631) + LGALS9*(-0.5509) + S100A9*(0.1446) + SEC14L2*(-0.2186) + NFE2L3*(0.0596) + FABP3*(-0.1861) + CD2*(-0.7902) + BATF*(0.7061) + PON1*(-0.336) + FMO3*(-0.3352) + CSF1*(0.2407) + CYP4A11*(0.1653) + ATP2A3*(-5.13) + MUC1*(-0.0063) + LAIR1*(-1.5441) + ATP1B3*(1.2079) + AZGP1*(0.0064) + NCF2*(1.0057) + CYP2C9*(-0.1876) + TMSB10*(-0.0666) + TACC3*(1.2139). Patients were divided into the high- and low-risk groups according to the median value of risk score. The transcriptomic profile of the selected 24 genes in the training set is illustrated in Figure 3D . Patients in the high-risk group had shorter OS compared with those in the low-risk group (HR=3.06, 95% CI 2.22-4.24, p<0.001, Figure 3E ), and AUC value of 2-year survival was 0.77 in the training cohort ( Figure 3F ).

We further tested whether the risk score was an independent prognostic factor of HCC. In the univariable analysis, besides risk score, several other variables such as TNM stage, HBV or HCV infection, and TMB were also associated with the OS with HRs (95% CI) of 2.57 (1.90-3.47), 2.49 (1.82-3.40), 1.82 (1.36-2.44) and 1.91 (1.39-2.62), respectively ( Supplementary Tables 3, 4 ). In the multivariable analysis, the association between the risk score and OS remained significant (HR=2.82, 95% CI 1.89-4.20, p<0.001, Supplementary Tables 3, 4 ). These results indicated that the risk score was associated with a better prognosis independent of TNM stage, HBV or HCV infection, and TMB.

To validate the performance of the 24-gene risk model, public datasets in GEO (GSE14520 and GSE10143) were collected together as the validation cohort after removing batch effect. The patients were also stratified into the high-risk and low-risk groups using the median value in the validation set. The transcriptomic profile of the selected 24 genes in the validation set is illustrated in Figure 3G . Consistently, patients with the high-risk score had shorter OS (HR=1.61, 95% CI 1.13-2.28, p=0.008, Figure 3H ) in the validation set, and the AUC value of 2-year survival was 0.65 ( Figure 3I ). In the multivariable analysis, the risk score was also associated with OS (HR=1.70, 95% CI 1.08-2.69, p<0.001, Supplementary Table 5 ), which indicated to be an independent factor of prognosis in HCC.

Next, a nomogram was constructed to predict the 1-, 2-, and 3-year OS rates in HCC patients integrating the risk score, tumor stage, HBV or HCV and TMB ( Figure 4A ). Calibration plots indicated that the nomogram had a good predictive power for 1-year and 3-year survival rates ( Figures 4B, C ). Compared to 1-year AUC value of 0.79 of pyroptosis risk score, the AUC value of the nomogram was up to 0.82 ( Figure 4D ). The nomogram further improved the performance and facilitated the clinical practice of the risk model.

We further validated the prognostic value of the risk score in the CPTAC dataset at mRNA and protein levels, which consisted of the transcriptomic and proteomic data of 159 patients with HCC and corresponding clinical information. Consistent results were observed at the mRNA level that patients with the high-risk score had worse OS compared to those with the low-risk score (HR=2.99, 95% CI 1.67-5.35, p<0.001, Figure 5A ) and AUC of 2-year survival was 0.69 ( Figure 5B ). Additionally, the robust performance was also observed in the multivariable analysis (HR=1.59, 95% CI 1.02-2.49, p=0.042, Supplementary Table 6 ).

To explore whether the risk models were also applicable at the protein level, the correlations between mRNA expression and protein expression of these candidate markers were examined. Finally, protein levels of 16 genes were significantly correlated with mRNA expression (p<0.05, Figure 5C ). Using the same regression coefficients as the mRNA model, these proteins were also calculated for risk scores in the CPTAC cohort. Patients with the high-risk score tended to have a relatively shorter OS (HR=1.57, 95% CI 0.92-2.68, p=0.094) ( Figure 5D ). We further selected the top five most related proteins, characterized as Pearson correlation coefficient (R) above 0.8, including CYP2A6, CYP2C9, G6PD, FMO3, and SEC14L2. Unsurprisingly, the high-risk patients in the 5-proteins model also had worse OS compared with the low-risk group (HR=2.97, 95% CI 1.66-5.31, p<0.001, Figure 5E ) and the model had a better predictive ability with a 2-year AUC of 0.68 ( Figure 5F ).

To further confirm the importance of CYP2A6, CYP2C9, G6PD, FMO3, and SEC14L2, the Human Protein Atlas (HPA) database was used to compare their protein expression in normal and HCC tissues. As demonstrated in Figure 6A , the expression of CYP2A6 and CYP2C9 were relatively lower in tumor tissues, while the expression of G6PD in HCC tumors was higher in tumor tissues than in normal tissues. In the proteomic data of the CPTAC cohort, G6PD was also significantly highly expressed in HCC tumor tissues compared to normal tissues, while the remaining four proteins were significantly less expressed, consistent with the results of immunohistochemistry ( Figure 6B ). KM survival curves showed that HCC patients with higher CYP2A6, CYP2C9, FMO3, and SEC14L2 protein expressions had longer OS and those of higher G6PD had shorter OS in the CPTAC cohort ( Figures 6C–G ).

To investigate the differences in physiological function between the high- and low-risk groups of patients, we compared the molecular mechanisms in immunity and angiogenesis. As shown in Figure 7A , blood vessel morphogenesis, sprouting angiogenesis and nine immunity-related pathways were significantly enriched in the patients with the low-risk score via GSEA. We next investigated the correlation between the expression of the 24 genes in the risk model and immune cell infiltration abundances and angiogenesis-related genes. The infiltration abundance of CD8-positive T cells was found to be significantly positively correlated with the expression of ATP2A3, CD2, LGALS9, and NFE2L3 and negatively correlated with the expression of FABP3 (p<0.05; Figure 7B ). Activated natural killing cell infiltration was also correlated significantly positively with CD2 and TMSB10 expression and negatively with SEC14L2 and G6PD (p<0.05; Figure 7B ). Twenty-four angiogenesis-related genes were also applied in correlation analysis with genes in the risk model ( Figure 7C ). In terms of gene expression, HIF1A and MMP9 were significantly correlated with all genes in the risk model (p<0.05; Figure 7C ). In terms of VEGFA, one of the core targets of anti-angiogenic drugs was also found to be significantly positively correlated with ATP2A3, G6PD, TACC3, and MUC1 and negatively correlated with PPT1, S100A9, TMSB10, FABP3, and BATF (p<0.05; Figure 7C ).

At last, we applied the risk model to the GSE109211 cohort containing 67 cases with advanced HCC treated with sorafenib. In the low-risk group, the proportion who responded to sorafenib was 45.5%, higher than that (17.6%) in the high-risk group (p=0.02, Figure 7D ). To compensate for the absence of other anti-angiogenic drugs in the treatment cohort, drug sensitivities were predicted by computation algorithm. The predicted IC₅₀ of forentinib, Lenvatinib, linifanib, and sunitinib were all significantly lower in the low-risk group ( Figure 7E ).