A Pyroptosis-Related Gene Signature to Predict Patients' Prognosis and Immune Landscape in Liver Hepatocellular Carcinoma

Background Liver hepatocellular carcinoma (LIHC) is a malignance with high incidence and recurrence. Pyroptosis is a programed cell death pattern which both activates the effective immune response and causes cell damage. However, their potential applications of pyroptosis-related genes (PRGs) in the prognostic evaluation and immunotherapy of LIHC are still rarely discussed. Methods Comprehensive bioinformatics analyses based on TCGA-LIHC dataset were performed in the current study. Results A total of 33 PRGs were selected to perform the current study. Of these 33 PRGs, 26 PRGs with upregulation or downregulation in LIHC tissues were identified. We also summarized the related genetic mutation variation landscape. GO and KEGG pathway analysis demonstrated that these 26 PRGs were primarily associated with pyroptosis, positive regulation of interleukin-1 beta production, and NOD-like receptor signaling pathway. An unfavorable OS appeared in LIHC patients with high CASP3, CASP4, CASP6, CASP8, GPX4, GSDMA, GSDME, NLRP3, NLRP7, NOD1, NOD2, PLCG1, and SCAF11 expression and low NLRP6 expression. A prognostic signature constructed by the above 14 prognostic PRGs had moderate to high accuracy to predict LIHC patients' prognosis. And risk score was correlated with the expression of CASP6, CASP8, GPX4, GSDMA, GSDME, NLRP6, and NOD2. Of these 7 genes, CASP8 was identified as the core gene in PPI network. Moreover, lncRNA MIR17HG/hsa-miRNA-130b-3p/CASP8 regulatory axis in LIHC was also detected. Conclusions The current study indicated the crucial role of PRGs in the prognostic evaluation of LIHC patients and their correlations with tumor microenvironment in LIHC.


Introduction
Globally, liver cancer accounted for over 800000 new cases and caused over 700000 cancer-related death in 2018. Of these new cases, 85% were diagnosed as liver hepatocellular carcinoma (LIHC), and over 40% were diagnosed at advanced stage [1]. LIHC is the most common subtype of liver cancer. Although some risk factors have been identified, including hepatitis B infection, liver cirrhosis, and alcoholic and nonalcoholic fatty liver diseases, and many approaches have been utilized in clinical practice, mainly including surgical resection, liver transplantation, chemotherapy, targeted drug treatment, and immunotherapy, the 5-year survival rate of LIHC in some developing countries is still only 18% [2][3][4]. And the minority of LIHC patients are eligible for these approaches due to high costs of treatment, serious drug adverse reactions, and multidrug resistance of tumor [5]. In a word, the current developments in the diagnosis, treatment, and prognostic evaluation of LIHC cannot meet the demand of patients. In the past 30 years, aberrant expression of genes in tumor tissues and their potential applications in clinical practice have been focused by clinicians. Thus, exploring potential gene signatures for therapeutic and prognostic assessment of LIHC is significant clinically.
With the rapid development of molecular biology and the enrichment of genomic data, comprehensive analysis on PRG signatures becomes feasible. In the present study, bioinformatic data was utilized to explore the expression profiles, prognostic performance, and related regulation axis of PRG signatures in LIHC. Moreover, we also explored the association between immune cell infiltration and PRG signatures in LIHC. The above findings may offer novel insights on the prognostic evaluation and therapy of LIHC.

Materials and Methods
2.1. Datasets. The genomic data of LIHC patients and related clinical data of these patients, including gender, age, tumor grade, and survival outcome, were downloaded from the Cancer Genome Atlas (TCGA) database on August 1, 2021. The TCGA-LIHC dataset (N = 374) was selected to perform the analyses. The data of copy number variation (CNV) and somatic mutation in LIHC were also extracted from TCGA database and UCSC Xena, respectively. Using R software V4.0.3, statistical analyses were performed. In addition, the P value cutoff was set as 0.05.

Mutation Analysis and Functional Enrichment
Analysis of PRGs. Mutation categories and mutation frequency of the 33 PRGs as well as their waterfall plots were constructed by "maftools" package in R. And the location of variation of these genes on 23 chromosomes was shown by "RCircos" package. Functional enrichment analysis, including Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, was performed to explore the potential molecular mechanisms and biological functions of target genes by Metascape (http://metascape .org/gp/index.html#/main/step1) [17]. Moreover, GO analysis included biological process (BP), cellular component (CC), and molecular function (MF) analysis.

Consensus Cluster
Analysis. Setting the threshold as follows: iteration = 100 and resample rate = 80%, consensus cluster analysis was conducted by "ConsensusClusterPlus" package in R. "Survival" package was utilized to evaluate the difference of each cluster in LIHC. Next, the difference of clinical paraments in each cluster was explored by "pheatmap" package. Furthermore, the ImmuneScore, ESTIMA-TEScore, StromalScore, and abundance of immune cell in each cluster were calculated by the ESTIMATE algorithm [18]. The comparison of above indexes was visualized by the "vioplot" and "ggpubr" package.

Prognostic Analysis and Prognostic
Model Construction of PRGs. The prognostic value of differently expressed PRGs was assessed by the "survival" package in R and was expressed as prognostic forest maps. The differently expressed PRGs with OS difference in low-expression and high-expression group were defined as prognostic PRGs. According to the above prognostic PRGs, LASSO Cox regression analysis in "glmnet" package was used to construct prognostic model. Those PRGs that constituted prognostic model were selected for further analyses. According to the median risk score, patients were divided into low-and high-risk group, and the OS curves of the two groups were compared. The predictive accuracy of this prognostic model was assessed by time ROC analysis. The univariate and multivariate Cox analysis was performed to excavate the influence of clinical factors, including age, gender, tumor grade, clinical stage and TNM stage, and risk score on the prognosis of LIHC patients. Subgroup analyses on the influence of clinical factors on the prognosis of low-risk and high-risk group were also performed. We were also interested in the difference of risk score in different subgroups, including age, gender, cluster, tumor grade, clinical stage, Immune-Score, and TNM stage. Besides, the correlation between abundance of immune cells and risk score was also analyzed.
2.6. PPI Network Construction of PRGs in LIHC. The PRGs that made up risk score were selected via above analyses. Using STRING (https://string-db.org/), a PPI network was constructed to illustrate the association among these genes and seek for the core gene [19]. Further analyses were performed to indicate the potential mechanisms, biological functions, and applications of the core gene in LIHC.  [20] is a bioinformatics platform that aims to illustrate the correlation between target gene and various immune cells. The "Gene" module was used to validate the correlation between the expression of CASP8 and the abundance of B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell. Next, we investigated the correlation between somatic copy number alterations (SCNAs) of CASP8 and immune infiltration level of different immune cells by the "SCNA" module. Furthermore, tumor mutation burden (TMB) and microsatellite instability (MSI) analysis were used to assess 9 Type P R K A C A C A S P 9 C A S P 4 C A S P 6 C A S P 8 C A S P 5 C A S P 1     2.8. mRNA-miRNA-lncRNA Network Construction of CASP8 in LIHC. A mRNA-miRNA-lncRNA network was constructed for finding out the potential regulatory axis of CASP8 in LIHC. The miRNA targets binding to CASP8 in LIHC were excavated by miRDB (http://www.mirdb.org/) [21], miRWalk (http:// mirwalk.umm.uni-heidelberg.de/) [22], and StarBase (https:// starbase.sysu.edu.cn/starbase2/index.php) [23] database. The common ground in these three databases was identified as the most significantly connected miRNAs. Next, according to the selected miRNAs, the lncRNA targets interacting to miRNAs were explored by StarBase and LncBase (https://diana.e-ce.uth .gr/lncbasev3) [24], and the association between the miRNAs and the lncRNAs was visualized. The expression level and prognostic performance of the targeted miRNAs and lncRNAs were also evaluated using TCGA-LIHC dataset.

Computational and Mathematical Methods in Medicine
analysis revealed these 26 PRGs mainly participated in pyroptosis, response to bacterium, positive regulation of interleukin-1 beta production, and cysteine-type endopeptidase activity involved in apoptotic process. The results of KEGG pathway analysis indicated these 26 PRGs were mainly associated with NOD-like receptor signaling pathway, pertussis, Epstein-Barr virus infection, apoptosis, and NF-kappa B signaling pathway.  There is no statistical difference between the overall survival curve of clus-ter1 and that of cluster 2 (P value = 0.057). Moreover, we also studied the association between clinical parameters and gene expression of cluster 1 and cluster 2 and found that age (>60 years or not) is a factor that influence the gene expression of cluster 1 and cluster 2 ( Figure 3(g), P value < 0.01).

The Correlation between Cluster and Immune Cell
Infiltration. To investigate the difference of the 2 clusters in tumor microenvironment of LIHC, we explored the relationship between immune cell infiltration and cluster. The StromalScore in cluster 1 was significantly higher than that in cluster 2 ( Figure 4 (Figure 4(d)).

Prognostic Performance and Prognostic Model
Construction of PRGs in LIHC. We also investigated the prognostic performance of above 26 PRGs. LIHC patients with high CASP3 (P value = 0.005), CASP4 (P value = 0.037), CASP6 (P value = 0.003), CASP8 (P value < 0.001), GPX4 (P value = 0.019), GSDMA (P value = 0.022), GSDME (P value < 0.001), NLRP3 (P value = 0.014), NLRP7 (P value = 0.045), NOD1 (P value = 0.005), NOD2 (P value = 0.001), PLCG1 (P value < 0.001), and SCAF11 (P value < 0.001) had worse OS than those with low expression of these genes, while LIHC patients with high NLRP6 expression had better OS than those with low NLRP6 expression ( Figure 5, P value = 0.014). No significant difference was found in OS between LIHC patients with high and low PRKACA, GSDMB, PJVK, CASP9, NLRP1, TIRAP, GSDMD, GSDMC, PYCARD, IL-1B, AIM2, and IL-6 expression (P value > 0.05). Thus, a total of 14 genes were consider as prognostic PRGs and were selected for further analyses. For better predicting the prognosis of LIHC patients, the 14 prognostic PRGs were selected to construct prognostic model by LASSO Cox analysis. The partial likelihood deviance and coefficients of prognostic model were presented in Figures 6(a) and 6(b). Risk score = ð0:0024Þ * CASP6 expression + ð0:0468Þ * CASP8 expression + ð0:0011Þ * GPX4 expression + ð0:0821Þ * GSDMA expression + ð0:0132Þ * GSDME expression + ð− 0:0240Þ * NLRP6 expression + ð0:0830Þ * NOD2 expression. LIHC patients were divided into high-risk and low-risk group according to the risk score. The OS plots of all LIHC cohort (P value < 0.001, Figure 6(c)), test cohort (P value =  14 14 14 14 14 13 13 12 8 7 7 8 7 7 4 Figure 6(d)), and training cohort (P value = 0.006, Figure 6(e)) revealed that high-risk group had worse OS than low-risk group with AUC of 0.738, 0.742, and 0.738, respectively. Figures 6(f)-6(h) illustrated the risk score distribution, patients' survival status, and gene expression level of low-risk and high-risk group in all LIHC cohort, test cohort, and training cohort, respectively. Next, univariate and multivariate analyses were performed to find out the potential factors that influenced the prognosis of LIHC patients. Evidently, risk score was considered as the independent factor affecting LIHC patients' prognosis in all LIHC cohort (Figures 7(a) and 7(b)), test cohort  3.6. Risk Score Associated with Clinical Parameters and Immune Cell Infiltration. We also focused on the influence of clinical parameters on risk score. The expression of the genes that constituted risk score and distribution of clinical characteristics in low-risk and high-risk group is summarized in Figure 9(a). Compared with cluster 1, LIHC patients in cluster 2 had higher risk score (Figure 9(d)

10
Computational and Mathematical Methods in Medicine risk score than those with grades I-II (Figure 9(e), P value = 0.0031), and LIHC patients with high ImmuneScore showed higher risk score than those with low Immune-Score (Figure 9(g), P value = 0.043). However, no significant difference was detected in the subgroups of gender, age, clinical stage, T stage, N stage, and M stage (Figures 9(b), 9(c), 9(f), and 9(h)-9(j), all P value > 0.05). In addition, the association between immune cell infiltration and risk score was also explored. However, only the correlation between abundance of activated memory CD4+ T cells and risk score was detected with a P value of 0.27 ( Supplementary Figure 1(a)). Thus, more studies are urgently needed to fill the gaps.

Immune Cell Infiltration, MSI, and TMB Analysis of CASP8 in LIHC.
We also constructed PPI network to visualize the association among the genes that constituted risk score ( Supplementary Figure 1(b)), which indicated that CASP8 was the core gene in the PPI network. Therefore, further analyses focusing on CASP8 in LIHC were performed.   (Figure 10(a)). Next, we explored the correlation between SCNAs of CASP8 and the infiltration levels of six immune cells in LIHC tissues. The infiltration levels of B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell were generally decreased when SCNAs appeared ( Figure 10(b)). Tumor mutation burden (TMB) and microsatellite instability (MSI) could assist clinicians to predict the efficacy of cancer immunotherapy. However, no significant correlation was detected between MSI score and CASP8 expression (Figure 10(c), P value = 0.509). In TMB analysis, similar result was presented in Figure 10(d) (P value = 0.745).
3.8. Construction of a mRNA-miRNA-lncRNA Network. The above analyses indicated that CASP8 was correlated with immune cell infiltration in LIHC. And CASP8 acted as the core gene in the PPI network constructed by PRGs associated with risk score. Thus, CASP8 was selected to construct a mRNA-miRNA-lncRNA network, which might uncover the potential CASP8-related regulatory axis in LIHC. First of all, CASP8 targeting miRNAs in miRDB, StarBase, and miRWalk database were compared to identify the common ground, which included hsa-miR-519a-3p, hsa-miR-105-5p, and hsa-miR-130b-3p (Figure 11(a)). Next, we explored the expression and prognostic value of these three genes in LIHC. The result indicated that only hsa-miRNA-130b-3p was differentially expressed in tumors and is significantly associated with patient prognosis.
To be more specific, the expression level of hsa-miRNA-130b-3p was significantly elevated in LIHC tissues compared with normal liver tissues ( Figure 11(b), P value = 7E-08), and LIHC patients with high hsa-miRNA-130b-3p expression had worse OS than those with low hsa-miRNA-130b-3p expression ( Figure 11(c), P value = 0.045). Thus, hsa-miRNA-130b-3p was considered as the most promising miRNA target of CASP8 in LIHC. The upstream lncRNA targets interacting to the promising miRNA targets were also explored via LncBase and StarBase database. A total of 3 lncRNA targets were identified based on the data in above databases, including H19, KCNQ1OT1, and MIR17HG ( Figure 11(d)). The expression and prognostic performance of these three genes in LIHC were also analyzed. Compared to normal liver tissues, the

Discussion
Pyroptosis is defined as a type of programed cell death that is dependent on the gasdermin protein family, and its occurrence is often as a result of activation of inflammatory CASP protein [25]. With the development of molecular biology, genetics, and medicine, the biological functions, mechanisms, and potential applications of pyroptosis are gradually uncovered. In the field of cancer research, some studies have summarized the important role of pyroptosis, which indicates that pyroptosis may influence the occurrence and progression of cancer in all stage [26,27]. And pyroptosis is considered as a promising mechanism in the cancer treatment due to the antiapoptotic effect of cancer cells [28]. However, the topic focusing on the role of PRGs in the prognosis and immune microenvironment of tumor is limited. Thus, we performed this bioinformatics analysis to systematically enucleate the expression, prognostic performance, and the effects on the immune microenvironment of PRGs in LIHC.
We firstly explored the expression level of 33 PRGs in LIHC tissues and normal liver tissues. Of these 33 PRGs, 26 genes with aberrant expression in LIHC tissues were selected to perform functional enrichment analysis and consensus cluster analysis. Programed cell death pathway, mainly including apoptosis and pyroptosis, protects mammals from infection. Outer membrane vesicles of Gramnegative bacteria delivered various bacterial molecules to host cells; outer membrane vesicle-associated molecules participated in the activation of apoptosis and pyroptosis [29]. CASP1 could cleave inactive IL-1 family to generate mature IL-1 family, such as IL-1β and IL-18, and GSDMD pore could mediate the release of IL-1 via electrostatic filtering [30]. He et al. reported that gene deletion of GSDMD interdicted the occurrence of pyroptosis and the secretion of IL-1β [31]. And the process of pyroptosis in various human cancer types was associated with NLRP3-related signaling pathways and NF-κB signaling pathway [12,[32][33][34]. To some extent, the above studies supported the result of functional enrichment analysis. However, more studies are required to validate the biological functions and potential mechanisms of PRGs in LIHC.
Prognostic model was widely utilized in prognostic evaluation of cancer patients, which provided epidemiological evidence for clinicians. In the current study, a total of 14 prognostic PRGs were enrolled to build up a prognostic model. Through further screening, the risk score consisted of the seven PRGs, including CASP6, CASP8, GPX4, GSDMA, GSDME, NLRP6, and NOD2. And the LIHC patients in high-risk score group generally showed worse OS compared to low-risk score group. Other three bioinformatics analyses, respectively, constructed prognostic model to explore the correlation between prognosis of lung  adenocarcinoma, skin cutaneous melanoma and glioblastoma, and PRGs, which indicated a worse prognosis in high-risk group [14][15][16]. The above results were similar with ours; nevertheless, large-sample multicenter studies are urgently needed to validate these results.
In the current study, CASP8 was detected as the core gene in the PPI network constructed by PRGs associated with risk score. Thus, we were interested on the role of CASP8 in LIHC. Previous studies reported that CASP8 not only induced apoptosis of LIHC cells but also had a nonapoptotic function in proliferation-related DNA damage response through CASP8/RIPK1/FADD/cFLIP complex in LIHC cells [35]. And Koschny et al. found that a stronger nuclear stain of CASP8 in LIHC cells and high nuclear CASP8 stain was associated with unfavorable prognosis after surgery and high tumor cell proliferation [36]. Thus, to find out the balance point of CASP8 in anticancer function and tumor-promoting action is of clinical significance.
A potential mRNA-miRNA-lncRNA regulatory axis in LIHC, lncRNA MIR17HG/hsa-miRNA-130b-3p/CASP8 regulatory axis, was also detected in the present study. MIR17HG participated in the generation of miRNA-17-92 cluster, which included miRNA-17, miRNA-18a, miRNA-19a, miRNA-20a, miRNA-19b, and miRNA-92a [37]. These six downstream miRNAs were detected increased expression levels in LIHC tissues and could promote invasion and proliferation of LIHC proliferation [38]. Wang et al. reported that miRNA-130b was upregulated in LIHC tissues in comparison with normal liver tissues, and high miRNA-130b was correlated with decreased survival, and overexpression of miR-130b promoted proliferation and metastasis of LIHC cells [39,40]. The above studies revealed the potential of this regulatory axis in LIHC cells. However, related in vivo and in vitro researches are urgently required.
Some limitations still exist in the current study. Firstly, most of the analyses in the current study are conducted at the transcription level; some results may not apply to the studies based on the protein level. Secondly, related fundamental and clinical studies that focus on the clinical value and molecule mechanisms of PRGs in LIHC are rare. Third, genetic background and etiology of LIHC patients are influenced by many factors, such as patients' race, patients' gender, and patients' age. Thus, more in-depth studies are necessary to validate the results.

Conclusion
In summary, the present study built up a pyroptosis-related gene signature to predict LIHC patients' prognosis and their correlations with immune infiltration, which indicated that pyroptosis-related genes were of great significance in the LIHC patients' prognosis and microenvironment of LIHC. However, in-depth studies are needed to validate our results.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.