Development and validation of a prognostic model related to pyroptosis-related genes for esophageal squamous cell carcinoma using bioinformatics analysis.

Background: Esophageal squamous cell carcinoma (ESCC) is one of the most lethal malignant tumors worldwide, and a larger number of ESCC patients have unsatisfactory overall survival (OS) rates. While pyroptosis participates in the development of a variety of malignancies, the function of pyroptosis-related genes (PRGs) in ESCC is still obscure. The aim of this study was to construct the pyroptosis-related prognostic model for ESCC, which will be developed to stratify the risk hazards of ESCC patients and to provide theoretical evidence for individualized treatment.
Methods: RNA-seq data of ESCC were download from the NCBI Gene Expression Omnibus (GEO) database. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to explore the potential biological functions or pathways. OS was considered as the primary prognosis outcome in this study. The riskscore was constructed by Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis. The pyroptosis-related prognostic model was constructed based on all independent prognostic factors and verified by C-index, Receiver operating characteristic (ROC) curves, and Calibration curves, and the role of the riskscore in ESCC immunotherapy was evaluated by the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm.
Results: The current study found 31 differentially expressed PRGs (P<0.001), and functional enrichment analysis showed these PRGs were enriched in positive regulation of cytokine production, interleukin-1 beta production. Univariate and multivariate Cox regression analysis were applied to validate that the riskscore based on four prognostic PRGs (HMGB1, IL-18, NLRP7, and PLCG1) was an independent prognostic factor for ESCC, and the C-index of prognostic model related to the riskscore (C-index =0.705) was higher than that of tumor node metastasis (TNM) stage (0.620). The low-risk group showed a better efficacy of immune checkpoint inhibitors. Conclusions: The riskscore related to PRGs was one of the independent prognostic factors for ESCC. Moreover, the prognostic model related to the riskscore could be used to predict the OS of ESCC patients effectively. However, there still were several limitations in this study, such as no external validation sample. In summary, our data provides a novel perspective in exploring the potential prognostic biomarkers of ESCC. 2022 Journal of Thoracic Disease. All rights reserved.

Introduction

Esophageal cancer (EC) is the 7th most prevalent malignant neoplasm and 6th leading cause of cancer-related death worldwide (1). Pathological types of EC are classified as esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC), of which ESCC accounts for 90% of all cases of EC and mostly occurs in the Eastern Asian and Eastern and Southern African regions (2,3). The etiology of ESCC is incompletely understood, and may include genetic factors (e.g., mutation), environmental factors (e.g., alcohol and tobacco), viral factors [e.g., papillomavirus (HPV) infection], or the coexistence of these factors (4). When the esophageal mucosa is exposed to repeated damage from mechanical injury or carcinogens, abnormal epithelium cell proliferation occurs, and the abnormal cell will eventually develop into invasive cancer. Owing to the large population and high incidence of ESCC, the number of cases in China makes up more than half of those worldwide (5). Although patients with ESCC are treated according to current standard guidelines, some with locally advanced ESCC have a unsatisfactory overall survival (OS) (6), with the global 5-year OS rate around 15–25% (1,7). Current treatment options for ESCC generally depend on its stage (8-10), with endoscopic resection for mucosal lesions versus esophagectomy for submucosal lesions the primary treatment option for early disease. Extensive studies have recorded that neoadjuvant therapy combined surgery is beneficial in improving the prognosis of locally advanced ESCC (11). Over the past years, technology advancements have comprehensively characterized the genomic landscape in ESCC, and targeted therapy has attracted increased attention (12). To the best of our knowledge, pembrolizumab, ramucirumab, and trastuzumab have been approved by the U.S. Food and Drug Administration (FDA) for the treatment of the advanced EC (13). Despite this, a minority of targeted therapies have prognostic and therapeutic implications for patients with ESCC, and the development of biomarkers and therapeutic targets is required.

Pyroptosis is triggered by caspase-1/4/5/11 and is one type of lytic programmed cell death (14). Differences between pyroptosis and apoptosis in morphology were observed by Zychlinsky et al. in 1992 (15), and Cookson et al. created the term “pyroptosis” after observing macrophages infected with bacteria in 2001 (16). Pyroptosis can cause cell swelling, plasma membrane lysis, and chromatin fragmentation, which lead to osmotic swelling and cell death (17). Moreover, it can cause the release of intracellular pro-inflammatory factors, including IL-18 and IL-1β, into the extracellular space (18). Pyroptosis shows several similar features to apoptosis, including DNA damage, nuclear pyknosis, and TUNEL-positive staining, but with several unique characteristics (19,20). The mechanism of pyroptosis can be summarized into the following three pathways: the canonical inflammasome pathway, the non-canonical inflammasome pathway, and the caspase-3-dependent pyroptotic pathway (21,22). Triggered caspase 1/4/5/11/3 lead to gasdermin D (GSDMD) cleavage into the gasdermin-C domain and gasdermin-D domain, and the gasdermin-D domain then combines with the acidic phospholipids of the membrane to form oligomeric death-inducing pores, called gasdermin pores (23-26). This type of pore in membrane leads to cell swelling, rupture of the membrane, and eventual cell death (27). Pyroptosis is considered a general innate immune response (18). It can defend against intracellular infection, lead to the death of intracellular bacteria through pore-induced intracellular traps (28,29), and eliminate the pathogen protective niche (27). Pyroptosis is also significantly linked to the tumorigenesis and development of various cancers. The NLR family pyrin domain containing 3 (NLRP3), one pyroptosis-related molecule, was significantly down-regulated in hepatocellular carcinoma (HCC) tissue and associated with prognosis of HCC patients (30). According to a published report, the inhibition of caspase-1 interfered with the anticancer activity of euxanthone in HCC cells (31). Dihlmann et al. have demonstrated that a low expression of the interferon-inducible protein AIM2 was found in colorectal cancer (CRC) cells, and the mortality of patients with low AIM2 expression levels was higher than patients with normal levels (32). A study has shown docosahexaenoic acid (DHA) promoted the activation of caspase-1 and gasdermin D in breast cancer cells, suggesting it induces pyroptosis, and the anti-cancer effect of DHA in breast cancer was associated with it (33).

Pyroptosis-related molecules have also been confirmed as potential biomarkers in various kinds of human malignant neoplasms, including cervical cancer, ovarian carcinoma (OC), non-small cell lung cancer (NSCLC), and gastric cancer (GC) (34-37). Ye et al. developed a prognostic model related to pyroptosis for OC for the first time, which demonstrated good OS prediction (38), and while a LUAD prognostic model constructed using five pyroptosis-related genes (PRGs) was able to predict the OS of LUAD patients with moderate to high accuracy (39). To date, several prognostic biomarkers and associated prognostic model have been suggested in ESCC. A prognostic model based on six N6-methyladenosine (m6A) RNA methylation regulator (METTL3, WTAP, IGF2BP3, YTHDF1, HNRNPA2B1 and HNRNPC) was able to predict the OS of ESCC patients (40). Yao et al. developed a prognostic model of ESCC associated with 10 M2 macrophage genes, which also had good predictive power for OS (41). However, there were several limitations for these studies. These prognostic models need to be further validated in several large cohorts and multi-center clinical trials. In addition, further experiments are needed to elucidate the potential mechanisms of these associated genes in the progression of ESCC.

To the best of our knowledge, no prognostic model related to ESCC has been developed. Herein, we performed a bioinformatics analysis to identify the four prognostic PRGs in ESCC. The riskscore based on the PRGs was constructed and validated in a validation set, and the prognostic model related to the riskscore was also constructed and validated. Finally, we explored the potential role of the riskscore in guiding the immunotherapy of ESCC patients. Our data provided a novel prognostic model that can be used to stratify the risk hazards of patients with ESCC, and to offer theoretical evidence for individualized therapy of ESCC patients. We present the following article in accordance with the TRIPOD reporting checklist (available at https://jtd.amegroups.com/article/view/10.21037/jtd-22-948/rc).

Methods

Research design

All analyses performed are shown in . Firstly, we selected expression profiles data of PRGs and performed differential analysis to obtain differential PRGs. Gene Ontology (GO)-terms Semantic Similarity Measures (GSSM), GO, Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses and protein-protein interaction (PPI) network construction were performed based on differential PRGs. Next, we obtained the training set and the validation set after randomly grouping the ESCC cases. In the training set, we found four prognostic PRGs by cox regression analysis and constructed the riskscore based on the four prognostic PRGs via Least Absolute Shrinkage and Selection Operator (LASSO) cox regression analysis. Subsequently, the riskscore was validated in the validation set, and the nomogram was established based on the riskscore and other independent prognostic factors. Lastly, the correlation between riskscore and the response to immunotherapy was explored via Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Figure 1

Flow diagram of all data analysis process in this study. GSSM, GO-terms Semantic Similarity Measures; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; LASSO, Least Absolute Shrinkage and Selection Operator.

Acquisition and processing of gene expression profile and clinical data

Expression profiles of mRNAs of 179 ESCC patients and matched clinical information were downloaded from the NCBI Gene Expression Omnibus (GEO) database (GSE53625) for differential expression and survival analysis. The clinical data, including eleven clinical variables, is shown in website: https://cdn.amegroups.cn/static/public/jtd-22-948-01.xlsx. The baseline characteristics of patients with ESCC in the GSE53625 was generated using the “tableone” R package, and the RNA-seq dataset was normalized via log2 conversion before further exploration.

Identification of differentially expressed PRGs in ESCC

The PRGs, including fifty-two genes, were obtained from previous studies and the Molecular Signatures Database (MSigDB database, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), and are listed in Table S1 (27,42-46). The PRGs-expression profiling data were extracted from the GSE53625 dataset, which is shown in website: https://cdn.amegroups.cn/static/public/jtd-22-948-02.xlsx. The “limma” and “reshape2” R packages were used to identify differentially expressed PRGs in ESCC tissues in comparison to the normal tissues, and the PPI network was established by STRING database (https://string-db.org/). The correlation network between the differentially expressed genes (DEGs) was built using the “igraph” R package and was visualized by the “reshape2” R package.

Functional enrichment analysis

GO analysis and KEGG analysis were used to explore the potential biological functions or pathways of DEGs. GO analysis of DEGs, including cellular component (CC) biological process (BP), and molecular function (MF) categories, was performed by the “clusterProfiler”, “org.Hs.eg.db” and “enrichplot” R packages. KEGG analysis of DEGs was conducted using the same R package, and the top 30 pathways or biological function were visualized using the “ggplot2” package.

GSSM analysis

According to a previous study, GSSM analysis by “GOSemSim” R package is a quantitative way to calculate the similarities between gene groups, including semantic similarity of GO terms (47). Therefore, the semantic similarity of GO terms between DEGs was measured using this package. The score of GO semantic similarity can be considered as functional similarity between genes.

Construction of the PRG prognostic model

The outcome of interest in this study was OS of ESCC. To construct an ESCC prognostic model related to PRGs, the GSE53625 dataset was split randomly into a training set and validation set via the “caret” R package, and the group allocation was performed in a 1:1 ratio (Table S2). Cox regression analysis was used to identify the prognostic PRGs in ESCC, and results were visualized by forest plots plotted via the “forestplot” R package. Survival analysis was performed via the “survival” R package, and Kaplan-Meier survival curves were established using the “survminer” R package. Log-rank tests were performed to obtain the P value of the Kaplan-Meier survival curves. LASSO-Cox regression analysis was then performed to construct a riskscore based on the prognostic PRGs, with the riskscore = (0.4714) × HMGB1 + (−0.1881) × IL-18 + (−0.0886) × NLRP7 + (0.2458) × PLCG1. The all patients with ESCC in the GSE53625 dataset were distributed into a high-risk group and low-risk group by the median of the riskscore, and the OS prediction ability of the riskscore was assessed via the Kaplan-Meier curves and the risk curves. When we identified riskscore as an independent prognostic factor, nine other clinical variables were included for analysis, including age, gender, tobacco use, alcohol use, tumor location, grade, tumor node metastasis (TNM) stage, anastomotic leak, pneumonia. To construct the prognostic model, all independent prognostic factors, including the riskscore, were used to construct the nomogram using a multivariable regression model, and for assessing the accuracy of the nomogram, the C-index, Receiver operating characteristic (ROC) curves, and Calibration curves were used. The model has no predictive power when the C-index is equal to 0.5, low accuracy when the C-index is 0.51–0.70, moderate accuracy when the C-index is 0.71–0.90, and high accuracy when the C-index is greater than 0.90 (48). The C-index of the nomogram and that of the TNM stage were compared to evaluate the accuracy of the nomogram further. All results of training set were verified by the validation set or all set.

Exploration of the role of the riskscore in ESCC immunotherapy

The TIDE algorithm (http://tide.dfci.harvard.edu), developed by Jiang et al., is based on integration and modelling of data from 189 human cancer research studies and can be used to evaluate the clinical response to immunotherapy (49). In the present study, the TIDE algorithm was utilized to assess the efficacy of immunotherapy in each ESCC patient. Subsequently, the Wilcoxon test was applied to evaluate the relationship between the riskscore and the efficacy of immunotherapy. A P value of less than 0.05 indicated the difference was significant.

Immune infiltration analysis

The activity or infiltration levels of 19 immunocytes or immune-related function was calculated by the single-sample gene set enrichment analysis (ssGSEA), which was run in the “GSEABase” and “GSVA” R packages (50). The correlation between the riskscore and the infiltration levels of the different immune cells or functions was analyzed via the Wilcoxon test. All results above were considered as significance if the P value were less than 0.05.

Statistical analysis

The R software (version 4.1.0) and various R packages were used for conducting the data analysis and visualizing the data. Statistical significance for all analyses was defined by P<0.05.

Results

Identification of the differentially expressed genes in ESCC

The transcriptome dataset of ESCC was obtained from the GEO databse (GSE53625), which consists of 179 ESCC samples and 179 normal samples. As shown in , the RNA-level expression of 52 PRGs in ESCC tissue and normal tissue were identified. We identified the differentially expressed PRGs via the differential expression analysis between ESCC samples and normal samples, and the results show there were 31 DEGs in ESCC (FDR <0.001) (). Of these, 16 DEGs were up-regulated in ESCC tissues, and the remainder were down-regulated.

Figure 2

Identification of differentially expressed PRGs in ESCC. (A) The mRNA expression of 52 PRGs in ESCC according to the whole genome expression profile data of ESCC from GEO database. [Red: tumor (N=179); Blue: normal (N=179)] (B) The logFC of 31 differentially expressed PRGs in ESCC. Red bars represent the up-regulated DEGs; blue bars represent the down-regulated DEGs. *, **, ***, suggest significance at P<0.05, P<0.01 and P<0.001, respectively. PRGs, pyroptosis-related genes; ESCC, esophageal squamous cell carcinoma; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes.

Functional enrichment analysis and the exploration of semantic similarity

To predict the CC, BP, and MF, GO enrichment analysis was performed, and the results showed the pyroptosis DEGs were enriched in positive regulation of cytokine production, interleukin-1 beta production, pyroptosis, inflammasome complex, inflammasome complex, late endosome, phospholipid binding, and cytokine receptor binding (). In addition, KEGG analysis suggested DEGs were primarily involved in the NOD−like receptor signaling pathway and necroptosis ().

Figure 3

Results of functional enrichment and GSSM analysis of DEGs. (A) GO analysis showing top-10 ranked BP, CC, and MF. (B) KEGG analysis showing top-30 ranked pathways. (C) GSSM analysis explored that the semantic similarity of GO terms of DEGs. The higher score of GO semantic similarity, the higher functional similarity. (D) The PPI network of DEGs. BP, biological process; CC, cellular component; MF, molecular function; GO, Gene Ontology; GSSM, GO-terms Semantic Similarity Measures; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.

To explore the semantic similarity of GO terms among the differentially expressed PRGs, the “GOSemSim” R package was used. The results suggested PYCARD, AIM2, CASP4, NOD1, and BAK1 had a higher semantic similarity of GO terms, which meant these genes had more interactions with other genes on the pathway and were more likely to be hub genes (). PPI network results of the DEGs are shown in , and the co-expression network of the DEGs is clearly documented in Figure S1.

Construction of a riskscore based on prognostic PRGs

The GSE53625 dataset was split randomly into a training set and validation set, and the characteristics of patients with ESCC in the two sets is shown in . To determine the PRGs with an impact on the prognosis of ESCC, Cox regression analysis was applied in the training set, and the result showed the expression levels of the four PRGs respectively had significant relations to the OS of ESCC patients (). Kaplan-Meier survival curves suggested the expression of HMGB1 (, P=0.013) and PLCG1 (, P=0.00069) were associated with poorer prognosis, and the expression of IL-18 ( P=0.045) and NLRP7 (, P=0.0045) were associated with a better prognosis. The riskscore based on the four prognostic PRGs mentioned above was constructed through the LASSO-Cox regression analysis (), with the riskscore = (0.4714) × HMGB1 + (−0.1881) × IL-18 + (−0.0886) × NLRP7 + (0.2458) × PLCG1. The training set was classified into a low-risk group and a high-risk group according to the median of the riskscore, and the survival analysis of the training set indicated the OS of the former was worse than that of the latter (, P<0.001). Similarly, the Kaplan-Meier survival curves of the validation set showed the low-risk group had a better OS (, P=0.03). In addition, patients from the training set were ordered (from left to right) according to the increasing riskscore, and the results showed the mortality risk increased and survival time decreased with the increase (). In the validation set, the trends of the distribution of survival states and survival time of patients were consistent with the trends in the training set, and the increase of deaths and decrease of survival time were significantly associated with the increased riskscore (). The results of principal component analysis (PCA) indicated ESCC patients could be significantly distinguished into high-risk and low-risk groups by the riskscore based upon the four prognostic PRGs ().

Table 1

Baseline characteristics of patients with ESCC in the GSE53625

Characteristics	All set (N=179)	Training set (N=90)	Validation set (N=89)
Status, n (%)
Alive	73 (40.8)	34 (37.8)	39 (43.8)
Death	106 (59.2)	56 (62.2)	50 (56.2)
Age (years), n (%)
≤65	130 (72.6)	69 (76.7)	61 (68.5)
>65	49 (27.4)	21 (23.3)	28 (31.5)
Gender, n (%)
Female	33 (18.4)	15 (16.7)	18 (20.2)
Male	146 (81.6)	75 (83.3)	71 (79.8)
Tumor location, n (%)
Upper	20 (11.2)	13 (14.4)	7 (7.9)
Middle	97 (54.2)	43 (47.8)	54 (60.7)
Lower	62 (34.6)	34 (37.8)	28 (31.5)
Grade, n (%)
I	32 (17.9)	16 (17.8)	16 (18.0)
II	98 (54.7)	52 (57.8)	46 (51.7)
III	49 (27.4)	22 (24.4)	27 (30.3)
TNM stage, n (%)
I	10 (5.6)	5 (5.6)	5 (5.6)
II	77 (43.0)	41 (45.6)	36 (40.4)
III	92 (51.4)	44 (48.9)	48 (53.9)
T, n (%)
T1	12 (6.7)	7 (7.8)	5 (5.6)
T2	27 (15.1)	14 (15.6)	13 (14.6)
T3	110 (61.5)	56 (62.2)	54 (60.7)
T4	30 (16.8)	13 (14.4)	17 (19.1)
N, n (%)
N0	83 (46.4)	42 (46.7)	41 (46.1)
N1	62 (34.6)	34 (37.8)	28 (31.5)
N2	22 (12.3)	8 (8.9)	14 (15.7)
N3	12 (6.7)	6 (6.7)	6 (6.7)
Anastomotic leak, n (%)
Yes	12 (6.7)	3 (3.3)	9 (10.1)
No	167 (93.3)	87 (96.7)	80 (89.9)
Pneumonia, n (%)
Yes	15 (8.4)	4 (4.4)	11 (12.4)
No	164 (91.6)	86 (95.6)	78 (87.6)
Tobacco use, n (%)
Yes	114 (63.7)	60 (66.7)	54 (60.7)
No	65 (36.3)	30 (33.3)	35 (39.3)
Alcohol use, n (%)
Yes	106 (59.2)	54 (60.0)	52 (58.4)
No	73 (40.8)	36 (40.0)	37 (41.6)

ESCC, esophageal squamous cell carcinoma; TNM, tumor node metastasis.

Figure 4

Identification of the prognostic PRGs. (A) Univariate Cox regression analysis was applied to determine the prognostic PRGs (P<0.05). Kaplan-Meier survival curves for the patients whose ESCC tissues expressed high and low level of HMGB1 (B), PLCG1 (C), IL-18 (D), and NLRP7 (E). CI, confidence interval; PRGs, pyroptosis-related genes; ESCC, esophageal squamous cell carcinoma.

Figure 5

Construction of the riskscore based on PRGs. (A) LASSO coefficient profiles of the four prognostic PRGs. (B) Ten-fold cross-validation of the riskscore. Kaplan-Meier survival curves of ESCC patients in the training set (C) and the validation set (D) showed the high-risk group (red) represented poor OS and low-risk group (blue) represented good OS. Risk score curve and scatter plots of the training set (E) and the validation set (F) indicated the number of death cases was increasing with the increase of the riskscore. The PCA of the training set (G) and the validation set (H). PRGs, pyroptosis-related genes; LASSO, Least Absolute Shrinkage and Selection Operator; ESCC, esophageal squamous cell carcinoma; OS, overall survival; PCA, principal component analysis.

Construction of a prognostic model related to the prognostic PRGs

We extracted total relevant clinical information from the GEO database (GSE53625). Subsequently, univariate () and multivariate () Cox regression analysis were applied in the training set, and the results showed the riskscore (P<0.05), TNM stage (P<0.05), and grade (P<0.05) were independent prognostic factors. At the same time, the riskscore, TNM stage, and grade in the training set were used to construct the nomogram for OS prediction (). The time-dependent ROC curves indicated the nomogram based on the three independent prognostic factors had excellent accuracy for the prognostic prediction of ESCC, as the Area under the ROC Curve (AUC) were 0.715, 0.820, and 0.728 for 1-, 3-, and 5-year (). With the aim of assessing the Calibration of the nomogram, Calibration curves predicting 1-, 3-, and 5-year OS were generated, and the results showed the actual OS rates were close to the nomogram-predicted OS rates (), while the C-index for the nomogram (0.705) was better than that for the TNM stage (0.620) (). In the all set () and the validation set (), the independent prognostic factors were also identified, and the results showed the riskscore (P<0.05), TNM stage (P<0.05), and age (P<0.05) were independent prognostic factors. The nomogram based on independent prognostic factors was also constructed in the all set and is shown in Figure S2A. The AUSs of the nomogram for predicting 1-, 3-, and 5-year OS in the all set were 0.701, 0.752, and 0.675, respectively (), and the Calibration curves for this nomogram also indicated the nomogram had a good prediction ability for OS rates (). In terms of the C-index, the nomogram (C-index =0.676) was also higher than the TNM stage (C-index =0.596) (). The nomogram of the validation set was also constructed (Figure S2B) and the time-dependent ROC curves analysis evaluated its accuracy. The nomogram of the validation set was verified by the ROC curves (), Calibration curves (), and C-index () and the results demonstrated the ability to predict the OS of ESCC patients was better than the TNM stage. In conclusion, the data above demonstrated the prognostic model related to PRGs had a satisfactory performance for predicting OS of ESCC.

Figure 6

Construction of a prognostic model related to prognostic PRGs in the training set. Univariate (A) and multivariate (B) Cox regression analyses of OS according to the riskscore and clinical factors. (C) The nomogram related to the riskscore for predicting the 1-, 3- , and 5-year OS of patients with ESCC. ROC curves (D) and the Calibration curves (E) for the nomogram for evaluating the predictive value. (F) Comparison between the C-index of the nomogram with that of the TNM stage suggested the prognostic model had a better OS prediction ability. CI, confidence interval; TNM, tumor node metastasis; AUC, area under the curve; PRGs, pyroptosis-related genes; OS, overall survival; ESCC, esophageal squamous cell carcinoma; ROC, receiver operating characteristic.

Figure 7

Validation set of the prognostic model related to the prognostic PRGs in all set and the validation set. Univariate and multivariate Cox regression analyses showed the riskscore was also an independent prognostic factor in all set (A,B) and the validation set (C,D). ROC curves and the Calibration curves for the nomogram of all set (E,F) and the validation group (H,I). The C-index of the nomogram of all set (G) and the validation set (J) were also higher than the TNM stage. CI, confidence interval; TNM, tumor node metastasis; AUC, area under the curve; PRGs, pyroptosis-related genes; ROC, receiver operating characteristic.

Role of the riskscore based on the PRGs in tumor immunotherapy

Pyroptosis is closely correlated with immunity, and the riskscore is based on the PRGs. Hence, we speculated the riskscore might be associated with tumor immunotherapy, and the relationships between the expression of immunotherapy-related mRNAs (PD-L1, CD23, CTLA4, IDO1, IFNG, IL-2, LAG3) and the riskscore were explored. The findings showed PD-L1 (, P=0.000046) and IL-2 (, P=0.0019) expression levels in the low-risk group were higher than in the high-risk group, and the high-risk group had a higher expression level of CD23 (, P=0.0023). In addition, the TIDE score was calculated using the TIDE web tool and patients were divided into two groups by its median. Patients with low TIDE scores showed better OS (, P=0.021). As described for , the high-risk group had a higher TIDE score, which meant that the effect of immunotherapy was significantly related with the riskscore, and the low-risk group had a better efficacy of immunotherapy (P=0.00034). In addition, patients in the high-risk group were more likely to experience immune evasion (Figure S3, P<0.001). Furthermore, the relationship between the riskscore and immune cells and function were analyzed using the ssGSEA R package, and the findings suggested significant differences were found in 10 kinds of immune cells, including Regulatory T Cells (Treg), tumor-infiltrating lymphocyte (TIL), T helper 2 cells (Th2 cells), T helper 1 cells (Th1 cells), T helper cells, Neutrophils, Neutrophils, Mast cells, Macrophages, Dendritic cells (DCs), and B lymphocytes (B cells), between the high-risk group and low-risk group (, P<0.05). Moreover, the scores of the immune functions in the low-risk group, including Antigen-presenting cell (APC)_co_inhibition, CC chemokine receptor (CCR), Parainflammation, and Type_II_interferons (IFN)_reponse, were higher than those in the high-risk group. Therefore, those immune functions played a greater role in low-risk group.

Figure 8

Exploration of the role of the riskscore in ESCC immunotherapy and the correlation between the riskscore and the immunotherapy-related mRNAs, including (A) PD-L1 (P=0.000046), (B) IL-2 (P=0.0019), and (C) CD23 (P=0.0023). (D) Kaplan-Meier survival curves for the patients with high and low level of TIDE score (P=0.021) (E) Relationship between the riskscore and TIDE score indicated the riskscore could predict the efficacy of immunotherapy (P=0.00034). (F) Differences in infiltration level of 19 immune cell types or the immune functions in the high-risk group and low-risk group. *, **, ***, stand for significance at P<0.05, P<0.01, and P<0.001, respectively, while ns, suggests no statistical significance. ESCC, esophageal squamous cell carcinoma; TIDE, Tumor Immune Dysfunction and Exclusion.

Discussion

Pyroptosis is a form of programmed cell death different from apoptosis, and is regulated by inflammasome and involved in tumorigenesis and progression of various malignancies (51,52). Accumulating evidence suggests pyroptosis has a dual effect on tumor growth and progression. It can promote tumor progression by releasing inflammatory cytokines and contributing to inflammation (53), and has the ability to inhibit tumorigenesis and tumor development caused by triggering tumor cell death (54). It has previously been suggested that a high expression level of Gasdermin B has a significant relation with the lower survival and increased metastasis in patients with breast cancer (55). NALP1 inflammasome has been proven to be a prognostic marker related to pyroptosis (56). However, to date there are no reports on prognostic signatures related to pyroptosis in ESCC, which was the objective of this study.

The first objective of this study was to determine the differentially expressed genes among PRGs in ESCC, and the results showed over half of all PRGs were DEGs in ESCC, indicating there was a relationship between ESCC and pyroptosis (, FDR <0.001). In accordance with our results, previous study have shown GSDME was up-regulated in ESCC tissues and a high expression level meant a better 5-year survival in patients with the disease (57). Additionally, it has been previously suggested that IL-18 expression of ESCC tissues was down-regulated and the deficiency of IL-18 enhanced the progression of ESCC in vivo and in vitro, which also was consistent with our findings (58,59). IL-18 was also one of the factors used to construct the riskscore related to pyroptosis.

The current study found DEGs most enriched in the positive regulation of cytokine production, interleukin-1 beta production, pyroptosis, inflammasome complex, late endosome, phospholipid binding, cytokine receptor binding, response to molecule of bacterial origin, and interleukin-8 production (). A previous study revealed IL-1 beta could promote the migration and invasion of ESCC cells, which was in accord with our results about functional enrichment (60). This study also found the five genes in the DEGs, including PYCARD, AIM2, CASP4, NOD1, and BAK1, had a higher semantic similarity of GO terms (). A literature review concluded PYCARD had dual role in cancers (61) by inhibiting tumors by mediating tumor cell death on the one hand, and promoting tumors by triggering the release of inflammatory cytokines and changing the tumor microenvironment on the other. This mechanism is similar to the dual effect mechanism on tumors related to pyroptosis. Moreover, PYCARD expression was usually silenced by methylation, and its methylation degree of PYCARD correlated with the depth of ESCC invasion (61). It has been previously confirmed that AIM2 expression was up-related in human cutaneous squamous cell carcinoma (cSCC) and its degree became higher with the grade of cSCC increasing (62). Similar histological features are seen in sSCC and ESCC, and previous work has shown this suggests a similar pathogenesis (63). Based on the results of these results and our data, we hypothesized AIM2 might play a role in ESCC.

In this study, we observed a significant relationship between the prognosis of ESCC and the four PRGs; PLCG1, HMGB1, IL-18, and NLRP7 (), and the riskscore was constructed based on these. Kaplan-Meier survival curves and the risk curves of the training set showed the riskscore could predict OS of ESCC patients with great accuracy (). The prognostic role of the riskscore was confirmed in the validation set, the related results of validation set being similar to those of the training set (). Prior study has revealed a signature based on extracellular matrix (ECM) which could also be used to predict the OS of patients with ESCC (64). Micro(mi)RNAs are non-coding RNAs and are significantly associated with malignant initiation and progression (65). A riskscore based on five prognostic miRNA was constructed in ESCC, and the high-risk group also had a worse survival (66). The present study is the first to establish a riskscore based on pyroptosis. We found it had a meaningful effect on the prognosis of ESCC patients and was an independent prognostic factor (, P<0.001). The riskscore and the other independent prognostic factors (TNM stage and grade) were used to construct the prognostic model related to pyroptosis in the training set. (), and the nomogram constructed in the training set could accurately and consistently predict the 1-, 3-, and 5-year OS (). Furthermore, the C-index of the prognostic model we constructed was significantly higher than that of the TNM stage (), and these results were validated in the validation set and all set. Our data also showed the riskscore was an independent prognostic factor in the validation set and all set (). The nomogram related to the riskscore obtained better OS prediction than the TNM stage, evidenced by the C-index of the nomograms being significantly higher than those of the TNM stage (). These findings clearly indicate our prognostic model related to PRGs had a better OS prediction ability because of the prognostic prediction value of the riskscore.

Pyroptosis is mediated by the gasdermin family of genes and is accompanied by immune and inflammatory responses (42). Caspase-1 has an important effect on innate immunity by triggering the release of the pro-inflammatory cytokines (IL-β and IL-18) and activating pyroptosis (18). Previous study has found a deficiency of IL-18 contributes to the impairment of anti-tumor immune response and the immune escape in ESCC (60). The results of functional enrichment also showed significant enrichment for the immune response (), suggesting the riskscore based on PRGs had some relation to immunotherapy in ESCC. We found the expressions of PD-L1 (), IL-2 (), and CD23 () were significantly associated with the riskscore. PD-1 is mainly expressed on activated T cells and is a significant immune checkpoint (67). It often combines with PD-L1, triggering tumor immune escape (68), and many studies have suggested that PD-L1 could be used to predict the efficacy for PD1/PD-L1 targeted therapy in a variety of tumor types (69). In 2019, an immune checkpoint inhibitor was approved by the FDA for use in the treatment of advanced ESCC patients with positive PD-L1 expression (70). In the present study, PD-L1 expression level in the low-risk group was higher than the high-risk group () and the high-risk group had a higher TIDE score, which meant immunotherapy showed better efficacy in the low-risk group (). These results indicate the riskscore based on PRGs plays a significant role in predicting the effect of immunotherapy.

In conclusion, we identified a riskscore based on four PRGs (PLCG1, HMGB1, IL-18, and NLRP7), which is an independent prognostic factor in ESCC. The prognostic model related to the riskscore and other independent factors was constructed, and the data demonstrated the prognostic model had a better OS prediction ability than the TNM stage. The riskscore could be taken as a reliable factor to guide the immunotherapy of ESCC. This study has several limitations. Firstly, the lack of an external cohort to validate the prediction accuracy of the prediction model is a major limitation of this study. Secondly, the number of cases included was only 179, which is small compared to other studies. Lastly, the present study is the only investigation into the prognostic model related to PRGs. Therefore, further studies should focus on in vitro and in vivo experiments to confirm our findings.

The article’s supplementary files as