A Novel Pyroptosis-related Genes Prognostic Signature Affects the Immune Inltration Patterns in Hepatocellular Carcinoma

Purpose Pyroptosis is an inammatory form of cell death associated with tumorigenesis and progression. However, the prognostic value of pyroptosis-related genes (PRGs) in hepatocellular carcinoma (HCC) have not been elucidated. Methods We downloaded mRNA expression proles and clinical information from TCGA and ICGC database. Then, differently expressed PRGs were screened to construct a multigene prognostic signature by least absolute contraction and selection operator (LASSO) Cox regression method in TCGA cohort. Date from ICGC was used to validate the robustness of this signature. Kaplan-Meier analysis was used to compare overall survival (OS) between high- and low-risk group. Univariate and multivariate Cox analysis were performed to identify the independent prognostic value of the signature. Gene set enrichment analysis (GSEA) was utilized to conduct GO and KEGG analysis. Single-sample gene set enrichment analysis was implemented to assess the immune cell inltration and immune-related function. TIDE algorithm evaluated the signicance of this signature in predicting immunotherapeutic sensitivity.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common serious life-threatening malignancies with a high cancer-related mortality rate throughout the world [1] [2] . Despite the continuous development of treatments for HCC, the improvement of survival rate is still limited due to the di culty in early diagnosis and high recurrence rate [3] . Therefore, it is important to explore better methods for early diagnosis and prognostic evaluation to geared towards ameliorating the treatment and prognosis of HCC patients.
Pyroptosis is a new form of programmed cell death, also known as cellular in ammatory necrosis, characterized with release of in ammatory cytokines and characteristic bubble-like protrusions [4] [5] . Gasdermin family proteins, including gasdermin-A to gasdermin-E and PJVK are the main executors of pyroptosis. Caspases play a central role in the assembly of in ammasomes and the initiation of pyroptosis. Through the canonical and non-canonical pathways, speci c caspase family proteins are activated, and cleave GSDMD, GSDMB or GSDME at the intramolecular binding between the N-terminal effector domain and C-terminal inhibitory domain [6] [7] . The cytotoxic N-terminal of gasdermin then oligomerizes and forms a 10 to 20 nm transmembrane pore, which promotes the maturation and release of pro-in ammatory cytokines (such as IL-1β and IL-18). Then it causes the in ltration of immune cells and in ammation in the microenvironment, which in turn leads to cell swelling and plasma-membrane rupture [8] . In addition to caspase, granzyme proteases from lymphocytes and natural killer cells can also mediate the cleavage of gasdermin and induce the pyroptosis of cells [9] .
Pyroptosis was rst studied in the immune defense against pathogen infection [10] . Recently, increasing evidence indicated pyroptosis a dual role of promoting and inhibiting in the pathogenesis of tumors [11] . On the one hand, it can effectively kill tumor cells by inducing apoptosis. On the other hand, with arousing in ammatory response, pyroptosis cause a suitable microenvironment for the malignant transformation and growth of tumor [5] . However, the role of pyroptosis in HCC and its rele in prognosis of HCC deserve further study.
In this study, we constructed a prognostic model based on pyroptosis related genes (PRGs). Then, we performed functional enrichment analysis to explore its underlying mechanism and the correlation with the tumor immune microenvironment. Collectively, our study aims to provide a new understanding of the role of pyroptosis in HCC, and explore the prognostic value of PRGs for HCC patients.

Data and information collection
We downloaded the total mRNA sequencing data and corresponding clinicopathologic information of 374 HCC and 50 normal tissues from the The Cancer Genome Atlas (TCGA) database on July 28, 2021 (https://portal.gdc.cancer.gov/). Another 231 patients' mRNA sequencing data and corresponding clinical information was obtained from the International Cancer Genome Consortium (ICGC) database for external validation cohort (https://dcc.icgc.org/projects/LIRI-JP). The involved data from these databases are open-access and our study abide by the data acquisition policies and publication guidelines of them.

Identifying the differentially expressed PRGs
We collected 90 genes related to pyroptosis from previously published reviews [4][12]- [24] , as shown in Supplementary Table 1. The mRNA sequencing data for the differentially expressed PRGs between tumor and normal tissues were identi ed using "limma" R package in R software 4.1.0 with a p value <0.05 for subsequent analysis. The "pheatmap" and "ggplot2" R packages were employed to drawn heatmaps, volcano plots and box plots.

Construction and validation of prognostic risk model based on PRGs
Cox regression analysis was used to assess the association between each PRG and survival status to assess the prognostic value of PRGs. Further, the least absolute shrinkage and selection operator (LASSO) -penalized Cox regression analysis were conducted using "GLMnet" R package to narrow the range of candidate genes and construct a prognostic model. The penalty parameter (λ) was determined and followed the minimum criterion (in other words, the λ value corresponds to the lowest partial likelihood deviation). Calculate the risk score using the following formula: risk score = where X represented the expression level of each retained PRGs and Y represented the corresponding regression coe cient. The mutation pro les of each gene involved in the risk model were analyzed in HCC via cBioPortal database (http://cbioportal.org) [25] . The median value of risk score was used to divide the HCC patients into high-risk and low-risk groups. "Rtsne" and "ggplot2" R packages were utilized to perform principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis for dimensionality reduction and verifying the distribution of different groups. Kaplan-Meier survival analysis was conducted to compare the overall survival (OS) between two subgroups. R packages "survival", "survminer", "survivalROC" and "timeROC" were employed to generate time-dependent ROC curve analysis. Univariate and multivariate Cox regression analysises were used for determination the independent prognostic value of the 8-gene model.

Gene set enrichment Analysis and Immune In ltration Analysis
To explore biological functions and pathways affected by the risk scores, gene sets database "c2.cp.kegg.v7.4.symbols.gmt" and "c5.go.v7.4.symbols.gmt" were used to conduct gene set enrichment analysis via the software Gene Set Enrichment Analysis (GSEA, version 4.0.1) [26] [27] . The FDR value q and normalized value p<0.05 were considered statistically signi cant.

Immune Function Analysis and Immunotherapy prodiction
Single-sample gene set enrichment analysis (ssGSEA) was performed with "GSVA", "ESTIMATE" and "GSEABase" R packages to score the immune functions, pathways, as well immune cells and stromal cells in ltration level, for each sample. Tumor immune dysfunction and rejection algorithm (TIDE, http://tide.dfci.harvard.edu) was used to compare the therapeutic effect evaluation of immune checkpoint blockade (ICB) in different groups.

Statistical analysis
The WilCoxon test were used to compare differences in gene expression between tumor and normal tissue and immune scores between high-and low risk groups. The two-sided log-rank test was utilized to perform Kaplan-Meier survival analysis for compare the OS of patients between two subgroups. Spearman correlation analysis was used to test the correlation between immune scores and risk scores of the prognosis risk model. All statistical analyses were performed with R software (V4.1.0) Results 1.Identi cation of differentially expressed PRGs.
First, the transcriptional pro les of 90 PRGs were investigated in the TCGA mRNA-seq data for 50 normal tissues and 374 HCC tissues. A total of 65 differentially expressed PRGs were screened out (All p<0.05). Among them, 58 genes were enriched in HCC tissues while the 7 other genes were downregulated compared to normal tissues. The mRNA expression levels of these genes are presented in the form of volcano map and heatmap ( Figure 1A, B). The correlation network of these differentially expressed PRGs was shown in the Figure 1C.

2.Classi cation of HCC based on differentially expressed PRGs
A consistent cluster analysis using the data of 374 HCC cases from TCGA database was preformed to explore the association between the transcriptional pro les of 65 PRGs and HCC subtypes. We increased the clustering variable (k) from 2 to 9 and found that the intra-group correlations were highest and the inter-group correlations were low when k=2, suggesting that 374 HCC patients could be divided into two clusters according to the expression of 65 PRGs (Figure 2A). Then we preformed PCA analysis to verify the discrete distribution of the two clusters (Supplementary Figure 1). Kaplan-Meier curves showed that OS of cluster 1 was signi cantly better than that of cluster 2 ( Figure 2B p=0.004872). Heatmap was drawn to show the relationship between the expression pro les of 65 PRGs and clinical features, such as tumor size (T1-4), lymph node metastasis(N0-1), distant organ metastasis(M0-1), clinical stage (stage I-IV), histologic grading (G1-G4), age (≤60 or >60 years), gender (female or male) and survival status (alive or dead) ( Figure 2C). It could be observed that there were signi cant differences in tumor size and clinical stage between the two clusters ( Figure 2D, Figure 2E, both p<0.001).

3.Construction of prognostic risk model according to PRGs
HCC samples with incomplete clinical information were removed and survival-related PRGs were screened by univariate COX regression analysis. The results revealed that a total of 11 PRGs (CASP8, CHMP3, HSP90AA1, BAK1, HSP90AB1, GSDMC, GSDME, STK4, CGAS, NOD2) were identi ed correlated with OS (p<0.01) and retained for further analysis ( Figure 3A). LASSO-Cox regression analysis was preformed to construct an 8 genes prognostic model (HSP90AA1, CHMP4B, BAK1, HSP90AB1, GSDMC, GSDME, STK4, NOD2) based on the optimal value of penalty parameter (λ) ( Figure  ). Based on the median value of risk score calculated by the formula, HCC patients in TCGA cohort were divided into low-risk and high-risk subgroups ( Figure 3D). The t-SNE analysis and PCA showed that patients with different risks were well distributed into two directions ( Figure 3E). In addition, the scatter plot shows that high-risk patients had more deaths and shorter survival time than low-risk patients ( Figure 3F). Consistently, the Kaplan-Meier curve showed that high-risk patients had signi cantly poorer OS than low-risk patients ( Figure 3G). Time-dependent ROC curves were generated for evaluation of the sensitivity and speci city of the prognostic model and the area under the curve (AUC) reached 0.717 at 1 year, 0.665 at 2 years, and 0.658 at 3 years ( Figure 3H).

4.External Validation of the 8-Gene risk Signature
As the external validation set for the prognostic model, patients in the ICGC cohort were also divided into high-risk or low-risk groups based on the median risk score calculated by the same equation ( Figure 4A). PCA and t-SNE analyses con rmed a satisfactory separation of patients between the two subgroups ( Figure 4B, C). As shown in Figure 4D to 4F, the OS and ROC analyses of these two subgroups showed similar results to the TCGA cohort.

Prognostic Value of the Prognostic Risk Model
Univariate and multivariate Cox regression analyses were preformed to assess independent predictive value of the risk scores derived from 8-genes model for OS of HCC patients. Univariate Cox regression analysis revealed that risk score was negatively correlated with prognosis, with a hazard ratio 3.303 Immunohistochemistry from the Human Protein Atlas (HPA) database for the prognostic PRGs showed that the expression in cancer tissues were mostly higher than that of normal tissues ( Figure 6A). Furthermore, we analyzed the expression and mutation pro les of the 8 PRGs based on the HPA database and cBioProtal database. Data form the 353 samples with mutation and Copy Number Alterations (CNA) showed that 67 patients (19%) had heterogeneous mutations ( Figure 6B). Thereinto, 0.57% patients had deep deletion, 0.85% patients had missense mutation and 0.28% patients had truncating mutation in HSP90AA1; 0.85% patients had ampli cation while 0.57% patients had structural variant in CHMP4B; 1.42% patients had ampli cation BAK1; 3.12% patients had ampli cation, 0.28% patients had truncating mutation and 0.57% patients had missense mutation in HSP90AB1; 10.2% patients had ampli cation, 1.13% patients had missense mutation and 0.28% patients had deep deletion in GSDMC; 0.28% patients had ampli cation and truncating mutation respectively, while 0.57% patients had missense mutation in GSDME; 0.28% patients had missense mutation STK4; 0.28% patients had ampli cation and deep deletion respectively, while 1.42% patients had missense mutation in NOD2.

7.Analysises of Biological Function and Pathway.
In order to further explore the biological processes and pathways through which the 8-gene signature affect the prognosis in HCC, GSEA were performed to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses between the high-and low-risk groups. 20 GO functions and 20 KEGG pathways were respectively enriched in the high-risk group with the FDR value q and normalized value p<0.05. The results turned out that some functions and pathway involved in regulation of cell cycle phase transition were enriched in the high-risk group. And in addition to this, several tumor-associated pathways were signi cantly enriched in high-risk group, including apoptosis, cell cycle, WNT, TGF-β, VEGF and MAPK signaling pathway. Interestingly, the high-risk group was signi cantly correlated with some pathways related with immune response, including FC-γ receptor mediated phagocytosis, nod like receptor signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway, toll like receptor signaling pathway, natural killer cell mediated cytotoxicity, mTor signaling pathway ( Figure 7A and 7B).

Comparison of Immune status between subgroups
Since the relevance between the risk score and immune response, we further compared the immune in ltration score and immune-related function between the high-and low-risk groups in both the TCGA and ICGC cohorts by using ssGSEA. The results showed that the high-risk group was associated with higher levels of immune cell in ltration, including aDCs, DCs, iDCs, pDCs, Macrophages, CD8+ T cells, Th1 cells, Th2 cells, TIL, Treg and Tfh in the TCGA cohort ( Figure 8A) as well Treg, aDCs, DCs, iDCs, pDCs, Figure 8B). In addition, data from the TCGA cohort suggested that the high-risk group had higher scores in terms of immune checkpoint, in ammation-promoting, APC co-inhibition, APC co-stimulation, T cell co-stimulation, T cell co-inhibition, HLA and MHC class I. Similar results were obtained for data from ICGC cohort. The correlation analysis suggested that risk score was also signi cantly positively correlated with immune score and stromal score in ICGC cohort (p<0.001), while correlated with immune score in TCGA cohort (p= 0.0011).

Macrophages, Neutrophils, CD8+ T cells, T helper cells, TIL, Treg, Th1 cells and Th2 cells in the ICGC cohort (
We further analyzed the correlation between the risk score and the expression level of immune checkpoints in TCGA cohort. The results suggested that risk score was signi cantly positively correlated with the expression levels of B7-H3, B7H4, TIM-3, LAG3, CTLA-4, BTLA, PD-L1, PD-L2 and PD-1 (Figure 9), which combined with TIDE scores implied that patients in high-risk group might be more sensitive to immune checkpoint blockade therapies (Supplementary Figure 4, p<0.0001).

Discussion
This study is the rst to construct a novel signature based on PRGs in HCC, which can be used to predict prognosis and re ect the immune status for HCC patients. The risk score of PRGs prognostic model combined with traditional TNM staging was proven to be better for predicting the prognosis of HCC in both the discovery TCGA cohort and the validation ICGC cohort.
Pyrolysis that occurs in tumor cells and tumor microenvironment (TME) causes different clinical outcomes through the interactions between pyroptotic cells and TME [28] . On the one hand, as one in ammatory form of programmed cell death, in addition to directly eliminating tumor cells, pyroptosis in tumor cells can also stimulate in ammatory response and reactivate anti-tumor immunity. Wang et al. found that active GSDMA3 could result in immune dependent tumor clearance. HMGB1 released by the pyroptosis of melanoma could activate dendritic cells to enhance the in ltration of CD8+ and CD4+ T cells to inhibit tumor progression and prolong survival [29] . The induction of pyroptosis could improve the e cacy of immune checkpoint therapy [30] . On the other hand, some of the key components in pyroptosis, such as NLRP3 in ammasome and Gasdermin family proteins, can promote the proliferation, migration and the epithelial-stromal transformation of tumor cells and correlate with poor prognosis [31] [32] . Various studies proved that IL-18/IL-1β released by pyroptotic tumor cells can induce the proliferation, immunosuppression, angiogenesis and metastasis of varieties of tumors [33] [34] . The proin ammatory factor HMGB1 produced by GSDME-mediated pyroptosis process promotes the proliferation of colon cancer via activating the ERK1/2 pathway [35] . Pyroptosis-released factors from tumor cells also can again trigger pyroptosis of macrophages, leading to cytokine release and subsequent cytokine release syndrome [36] . Further, adjacent normal cells stimulated by in ammatory factors released by pyroptosis may transform into tumor cells [37] . In recent years, some researches focus on the notion that chronic in ammation and tumor necrosis related to pyroptosis favor tumor progression. However, the PRGs signature as prognostic model for HCC has not been reported.
In this study, we rst compared the expression of 90 PRGs in HCC tissues and adjacent peritumoral tissue in the TCGA cohort. Then, 65 differentially expressed PRGs were screened out and analyzed for their relationship with overall survival. We also found that HCC patients could be clustered into two groups based on the 65 differentially expressed PRGs. In addition, there were signi cant differences in the clinicopathological characteristics between different pyroptosis clusters, such as tumor stage, tumor size, and survival, which suggesting that the heterogeneity of pyroptosis in tumors may results in different clinical outcomes. Next, we constructed a prognostic model that integrates 8 PRGs through LASSO regression analysis and found that the OS of patients in different risk subgroups is signi cantly different, both in the discovery set and external validation set. Furthermore, Multivariate Cox regression analysis showed that risk score was an independent predictor for OS.
The eight PRGs that constitute the prognostic model in our study include HSP90AA1, CHMP4B, BAK1, HSP90AB1, GSDMC, GSDME, STK4 and NOD2, all of which were expressed signi cantly higher in HCC tissues compared to normal tissues and associated with poor prognosis. Both HSP90AA1 and HSP90AB1 are members of HSP90 family, which has been considered promoting the formation and progression of cancer via varieties of carcinogenic signal pathways [38] . In HCC, HSP90 promotes glycolysis and proliferation, and inhibits apoptosis by mediating PKM2 phosphorylation [39] . Under hypoxia, the transcription of GSDMC was enhanced to switch TNFα-induced apoptosis to pyroptosis in cancer cells and correlated with poor survival [40] . It's been demonstrated that miltirone inhibited HCC cells growth through GSDME-dependent pyroptosis [41] , which seemed to contradict the results from TCGA data base. It may be explained by the reason that different levels of GSDME play different roles in HCC development, but the underlying mechanism remains to be addressed. BAK is a key regulator of mitochondria-mediated apoptosis and identi ed as an oncogenic regulator of HCC. Previous studies suggested BAK plays a promoting role in the stage of hepatocarcinogenesis. BAK de ciency inhibited hepatocyte apoptosis and reduce the tumorigenic rate of liver cancer [42] . CHMP4B, a subunit of the endosomal sorting complex required for transport-III, plays an important role in the budding and scission of membrane vesicles of the cell and the late stage of mitotic cell division [43] . NOD2, one of the nummber of intracellular PRRs, exaggerates the in ammatory response through RIP2 signaling and promotes hepatocarcinogenesis in mice. The overexpression of NOD2 in HCC also impairs DNA repair through lamin A/C degradation and correlated with poor prognosis of HCC patients [44] . STK4, as known as MST1, is a serine/threonine kinase in the Hippo signaling pathway, which role in cancer is yet well illuminated.
Recent studies show that functions of Hippo signaling changes with temporal and spatial contexts [45] . In our study, both the TCGA and ICGA database suggested that STK4 expression is higher in tumor tissues than that in normal tissues, and positively correlated with poor prognosis. It may be explained by the spatiotemporal heterogeneity of STK4, which has a dynamic expression pattern and play different roles at different stages of cancer.
We then preformed GSEA analysis to explore the potential biological processes and pathways associated with these PRGs. Some representing oncogenic pathways were signi cantly enriched in high-risk group, including MAPK, VEGF, WNT, TGF-β, NOTCH signalings. In addition, as expected, some immune responserelated biological functions and pathways were signi cantly enriched in high-risk group, suggesting pyroptosis process is accompanied with the secretion of in ammatory factors and the activation of in ammasomes.
The difference between high-and low-risk groups in immune in ltration and immune-related functions further validated that in ammation response in the tumor microenvironment might be affected by the prognostic risk model to make it predictive. As observed from the results, the in ltration levels of DCs and macrophages in tumor tissues in the high-risk group is higher than that in low-risk group. Previous studies have demonstrated that tumor-in ltrating DCs play an immunosuppressive role, and plasmacytoid dendritic cells pDCs produced low amounts of INF-α, which facilitates the expansion of regulatory T cells and immune tolerance [46] . There also evidence demonstrated that the intratumoral in ltration of macrophages promoted HCC growth, metastasis, and resistance to sorafenib [47] . The difference of Treg levels between low-and high-risk groups also indicates immune function suppression in the high-risk group, which may response to regulating the overactive in ammatory reactions caused by pyroptosis. Of note, the high-risk group has lower activity of type II IFN response with function of tumor immune surveillance and antitumor mechanisms than that of low-risk group [48] , While Type I IFNs principally produced by innate immune cells have no difference between the two groups [49] .
Furthermore, in our study, the expression of immune checkpoints, including PD1, PDL1, PDL2, CTLA4, B7-H3, B7-H4, TIM-3, LAG3 and BTLA in high-risk group were signi cantly higher than those of low-risk group, and positively correlated with the risk score [50] . This may also be one of the reasons for the impairment of immune functions in the high-risk group and can be used to predict the sensitivity to immunotherapy for HCC patients. Based on these ndings, it is reasonable to assume that the poor survival outcome of high-risk HCC group may be caused by attenuated anti-tumor immunity.

Conclusion
In conclusion, we evaluated differential expression of PRGs between normal and HCC tissues and de ned a novel 8-gene prognostic risk model based on the pyroptosis for HCC. Then, we veri ed their relationship with prognosis and clinical features both in TCGA and ICGC cohorts. Moreover, the model risk score was found to be associated with immune response and in ltration of immune cells in the tumor microenvironment. The speci c underlying mechanism of PRGs affecting prognosis and the interaction between pyroptosis and tumor microenvironment need further exploration. Our research established a new gene model for predicting the prognosis of HCC patients and revealed potential effect of PRGs on tumorigenesis and development of HCC, which also provided a valuable insight into the relationship between pyroptosis and immunotherapy in HCC.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
All the authors consented to the publication of this research.

Figure 1
Expression pro le of the candidate 90 PRGs and the interactions among them A. Volcano plot for identi cation of differentially expressed 65 pyroptosis-related genes from the total of 90 PRGs. Red represents high expression, green represents low expression, black represents no difference between HCC and normal tissues. B. Heatmap (green, low expression level; red, high expression level) of the differentially expressed pyroptosis-related genes between the normal (N, blue) and the tumour tissues (T, red). C. The correlation network of the differently expressed 65 genes.     and ICGC cohort (B). C, D Comparison of the ssGSEA scores of 13 immune-related pathways between low-(green) and high-risk (red) group in the TCGA cohort (C) and ICGC cohort (D). E, F. The correlation analysis between the risk score and immune score, as well as the risk score and stromal score in TCGA cohort (E) and ICGC cohort (F). The correlation analysis between risk score and the expression levels of immune checkpoints, such as B7-H3, B7H4, TIM-3, LAG3, CTLA-4, BTLA, PD-L1, PD-L2 and PD-1