Pyroptosis-related genes IRAK1, NOD2, POP1 and YWHAB impact the prognosis through immune functional pathways in hepatocellular carcinoma

Background: Hepatocellular carcinoma (HCC) is a highly heterogeneous disease, which makes the prognostic prediction challenging. Pyroptosis can inuence the proliferation, invasion and metastasis of tumors. However, the prognostic value of pyroptosis-related genes in HCC remains to be further elucidated. Methods: In this study, univariate Cox regression analysis and the least absolute shrinkage and selection operator Cox regression were performed to construct a four-gene signature. Kaplan-Meier curve, time-dependent receiver operating characteristic curve and univariate/multivariate Cox analysis were used to evaluate the prognostic model. Single-sample gene set enrichment analysis was performed to calculate the enrichment scores of diverse immune cell subpopulations and related functions. The correlation of the prognostic model with immune inltrate types and drug chemosensitivity were analyzed by Analysis of variance and Pearson correlation respectively. Spearman correlation was carried out to analyze the associations between the risk score and immune checkpoint-related genes, cell cycle-related genes, and multidrug resistance-related genes. Gene set enrichment analysis was conducted to further explore the underlying molecular mechanisms. Finally, real-time quantitative polymerase chain reaction and immunohistochemistry were applied to analyze the mRNA and protein expression levels between HCC tissues and adjacent non-tumorous tissues in an independent sample cohort. Results: A pyroptosis-related gene signature was constructed to classify patients into two risk groups. Compared with the low-risk group, patients in the high-risk group showed signicantly reduced overall survival (OS). The risk score was an independent predictor of OS. Tumor-related pathways were enriched in the high-risk group and patients with different risk scores showed different immune status. Immune checkpoint-related genes, cell cycle-related genes, and multidrug resistance-related genes were overexpressed in the high-risk group compared with the low-risk group. The expression level of each prognostic gene was negatively related to the anti-tumor drug sensitivity. Furthermore, the expression level of each prognostic gene in HCC tissues was higher than that in adjacent non-tumorous tissues in an independent sample cohort and similar results were found in most cancer types.


Background
Hepatocellular carcinoma (HCC) is one of the most frequently diagnosed malignant tumors and ranks the fth highest incidence and third highest mortality rate of all malignancies worldwide [1]. Of all liver cancer cases, more than 90% are HCCs [2]. Risk factors include chronic hepatitis B virus and hepatitis C virus infections, alcohol abuse, nonalcoholic fatty liver disease, and exposure to dietary toxins such as a atoxins [3]. Patients are often diagnosed with HCC in advanced stages, which is the reason for the poor prognosis [2]. Also, unclear molecular mechanisms lead to the poor prognosis of HCC. The overall survival (OS) of HCC varies across the world, with a 5-year survival rate of 18% in the United States, and only 12% in China [4]. Therefore, it is critically important to clarify the molecular mechanisms underlying the poor prognosis of HCC and explore new prognostic biomarkers for HCC.
Apoptosis as the signi cant mechanism of anti-cancer defense has been extensively and deeply studied, but the relationship between pyroptosis and cancer remains unclear [5]. From the current data, we learn that drug resistance to apoptosis is a possible cause of cancer treatment failure, and introducing nonapoptotic programmed cell death, such as pyroptosis, may be an effective way to overcome apoptosisresistant cancers [14]. Previous studies have found that pyroptosis can in uence the proliferation, invasion and metastasis of tumors [5]. Induced pyroptosis of cancer cells can eliminate malignant cells and suppress the occurrence and development of cancers [15]. Also, accumulating studies have proved that pyroptosis-related genes play an important role in HCC such as caspase-1 [16,17], caspase-3 [17] and GSDMD [17]. However, it is still largely unknown whether these pyroptosis-related genes are associated with the prognosis of HCC. Furthermore, the relationship between these genes and the immune status during the progression of HCC deserves further study.
In the current study, a prognostic four-gene signature was developed based on pyroptosis-related genes.
To gain insight into the four-gene signature, we discussed its relationship with survival and clinicalpathological characteristics, immune in ltration, immune pathways, immune checkpoint-related genes, cell cycle-related genes, multidrug resistance-related genes, and chemosensitivity. Besides, we examined the mRNA and protein expression levels of prognostic genes between HCC and adjacent non-tumorous tissues in an independent sample cohort. It is hoped that this robust model might provide new guidance and insights for the prognosis and treatment of HCC.

Data collection from TCGA and ICGC
The mRNA sequencing data of HCC and correlated clinical information were obtained from public databases. The data of 370 HCC patients were obtained from The Cancer Genome Atlas hepatocellular carcinoma (TCGA-LIHC) portal (https://portal.gdc.cancer. gov/repository). Additional 231 tumor samples were collected from the International Cancer Genome Consortium hepatocellular carcinoma (ICGC-LIRI-JP) portal (https://dcc.icgc.org/projects/LIRI-JP). We collected the gene expression and clinical characteristics of each HCC patient in TCGA and ICGC cohorts (see Original Data le 1-4). The present study follows the access policies and publication guidelines. Then, 69 pyroptosis-related genes were obtained from the previous literature [5,12,18,19,20,21] (Supplementary Table 1).
Identi cation of DEGs between HCC and adjacent non-tumorous tissues With the help of the "limma" R package, the differentially expressed genes (DEGs) between HCC and adjacent non-tumorous tissues were screened with a false discovery rate (FDR) < 0.05. PPI networks of DEGs were constructed with the STRING database.
Establishment and validation of a peroxisome-related prognostic model Univariate Cox analysis of OS was performed to identify pyroptosis-related genes among the DEGs with prognostic values. Subsequently, the prognostic gene signature was constructed by the least absolute shrinkage and selection operator (LASSO) penalized Cox regression analysis, and the optimal values of the penalty parameter were determined through 10-times cross-validations by using the R package "glmnet". The function of the LASSO was to select and shrink variables, and some regression coe cients were strictly equal to 0, thereby obtaining an interpretable model. The normalized expression matrix of candidate prognostic DEGs was the independent variable in the regression, and OS as well as the status of patients were the response variables. The risk score of each patient was calculated according to the expression level of each pyroptosis-related gene and its corresponding regression coe cients. The formula was established as follows: score = e sum (each gene' s expression × corresponding coe cient) . Then, we classi ed patients into high-and low-risk groups for OS according to the median value of the risk score determined by the "survminer" package. With the help of the "Rtsne" R package, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed to explore the distribution of different groups based on the expression levels of the prognostic genes. Receiver operating characteristic (ROC) curves were plotted and the area under the curve (AUC) values were calculated using the R package "survivalROC".

Identi cation of Independent prognostic factors for OS in HCC
To evaluate the predictive power of the risk score and other clinical risk factors, including age (< 60, ≥ 60), gender (male, female), tumor grade (G1/2, G3/4) and tumor stage ( / , / ). We rst performed univariate Cox regression analysis and multivariate Cox regression analysis to identi ed independent prognostic factors. Then we calculated the AUC of each independent prognostic factor to compare the prognostic values of clinical risk factors and the risk score.

Evaluation of immune status in different risk groups
With the help of the "GSVA" R package, single-sample gene set enrichment analysis (ssGSEA) was performed to quantify the enrichment scores of diverse immune cell subpopulations and related functions. Six immune subtypes were de ned to measure the type of immune in ltration in the tumor microenvironment and the correlation between the risk score and immune in ltration subtypes was performed by Analysis of Variance. Our study showed the immune subtypes of cancer patients in the TCGA cohort (see Original Data le 5). Besides, the Wilcoxon test was performed to analyze the expression levels of immune checkpoint-related genes (PD-L1, PD-L2 and Galectin-9) in different risk groups and the association between the risk score and immune checkpoint-related genes was analyzed by Spearman correlation.
Cell cycle-related genes, multidrug resistance-related genes and chemotherapy sensitivity analysis The NCI-60 database containing 60 different cancer cell lines from 9 different types of tumors was accessed through the CellMiner interface (https://discover.nci.nih.gov/cellminer/). Spearman correlation was used to analyze the correlation of the prognostic model with cell cycle-related genes and multidrug resistance-related genes. In our study, a number of 263 drugs, which were on clinical trials or approved by the FDA, were used to evaluate the correlation between prognostic gene expression and drug sensitivity by Pearson correlation (Supplementary Table 2).

KEGG Functional Enrichment Analyses
To detect potential biological functions and involved signaling pathways in different risk groups, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed with gene set enrichment analysis (GSEA) by GSEA software 4.1.
Real-time quantitative-polymerase chain reaction (qRT-PCR) 20 pairs of HCC and adjacent non-tumorous tissues were recruited from the First A liated Hospital of Wenzhou Medical University, and this study was approved by the Review of Ethics Committee in Clinical Research of the First A liated Hospital of Wenzhou Medical University. RNA was extracted from the HCC and adjacent normal liver tissue specimens using the Trizol reagent according to the manufacturer's instructions (Servicebio). Then, the total RNA was reverse-transcribed into cDNA with RevertAid First Strand cDNA Synthesis Kit (Thermo). GAPDH levels were used to normalize all gene expression levels. PCR was performed using FastStart Universal SYBR Green Master (Roche) on the ABI StepOne (Applied Biosystems). Primers are listed in Supplementary Table 3. Each PCR reaction was repeated in triplicate for stable results. Gene expression level was quanti ed using the 2 −ΔΔCt method. We obtained the expression data of each prognostic gene from qRT-PCR (see Original Data le 6) Immunohistochemistry (IHC) IHC experiments were performed to validate the protein expression levels of the four prognostic genes in 10 paired HCC and adjacent non-tumorous tissues. All specimens were xed in 10% neutral-buffered formalin at room temperature. After embedded in para n, a series of 4-μm sections were sliced. Sections were depara nized, followed by rehydration and then a 10-min boiling in 10 mmol/L citrate buffer (pH = 6.0). Endogenous peroxidase was blocked by immersing the sections in a 3% solution of hydrogen peroxide in methanol for 10 min. To block the unspeci c binding, sections were incubated with 1% bovine serum albumin in phosphate-buffered saline for 30 minutes. After that, slices were stained by primary antibodies and HRP-conjugated secondary antibodies. Antibodies are detailed in Supplementary Table 4. Then, the slices were stained with 3,3'-diaminobenzidine and counterstained with hematoxylin. In the end, after dehydration, the samples were sealed, observed and photographed. Our study showed the represent pictures of each prognostic gene from IHC (see Original Data le 7-14).

Statistical analysis
Wilcoxon test was performed to compare gene expression between tumor and adjacent non-tumorous tissues. The Chi-squared test was used to compare the differences between categorical variables. The difference in OS between the two groups was evaluated by Kaplan-Meier analysis and the log-rank test. Univariate and multivariate Cox regression analyses were implemented to screen the signi cant factors affecting OS. Mann-Whitney test was performed to compare the ssGSEA scores of immune cells and the activity of immune-related pathways in different groups. Spearman correlation was carried out to analyze the associations between the risk score and immune checkpoint-related genes, cell cycle-related genes, and multidrug resistance-related genes. Pearson correlation was utilized to analyze the relationship of the prognostic gene expression to drug sensitivity. R software (version 4.0.2) and SPSS software (version 23.0) were used for all statistical analyses and plots. If not speci ed above, a two-tailed P value < 0.05 was considered statistically signi cant in this study.

Results
The overall data processing ow is shown in Figure 1. In our study, 365 HCC patients were provided by the TCGA-LIHC cohort and 231 HCC patients were provided by the ICGC (LIRI-JP) cohort. Clinical data for these patients are presented in Table 1.
Identi cation of pyroptosis-related prognostic DEGs in the TCGA cohort More than half of the pyroptosis-related genes (47/69, 68.1%) were expressed differently between tumor and adjacent non-tumorous tissues. Among them, twenty-one genes were associated with OS in univariate Cox regression analysis ( Figure 2A). A total of 21 pyroptosis-related DEGs were preserved in our study (all FDR < 0.05) ( Figure 2B-C). The interaction network and correlation among these genes are presented in Figure 2D-E.
Construction of a four-gene signature in the TCGA cohort A total of 21 survival-related DEGs were mentioned above, then LASSO Cox regression analysis identi ed four genes to construct the prognostic model (Supplementary Figure 1). The four genes included in the model were IRAK1, NOD2, POP1 and YWHAB. The risk score = 0.101 * expression level of IRAK1 + 0.003 * expression level of NOD2 + 0.164 * expression level of POP1 + 0.438 * expression level of YWHAB. The patients were strati ed into a high-risk group and a low-risk group according to the median cut-off value ( Figure 3A). Based on PCA and t-SNE analyses, patients in different risk groups were distributed in two directions ( Figure 3E-F). In the TCGA cohort, we found that the risk scores were increased with the increase of tumor grade and TNM stage ( Table 2). Patients with high risk were more likely to die earlier than those with low risk ( Figure 3B). Conformably, the Kaplan-Meier curve revealed that the OS was signi cantly poorer in the high-risk group than in the low-risk group (P < 0.001) ( Figure 3I). Timedependent ROC analysis was used to evaluate the predictive power of the four-gene signature and the AUC for 1-, 2-, and 3-year OS were 0.735, 0.666, and 0.663, respectively ( Figure 3J). Survival analysis con rmed that the high expression of IRAK1, NOD2, POP1 and YWHAB was correlated with the poor prognosis of HCC (P < 0.05) (Supplementary Figure 2). Furthermore, the expression of each prognostic gene was higher in HCC tissues compared with adjacent non-tumorous tissues (Supplementary Figure 3).

Validation of the four-gene signature in the ICGC cohort
The prognostic model was validated in the ICGC cohort. Patients were divided into a high-risk group and a low-risk group based on the median value from the TCGA cohort ( Figure 3C). Similar to the above results, the patients in different risk groups were distributed in discrete directions based on PCA and t-SNE analyses ( Figure 3G-H). Besides, the risk scores were increased with the increase of the TNM stage (Table  2). Equally, compared with the low-risk group, patients in the high-risk group had a shorter survival time ( Figure 3D) and poorer OS (P < 0.05) ( Figure 3K). Finally, the AUC of the four-gene signature for 1-, 2-, and 3-year OS were 0.638, 0.636, and 0.638, respectively ( Figure 3L).
Furthermore, the combination of the risk score with the tumor stage could provide a more accurate prediction of 3-year OS for HCC (TCGA cohort: AUC = 0.718; ICGC cohort: AUC =0.766) ( Figure 5).

Risk scores and prognostic genes in different clinical characteristics groups
According to the risk score and clinical characteristics of HCC patients in the TCGA cohort, we discovered that the risk score was higher in the advanced tumor grade and advanced tumor stage ( Figure   6C-D). The results obtained from the ICGC cohort were similar to the TCGA cohort, but there was no data about the grade of LICH in the ICGC cohort ( Figure 6E-G). Combining the expression levels of prognostic genes with the clinical characteristics of HCC revealed that the expression levels of NOD2, POP1 and YWHAB were different in tumor grade (P < 0.05), and the expression levels of POP1 as well as YWHAB were different in the tumor stage (P < 0.05) in the TCGA cohort (Supplementary Figure 4C-D). Moreover, in the ICGC cohort, IRAK1 was distinguishingly expressed in the tumor stage (Supplementary Figure 4G).

Relationship between risk score and immune status
To determine the immune status in different risk groups, ssGSEA was performed to quantify the enrichment scores of diverse immune cell subpopulations and related functionsIn both the TCGA and the ICGC cohorts, the immune cell subpopulations including aDCs, iDCs, macrophages, Tfh as well as Treg showed high in ltration in the high-risk group, but NK cells were opposite (P < 0.05) ( Figure 7A-B). Moreover, the scores of immune-related functions including APC co-stimulation, CCR, Check point, HLA, MHC class I and parain ammation were signi cantly elevated in the high-risk group, but the score of type Interferon (IFN) response was opposite (P < 0.05) ( Figure 7C-D).
To further understand the role of the four-gene signature in HCC, we explored the expression levels of immune checkpoint-related genes (PD-L1, PD-L2 and Galectin-9) in different risk groups. We found that the expression levels of PD-L1, PD-L2 and Galectin-9 were signi cantly higher in the high-risk group than in the low-risk group and they were all positively correlated with the risk score (P < 0.05) ( Figure 8A-B).
The correlation between risk score and different immune in ltrate subtypes in HCC patients was conducted by the Analysis of Variance. Six types of immune in ltrates were found in human tumors. They are C1 (wound healing), C2 (IFN-r dominant), C3 (in ammatory), C4 (lymphocyte depleted), C5 (immunologically quiet) and C6 (TGF-β dominant) [22]. No patient sample belonged to the C5 immune subtype and only one patient sample belonged to the C6 immune subtype, so we excluded them. From C1 to C4 subtype, the value of the risk score showed a downward trend, and risk scores between C1 and C3, C4, also between C2 and C3, C4 were signi cantly different (P < 0.05) ( Figure 9). Besides, the expression levels of prognostic genes were signi cantly different in immune in ltrate subtypes. Most of the prognostic genes (IRAK1, POP1 and YWHAB) were highly expressed in the C1 subtype and NOD2 was highly expressed in the C2 subtype (P < 0.05) (Supplementary Figure5).

Differential expression of cell cycle-related genes in risk groups
To understand the correlations between the risk score and the expression levels of cell cycle-related genes, Spearman correlation was carried out. CDK2, CDK4, CDC25A and Cyclin E were overexpressed in the high-risk group compared with the low-risk group. Furthermore, the expression levels of CDK2, CDK4, CDC25A and Cyclin E were all positively correlated with the risk score (P < 0.001) ( Figure 10A-B).

Analysis of the correlation between the risk model and multidrug resistance-related genes and chemotherapeutics
To assess the impact of the four-gene signature on HCC treatment, we analyzed the correlations of the risk score with the expression levels of multidrug resistance-related genes. The results indicated that the expression levels of MRP1, MRP4 and MRP5 were higher in the high-risk group than in the low-risk group. Also, the expression levels of MRP1, MRP4 and MRP5 were positively correlated with the risk score ( Figure 11A-B). Furthermore, the correlation analysis was conducted to explore the correlation between the risk score and chemotherapeutics. We found that the expression levels of prognostic genes were negatively correlated with the sensitivity of some chemotherapy drugs ( Figure 12). For example, the sensitivity of Bosutinib was negatively correlated with the expression of IRAK. More details are listed in supplementary table 5.

KEGG functional enrichment analysis
The GSEA was used to perform KEGG analysis in different risk groups to explore the potential molecular mechanisms of the four-gene signature. According to the results, we found that some pathways were enriched in the high-risk group in both the TCGA and ICGC cohorts ( Figure 13). Of note, some pathways were related to cancer processes such as cell cycle, MAPK signaling pathway, mismatch repair, mTOR signaling pathway, Notch signaling pathway, pathways in cancer, VEGF signaling pathway and Wnt signaling pathway. The detailed information is displayed in supplementary gure 6-7.

Veri cation of prognostic gene expression levels between HCC and adjacent non-tumorous tissues by qRT-PCR and IHC
We analyzed mRNA and protein expression levels by qRT-PCR and IHC to validate the different expression levels of four prognostic genes (IRAK1, NOD2, POP1 and YWHAB) in HCC and adjacent non-tumorous tissues. The qRT-PCR results indicated that the expression levels of IRAK1, NOD2, POP1 and YWHAB in HCC tissues were higher than those in adjacent non-tumorous tissues ( Figure 14A). The same results were shown by IHC ( Figure 14B). Furthermore, the results were consistent with RNA sequencing results of the four prognostic genes in the TCGA cohort (Supplementary Figure 3). Consistently, we compared the expression levels of four prognostic genes in normal and tumor samples across different TCGA cancer types ( Figure 15). The result indicated that prognostic genes in tumor samples were highly expressed in most cancer types compared with normal samples.

Discussion
HCC, one of the most common cancers, is a leading cause of cancer-related deaths worldwide [23] and has a poor prognosis due to its underlying complex mechanisms. Consequently, it is critically important and urgent to nd an effective prognostic prediction method and treatment strategy for HCC.
In the current study, a total of 47 overlapping DEGs between HCC tissues and non-tumor tissues were identi ed from the TCGA cohort, of which 21 were associated with OS. Then, a novel four-gene signature was constructed by the LASSO algorithm in the TCGA cohort and validated in the ICGC cohort. The risk score was calculated based on the above model, which strati ed HCC patients into high-or low-risk groups. Furthermore, the risk score was an independent prognostic factor of OS.
The prognostic model was composed of four pyroptosis-related genes (IRAK1, NOD2, POP1 and YWHAB) and all of them were upregulated in HCC tumor tissues and were correlated with poor prognosis in our study. Previous studies have revealed that IRAK1 is closely related to a variety of malignant cancers [24], such as HCC [25] and lung cancer [26]. Also, Cheng et al. [25] claimed that IRAK1 augments HCC stemness and drug resistance via the AP-1/AKR1B10 signaling pathway [25]. NOD2 is a member of the NOD-like receptor family [27]. A study by Zhou et al. showed that hepatic NOD2 promotes hepatocarcinogenesis and it is with strong linkage to the poor prognosis of HCC [28]. POP1, as a new candidate gene, may be responsible for a peculiar type of skeletal dysplasia with severe short stature [29]. An additional study on POP1 showed that POP1 is associated with human prostate cancer [30]. YWHAB encodes a number of 14-3-3 family proteins, of which 14-3-3β regulates multiple signaling pathways in normal and cancer cells [31]. Moreover, up-regulation of the 14-3-3β enhances HCC cell migration and proliferation [32,33]. Collectively, most of the prognostic genes in the model have been reported to be involved in the occurrence and development of cancers and are associated with HCC, suggesting that these genes have the potential to be used as prognostic biomarkers for HCC. Nevertheless, whether these genes affect the prognosis of HCC remains to be elucidated.
Pathway regulation plays an important role in the occurrence and development of tumors [34,35]. Here, pathways such as MAPK signaling pathway, Notch signaling pathway, P53 signaling pathway, VEGF signaling pathway and Wnt signaling pathway were enriched in the high-risk group and all of them have been proved to be involved in the cancer process [36, 37,38]. Based on this evidence, one can speculate that four genes in our model may promote the carcinogenesis and development of HCC through the above biological pathways.
The cell cycle is a sophisticated and tightly controlled process, commonly divided into G0/G1, S, G2, and M phases. These ordered events are fundamentally governed by multiple CDKs and cyclins [39]. Errors occurring at any step of the cell cycle can lead to apoptosis or cause genetic instability promoting cancer aggressiveness [40]. Interestingly, cell cycle-related genes (CDK2, CDK4, CDC25A and Cyclin E) were aberrantly up-regulated in the high-risk group and aberrant up-regulation of CDK2, CDK4, CDC25A and Cyclin E was positively correlated with the risk score. Of note, the cell cycle pathway was also enriched in the high-risk group (Fig. 13). Hence, it is reasonable to assume that the four-gene signature may drive the occurrence and development of HCC by changing the cell cycle.
Due to the heterogeneity of the different epigenetic and genetic backgrounds in tumor subtypes, patients with similar clinical characteristics, such as tumor stage, may have different outcomes [41]. Our study identi ed that the risk score was an independent prognostic factor of OS according to univariate and multivariate Cox regression analyses, and it was comparable or superior to the clinical characteristics (age, gender, tumor grade, and tumor stage) of HCC at predicting OS. We also analyzed the relationship between the risk score and different clinical statuses and found the risk score was increased with the increase of tumor stage and grade. Combined with the tumor stage, the risk score performed better in predicting OS. Therefore, the combined use of the four-gene signature and clinical-pathological features might be bene cial to identify high-risk patients and predict the prognosis of HCC.
According to our ndings, most of the prognostic genes were highly expressed in the more aggressive subtype of immune in ltration (C1), suggesting that the immune in ltrate subtype may be one of the factors affecting the prognosis of HCC. Besides, the high-risk group has higher fractions of macrophages and Treg cells. Accumulating evidence indicated that macrophages and Treg cells promote the growth and metastasis of HCC [42,43]. Moreover, low levels of NK cells [44,45] and type II IFN response [46,47] are associated with impaired antitumor immunity. Further considering the results of our study, we suspected that attenuated antitumor immunity might be an important reason for the poor prognosis of HCC. Unexpectedly, we discovered that aDCs, iDCs, APC co-stimulation, HLA and MHC class I were enriched in the high-risk group, and all of them are related to antigen presentation, suggesting the prognostic genes in our study may change the tumor microenvironment and immune status through affecting antigen presentation, but speci c mechanisms warrant further investigation. It is worth noting that immune checkpoint-related genes (PD-L1, PD-L2 and Galectin-9) were aberrantly up-regulated in the high-risk group and the overexpression of PD-L1, PD-L2 and Galectin-9 was positively correlated with the risk score. During tumor occurrence and development, the effective anti-tumor immune surveillance in the liver microenvironment is impaired, and immune checkpoints, especially the PD-1/PD-L1 signaling pathway, play an important role in the said process [48]. Moreover, tumor-associated PD-L1 expression is predictive of clinical response to PD-1-directed immunotherapy [49]. Therefore, we can predict the expression levels of immune checkpoints in HCC patients through the four-gene signature, thereby judging the sensitivity of HCC patients to immunotherapy in clinical treatment.
During the long-term of traditional chemotherapy for HCC, multiple drug resistance occurs frequently, leading to the relapse of cancer and intractable tumor [50]. Our study indicated that multidrug resistancerelated genes, such as MRP1, MRP4 and MRP5, were highly expressed in the high-risk group compared with the low-risk group. Moreover, the expression levels of prognostic genes were negatively correlated with the sensitivity of certain chemotherapy drugs, such as epirubicin and uorouracil, suggesting such aberrant up-regulation of prognostic genes could promote HCC drug resistance. Given these ndings, it is reasonable to assume that the four-gene signature may contribute to the cancer cell resistance to chemotherapy treatment. Remarkably, some chemotherapy drugs, including epirubicin [51] and uorouracil [52], can be used in the clinical treatment of HCC, so we have reason to believe that the fourgene signature may have important clinical application value in the selection of chemotherapy drugs sensitive to HCC patients.
In addition, we found that the expression levels of prognostic genes in HCC tissues were higher than those in adjacent non-tumorous tissues in an independent sample cohort. Intriguingly, the prognostic genes were also signi cantly up-regulated in most cancer types, suggesting that the four-gene signature may have a universal role in predicting the prognosis of pan-cancers.
Nevertheless, the speci c underlying mechanisms between pyroptosis-related genes and tumor immunity in HCC remain poorly understood and need extensive study.

Conclusion
Our study established a four-gene signature, which can be used to predict the prognosis of HCC patients. In addition, we evaluated the four-gene signature under various clinical settings, such as survival and clinical-pathological characteristics, immune in ltration, immune pathways, immune checkpoint-related genes, cell cycle-related genes, multidrug resistance-related genes, and chemosensitivity. Finally, laboratory experiments were performed to examine the mRNA and protein expression levels of the four pyroptosis-related genes. In summary, the four-gene signature can provide new insights into the prognosis and immune of HCC and points out a possible direction for individualized treatment of HCC in the future.

Declarations
Ethics approval and consent to participate All procedures performed in this study involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. This study was approved by the Review of

Consent for publication
Not applicable Availability of data and materials All data generated or analyzed during this study are included in this published article and its supplementary information les.

Competing interests
The authors declare that they have no competing interests.         Comparison of the risk scores between different immune in ltrate subtypes