The Pyroptosis-Related Signature Predicts Prognosis and Indicates Immune Microenvironment Infiltration in Gastric Cancer

Gastric cancer (GC) is one of the leading causes of cancer-related deaths and shows high levels of heterogeneity. The development of a specific prognostic model is important if we are to improve treatment strategies. Pyroptosis can arise in response to H. pylori, a primary carcinogen, and also in response to chemotherapy drugs. However, the prognostic evaluation of GC to pyroptosis is insufficient. Consensus clustering by pyroptosis-related regulators was used to classify 618 patients with GC from four GEO cohorts. Following Cox regression with differentially expressed genes, our prognosis model (PS-score) was built by LASSO-Cox analysis. The TCGA-STAD cohort was used as the validation set. ESTIMATE, CIBERSORTx, and EPIC were used to investigate the tumor microenvironment (TME). Immunotherapy cohorts by blocking PD1/PD-L1 were used to investigate the treatment response. The subtyping of GC based on pyroptosis-related regulators was able to classify patients according to different clinical traits and TME. The difference between the two subtypes identified in this study was used to develop a prognosis model which we named “PS-score.” The PS-score could predict the prognosis of patients with GC and his/her overall survival time. A low PS-score implies greater inflammatory cell infiltration and better response of immunotherapy by PD1/PD-L1 blockers. Our findings provide a foundation for future research targeting pyroptosis and its immune microenvironment to improve prognosis and responses to immunotherapy.


INTRODUCTION
Gastric cancer (GC) is the world's third-highest cause of death by cancer (Smyth et al., 2020). Each year, at least 1 million people worldwide are diagnosed with GC (Thrift and El-Serag, 2020). This disease is mostly detected in its advanced stages and abnormalities in the tumor microenvironment (TME) may lead to widespread tumor heterogeneity. Furthermore, there is significant heterogeneity with regards to the response of GC patients to therapy. So, the prognosis has not been improved (Bray et al., 2018). Pyroptosis refers to the cleavage of gasdermins via classical and non-classical pathways and can lead to the continuous expansion of cells until the cell membrane ruptures and causes the release of the cell contents, thus triggering a strong inflammatory response Ding et al., 2016;Rogers et al., 2017;Xia et al., 2019;Broz et al., 2020;Zhang Z. et al., 2020;Zhou et al., 2020). Pyroptosis plays an important role in antagonizing infection and endogenous danger signals. Pathogens such as H. pylori or chemotherapy drugs can cause pyroptosis in patients with GC. Pyroptosis creates a tumor-suppressive environment by releasing inflammatory factors, however, it can also weaken the body's immune effect on tumor cells and accelerate tumor growth in different cancers (Zaki et al., 2010;Ellis et al., 2011;Chen et al., 2012;Ma et al., 2016). However, the effect of pyroptosis on the prognosis of GC is not clear.
The classification of GC patients by next-generation sequencing is a novel method that can quickly identify cancer characteristics and inform us about the most appropriate treatment strategies. Drug treatment already uses HER2 as a predictive biomarker (Smyth et al., 2020). But the value of HER2 in the prognosis of GC remains controversial (Tanner et al., 2005;Park et al., 2006;Gravalos and Jimeno, 2008;Ruschoff et al., 2010;Begnami et al., 2011;Yoon et al., 2012). PD-L1, as an immunotherapy index, also requires further verification (Shitara et al., 2018). Other biomarkers are currently being evaluated. Due to the lack of subgroup classifications, clinical practice cannot be guided by molecular subtypes. Therefore, there is an urgent need for the development of an effective gene signature to indicate prognosis and to guide clinical treatment, especially with regards to targeted therapy and immunotherapy.
In the present study, we aimed to build a scoring model (that produced the PS-score) by classifying GC patients based on pyroptosis-related regulators to predict prognosis and guide clinical treatment. We clustered 618 patients with GC according to pyroptosis-related genes and identified two types of pyroptosis-related subtypes that were related to prognosis and immune infiltration. On this basis, the PSscore can be determined by constructing a pyroptosis-related model using the LASSO-Cox method. This score is able to predict prognosis, immune infiltration, and immunotherapy response. Our findings indicate the potential connection between pyroptosis, prognosis, the immune microenvironment, and the response to immunotherapy of GC patients.

Sources of Gastric Cancer Datasets and Preprocessing
The workflow chart (Supplementary Figure 1) describes which samples were utilized at each stage of statistical analysis. Microarray data from Affymetrix R was obtained from Gene Expression Omnibus (GEO). Batch effects from non-biological technical biases were corrected by using the "ComBat" algorithm of the "SVA" package. All of the clinical information used in this study are publicly available in the GEO database. As to datasets in TCGA, RNA sequencing data (FPKM value) of gene expression were downloaded from UCSC. Patients without survival information were removed from further analysis.

Human Clinical Specimens
Twenty-two pairs of RNA samples of GC and adjacent nontumor tissues were obtained from Jinan Central Hospital, Shandong, P. R. China. The study protocol was approved by Shandong University Research Ethics Committee.

Consensus Clustering
Consensus clustering was applied to identify distinct pyroptosisrelated patterns relating to the expression of pyroptosis regulators by the k-means method. The number of clusters, and their stability, were determined by the consensus clustering algorithm using the "ConsensuClusterPlus" package (John Hartigan, 1979). We performed 1,000 times repetitions to guarantee the stability of our classification (Wilkerson and Hayes, 2010).

Differentially Expressed Genes (DEGs)
We used the empirical Bayesian approach of the "limma" package to obtain DEGs (Ritchie et al., 2015). The significance criteria for selecting DEGs was set as an adjusted P < 0·05 and an absolute value of Log2 FC ≥ 0·8.

TME Cell Infiltration
We used the CIBERSORTx algorithm and EPIC to quantify the proportions of immune cells. For CIBERSORTx, we uploaded the normalized gene expression data to the web portal using LM22 signature and 1,000 permutations (Newman et al., 2015). EPIC is a web-based analytical and discovery platform for analyzing mass cytometry data from immune cells in a standardized manner (Yeo et al., 2020). Tumor purity scores were estimated by the "ESTIMATE" package (Yoshihara et al., 2013;Yang et al., 2019).

Multivariate Cox Regression Analyses
We used univariate Cox regression and multivariate Cox regression analyses for overall survival (OS) in four GEO datasets described earlier. A false discovery rate (FDR) < 0.05 was used as a statistical boundary. The results of multivariate prognostic analysis for pyroptosis-related subgroups were acquired by application of the "forestplot" package.

The Establishment of a PS-Score Scoring Model and Prognostic Analysis
We established an efficient prediction model using LASSO−Cox analysis. OS was then used to derive the most useful predictive features from the training cohort (Lossos et al., 2004). PS-score = k i=1 βiPi where k, βi, Pi represented the number of signature genes, the coefficient index, and the gene expression level, respectively. The cut-off point was determined using the "survminer" package. We used Kaplan-Meier survival curves to identify the ability of the model to distinguish different subtypes of patients and time-dependent receiver operating characteristic curves (ROC curves) to determine the efficiency of the model. The C-index was calculated by the "survcomp" package and compared using the "cindex.comp" package.

RNA Extraction and Real-Time PCR
Total RNA was extracted by Trizol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's protocol. Subsequently, the extracted RNA was reversetranscribed using PrimeScript RT reagent Kit with gDNA Eraser (Takara, Japan). The cDNAs were subjected to SYBR Green-based real-time PCR analysis. The primers used in real-time PCR assays were listed in Supplementary Table 1.

Statistical Analysis and Cut-Off Value
Correlation coefficients were computed by Spearman's and distance correlation analyses. Log-rank tests were utilized to identify the significance of differences in survival curves. The cut-off value mentioned in this article was 1·258 as the best cut-off value from the "survminer" package. ROC curves, time-dependent ROC curves, and the area under curves (AUC) were derived using the "pROC" and the "timeROC" packages, respectively. Comparisons of the integrated area under the curves (IAUC) were carried out with the "iauc.comp" package. The "RCircos" package allowed us to plot the copy number variation landscape of pyroptosis regulators in 23 pairs of chromosomes (Mayakonda et al., 2018). Analyses between the two groups were performed using the Wilcox test. The Kruskal Wallis test was also used to compare three or more groups. Gene expression data, and all statistical analyses were carried out in R 4.0.0, GraphPad Prism 8, and SPSS26 software. Clinical features were compared by Chisquare tests or Fisher's exact tests where appropriate. All statistical P values are two-side and P < 0·05 represents statistical significance.

Overview of Genetic Changes and Expression Variations of Pyroptosis-Related Regulators in GC
We analyzed the network of potential biological functions associated with the 11 pyroptosis-related regulators by using the STRING platform ( Figure 1A). The regulators focused predominantly on the regulation of immune response and pyroptosis. At the genetic level, 59 of the 433 samples (about 13.63%) showed pyroptosis-related regulator mutations. Of these, CASP5 showed the highest frequency of mutations. We did not identify any GSDME mutations in any of the GC samples ( Figure 1B). We also found CNVs in 7 of the 11 pyroptosisrelated regulators in GSE62717; these were common changes and most were concentrated on copy number amplification ( Figure 1C). We identified the alterations of the 7 regulators featuring CNVs on the chromosome (Figure 1D). At the expression level, these 11 regulators were able to help us distinguish normal samples from tumor samples in GC patients ( Figure 1E). Compared with normal samples, except for GSDME and GZMA, the remaining regulators all showed increases in GC samples ( Figure 1F and Supplementary Figure 2A). To more effectively verify our findings, we tested the mRNA levels of 11 pyroptosis-related genes in 22 pairs of tumor tissues and normal adjacent tissue samples gathered from our hospital. The same finding was that these regulators were upregulated in tumors ( Figure 1G). In order to better explore the possible relationship between the genetic level and the expression level of pyroptosis-related regulators, we combined the CNV pattern of TCGA-STAD and its expression level for joint analysis. Interestingly, we found that the expression levels of CASP1, CASP3, CASP4, CASP5, and CASP8, all increased with increases in copy number; a similar relationship was identified for GSDMD (Supplementary Figures 2B-G). We hypothesized that changes in CNV may be an important factor that can lead to abnormal gene expression. Our analysis showed that the expression levels of pyroptosis-related regulators were related to GC, thus suggesting they may reflect different traits in patients. The expressions of pyroptosis-related regulators between normal tissues (n = 100) and gastric tissues (n = 300) in GSE66229 cohort (Wilcox test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001; ns, not statistically significant). (G) The mRNA levels of pyroptosis-related regulators in twenty-two pairs of GC and their paired adjacent normal tissues (n = 22) were measured by real-time PCR (Paired t-test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001).
Frontiers in Cell and Developmental Biology | www.frontiersin.org

Identification of a GC Classification Pattern Mediated by 11 Pyroptosis-Related Regulators
We created a queue using four GEO datasets from the same platform along with OS data and clinical information. Based on the expression levels of 11 pyroptosis-related regulators, we identified two different regulation patterns by using the unsupervised clustering method, including 267 cases in pyroptosis-related cluster 1 and 351 cases in pyroptosisrelated cluster 2 ( The survival advantage of cluster 1 was higher than that of cluster 2 ( Figure 2B). To explore the differences in biological behavior between these two patterns, we performed GSVA enrichment analysis ( Figure 2C and Supplementary Table 3). Cluster 1 showed enrichment in terms of pathways associated with immune activation, including antigen processing and presentation, TOLL, NOD, and RIG I-like receptor signaling pathways, B cells and activated T-cell receptor signaling pathways, the activated JAK-STAT signaling pathway, base excision repair, and H. pylori infection resistance. Pyroptosis is known to be converted from apoptosis with different bacterial or viral infections (Rogers et al., 2017;Zhang Z. et al., 2020;Zhou et al., 2020). Therefore, cluster 1 was also significantly enriched in the apoptotic pathway. Cluster 2 showed enrichment in carcinogenic activation pathways, such as the mTOR, TGFβ, NOTCH, and WNT signaling pathways. Subsequently, we also confirmed that the two regulatory patterns could be distinguished by the application of the 11 pyroptosisrelated regulators ( Figure 2D).

Differences in TME Infiltration and Clinical Characteristics Between Two Pyroptosis-Related Subtypes
Next, we analyzed cell infiltration data and found that activated innate immune cell infiltration was abundant in cluster 1, including the presence of dendritic cells, M1 macrophages, B cells, along with activated CD 4 and CD 8 T cells, thus conferring a significant survival advantage. Cluster 2 was enriched with endothelial cells, mast cells, M2 macrophages, and resting T4 memory cells (Figures 3A,B). Although some tumor tissues possessed a large number of immune cells, these immune cells could not penetrate the tumor but were forced to stay in the surrounding matrix. The activation of the matrix in a tumor microenvironment was therefore considered to be immunosuppressive (Chen and Mellman, 2017). In the present analysis, the ESTIMATE score also showed that the immune score of cluster 1 was higher than that of cluster 2 (Supplementary Figure 3J). We found that these two regulatory patterns had completely different TME cell infiltration characteristics. Cluster 1 was classified as an immune-inflamed phenotype while cluster 2 was classified as an immune-excluded phenotype (Chen and Mellman, 2017). We also found that the proportion of TME cell types in the different clusters was different while the composition types were the same (Supplementary Figure 3K). These suggested that although these regulatory patterns did not increase or decrease the types of immune infiltrating cells, it was likely to change their proportion and thus change the TME.
To further explore the clinical manifestations of these two different clusters, we focused on the ACRG cohort which represents the most comprehensive study of clinical information relating to our 618 patients ( Figure 3C). Most regulators were expressed at high levels in cluster 1. Patients with EMT molecular subtypes (Cristescu et al., 2015) were prominent in cluster 2, while patients with MSI subtypes were classified into cluster 1. Patient survival status also corresponded well to the survival advantages described in our previous research. Compared with other pathological types of GC, patients with signet ring cell carcinoma were mainly associated with the cluster 2 pattern. Moreover, we noticed that patients with early GC were associated with a cluster 1 regulatory pattern and that patients with advanced GC were mainly associated with the cluster 2 regulatory pattern. This explains why the cluster 1 regulatory pattern was associated with a better survival advantage. We conducted statistical analyses on 282 patients that possessed recurrence data in the ACRG cohort. The cluster 2 pattern was associated with a greater number of recurrence cases ( Figure 3D). These findings showed that different pyroptosis-related patterns represented different GC features and had different TME statuses.

Development and Validation of a Gene Signature Based on Pyroptosis-Related Clusters
To better apply these subtypes to the clinical treatment of GC and determine a specific score for every patient, we next explored differences between the two patterns and determined a specific gene signature. We also quantified the gene signature so that it could be applied to the diagnosis and treatment of each patient. First, we identified 113 DEGs with an absolute value of Log2 FC < 0·8 and p < 0·05 associated with the two different regulatory patterns ( Figure 4A and Supplementary  Table 4). Next, univariate and multivariate Cox regression analysis identified 22 genes that could be used as an independent prognostic signature ( Figure 4B and Supplementary Tables 5, 6). Ten genes, identified as the most immune response-related genes, appeared in cluster 1 while 12 genes encoding cancer occurrence proteins tended to be more prevalent with the cluster 2 pattern (Figure 4C). To build a model that would be able to quantify each patient, six of the 22 DEGs were retained by application of LASSO-Cox regression model with a minimum of λ. We used these to build a pyroptosis-related signature score which we named the "PS-score" (Figure 4D and Table 1). Next, we attempted to further determine the value of PS-score by predicting the prognosis of patients. We divided patients into high and low PS-score groups with the best cut-off value of 1·258. We found that the low group had a clear survival advantage over the high group (Supplementary Figures 4A,B).
To prove the universal indicative value of the PS-score, we also verified this score in more cohorts and obtained the same results (Figures 4E,F). We further proved that the PS-score was a good indicator for the 3-year survival and 5-year survival of GC patients (Figure 4G and Supplementary Figure 4C). Because The consensus score matrix of all samples when k = 2 in GEO cohorts (GSE15459, GSE34942, GSE57303, and GSE62254). Two samples were more likely to be grouped into the same cluster when there was a higher consensus score between them in different iterations. (B) OS curves for the two pyroptosis-related clusters based on 618 patients with gastric cancer from four GEO cohorts (GSE15459, GSE34942, GSE57303, and GSE62254) (Log-rank test, p < 0.0001). OS, Overall survival. (C) The heatmap was used to visualize biological processes analyzed by GSVA which showed the active biological pathways in distinct pyroptosis-related clusters (Bayes moderation, * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001). (D) Principal component analysis for the expression of pyroptosis-related regulators to distinguish cluster 1 (n = 267) from cluster 2 (n = 351) in GEO cohorts (GSE15459, GSE34942, GSE57303, and GSE62254).
the prognostic labels of GC have been discussed extensively in recent years and play an important role in early diagnosis and treatment, we also compared several other rigorous prognostic models to evaluate the important function of the PS-score for evaluating prognosis. The analysis showed that the PS-score was better than other models for predicting the prognosis of GC FIGURE 3 | Different pyroptosis-related clusters showed diverse clinical features and TME cell infiltration. (A,B) The abundance of every type of TME infiltrating cells between the two pyroptosis-related clusters analyzed, respectively, by CIBERSORTx and EPIC in GEO cohorts (GSE15459, GSE34942, GSE57303, and GSE62254) (Wilcox test, * P < 0.05; * * P < 0.01; * * * P < 0·001; * * * * P < 0.0001; ns, not statistically significant). (C) Consensus clustering of differential expression genes between the two pyroptosis-related clusters in GSE62254 cohort. Columns of the heatmap represented 300 gastric cancer samples (Chi-square test or Fisher's exact test, * * P < 0.01; * * * P < 0.001). (D) Unsupervised clustering of differential expression genes between the two pyroptosis-related clusters in GSE62254 cohort. Columns of the heatmap represented 282 gastric cancer patients with recurrence records (Chi-square test, * * * P < 0.001).
patients (Cho et al., 2011;Cai et al., 2020Wang et al., 2020: Figure 4H). Therefore, we included PS-score as an effective indicator and other clinical characteristics, into Cox regression analysis and found that PS-score and stage were both two factors that independently affected the prognosis of GC patients ( Table 2). Interestingly, we found that the combination of In the LASSO-Cox model of GSE62254, the minimum standard was adopted to obtain the value of the super parameter λ by 10-fold cross-validation. The λ value was confirmed as 0.07558 where the optimal lambda resulted in 6 non-zero coefficients. (E) OS curves for the different PS-score subgroups with the cut-off value 1,258 about 300 patients with gastric cancer from GSE62254 cohort (Log-rank test, p < 0.0001). (F) OS curves for the different PS-score subgroups with the cut-off value 1·258 among 618 gastric cancer samples from four GEO cohorts (GSE15459, GSE34942, GSE57303, and GSE62254) (Log-rank test, p < 0.0001). (G) The time-dependent receiver operating characteristic (ROC) analysis of the PS-score. The area under the curve (AUC) was 0.727, 0.738 at 3 years, and 5 years, respectively, in GSE62254 cohort. (H) ROC curves about PS-score ISSGC Score, Risk Score, and GPSGC Score in GSE62254 cohort (Mann Whitney tests; compared with PS-score, * * * * P < 0.0001). PS-score and GC stage could better predict the survival of patients (Supplementary Figure 4D). These findings indicated that the PS-score was a promising potential evaluation indicator for the prognosis of GC patients.

Low PS-Scores Identified the Alleviation of Clinical Characteristics and Immune Pathway Activation
We also found that the high-score group in the TCGA-STAD cohort predicted poor survival characteristics with the same cutoff value ( Figure 5A). And we were pleasantly surprised to find that PS-score not only predicted the overall survival of patients with GC but also could predict the disease-specific survival of patients specifically (Supplementary Figure 4E). Subsequently, to verify the universal applicability of PS-score for all digestive tract cancers, we collected data from seven digestive tract cancer samples including gastric cancer and divided these into two groups with the same cut-off value. We were surprised to find that the low-score group also showed a survival advantage ( Figure 5B). This implied that the PS-score may have an important value as a prognostic indicator not only in GC but also in other forms of gastrointestinal cancer.
To determine the specificity of the PS-score in patients with different clinical manifestations, we analyzed the relationship between clinical traits and PS-score in the ACRG cohort. We found that the cluster 1 pattern, with an obvious survival advantage, presented with obviously low scores; this was consistent with our previous research ( Figure 5C). Patients with high scores also experienced more relapses (Figure 5D). We were surprised to find that with an increase in patient survival stage, the PS-score showed a gradually increasing trend, thus showing that the score for patients with advanced gastric cancer was higher than that for patients with early gastric cancer ( Figure 5E). This suggested that PS-score had the potential as a clinical index to quantify the survival risk of GC patients. Similarly, the EMT subtype with a poor response to treatment (Cristescu et al., 2015) also exhibited high scores (Figure 5F). We also found that different Lauren classifications and pathological types showed completely different PS-score values (Supplementary  Figures 4F,G). Compared with other pathological types, patients with signet ring cell carcinoma which were mainly associated with the cluster 2 pattern, had high PS-scores. We also found that second/third-line treatment was the commonly used strategy in the clinic. However, this system lacks a suitable measurement standard. A low PS-score indicated a good response to followup treatment, thus indicating that PS-score may represent a good marker for GC treatment ( Figure 5G). As the main cause of gastric cancer, the screening of H. pylori has significantly reduced the prevalence of GC (Smyth et al., 2020). Therefore, we specifically analyzed the relationship between H. pylori infection and PS-score and found that patients with positive infection may have higher scores (Supplementary Figure 4H). Although the limitation of the sample size made the statistical significance less obvious, we still observed a trend consistent with the previous results. This gave us a hint, but more verification is still needed.
In addition to different clinical phenotypes, the expression of different pyroptosis-related regulators in patients would also affect the specificity of the PS-score. We found that regardless of whether we used the ACRG cohort or the TCGA-STAD cohort, except for GSDME, the expression levels of other regulators were significantly reduced in the high-score group (Figure 5H and Supplementary Figure 4I). GSDME was generally expressed at higher levels in the high-score group; this may be related to the cytokine release syndrome (CRS) induced by patients with high GSMDE expression during treatment (Liu et al., 2020a). CRS was characterized by fever, hypotension, and respiratory insufficiency associated with elevated serum cytokines (Davila et al., 2014), which was positively correlated with GSDME expression level (Liu et al., 2020a). Therefore, treatments for GSDME should be used with caution to avoid undesirable side effects (Ibrahim et al., 2021). Although GSDME is considered to be a probable tumor suppressor gene (Wang et al., 2017), its precise function needs to be explored further. We were unable to judge the prognosis of patients by considering only the level of GSDME expression. Similarly, it would be unscientific to predict the survival ability of GC patients solely by considering the expression of pyroptosisrelated regulators. Our research aimed to create a quantitative pyroptosis-related model to comprehensively predict prognosis and the value of treatment. Based on the specificity data supporting the PS-score, we explored the pathway enrichment associated with the PS-score to identify the internal mechanisms involved. We found that the low-score group was significantly enriched in inflammatory signaling pathways. The PS-score showed a positive response to base damage repair and response to H. pylori infection and was also related to ubiquitination modification and drug metabolism regulation. However, the high-score group was mainly enriched in signal pathways related to cancer development ( Figure 5I and Supplementary Table 7). These findings showed that PS-score exerted specificity in different patients associated with different patterns of signaling pathway activation and could be used to evaluate certain clinical characteristics and therapeutic effects in GC patients.

The PS-Score Could Predict Prognosis in Clinical Scenarios and Represent TME Differences
Given the importance of the PS-score in predicting the prognosis of GC patients, we next attempted to explore its value for clinical application. We constructed a nomograph featuring seven clinical features that were easily accessible and generally believed to have a certain impact on the prognosis of GC and the ability of the PS-score to predict the survival rates of GC patients at 3, 5, and 8 years ( Figure 6A). A C-index of 0·764 indicated that the nomogram had a good predictive value ( Figure 6B). When we used the PS-score as a separate indicator to distinguish patients with GC, we found the 300 patients of the ACRG cohort were divided into a high group and a low group (Figure 6C). An alluvial diagram was used to visualize the changes in patient characteristics ( Figure 6D).
Based on these results, we found that PS-score could play an important role in clinical prediction. Next, we investigated whether PS-score would have a guiding value for clinical treatment, especially immunotherapy. We analyzed the infiltration of TME cells with different PS-scores (Figures 6E-G). Immune activation-related cells, such as activated T cells, NK cells, M1 macrophages, dendritic cells, neutrophils, showed significant negative correlations with PS-score. The higher scores were closely related to resting memory cells, monocytes, mast cells, cancer-associated fibroblasts, and endothelial cells. The immunescore from ESTIMATE analysis decreased as the PSscore increased, while the stromalscore showed the opposite effect. These data showed that the low-score group had a stronger immune response than the high-score group. The differences in TME cells may be the main reason for the heterogeneity of PS-score. Tumors that could attract more T cell infiltration are referred to as "hot tumors" and are more sensitive to immunotherapy and showed better immunotherapy effects . Clinical trials and preclinical studies have also revealed that patients with higher somatic tumor mutational burden (TMB), when treated with immune checkpoint blockade therapy, were associated with enhanced responses, long-term survival, and lasting clinical benefits (Zhang B. et al., 2020).  (C-G) PS-score in different clinical trait constituents including pyroptosis-related cluster, recurrence status, TNM staging, ACRG subtype in GSE62254 cohort, and followup treatment response in TCGA-STAD (Wilcox test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001; ns, not statistically significant). (H) Differential expression of pyroptosis-related regulators in low PS-score subgroup (n = 209) and high PS-score subgroup (n = 91) of GSE62254 cohort (Wilcox test, * * * * P < 0.0001). (I) Visualization of biological processes analyzed by GSVA in distinct PS-score subgroups (Bayes moderation, * P < 0.05; * * P < 0.01; * * * * P < 0.001; * * * * P < 0.0001).
Fortunately, we explored the TMB of different PS-scores and found that the lower group had a higher TMB, which suggested a better immunotherapy response ( Figure 6H). These results proved that PS-score may have an application value for predicting the prognosis of GC patients, and would reflect the response of immunotherapy to a certain extent. The correlation between every type of TME infiltrating cells and PS-score analyzed, respectively, by CIBERSORTx, ESTIMATE score, and EPIC in cohorts (GSE15459, GSE34942, GSE57303, GSE62254, and TCGA-STAD; Consolidation was the cohort consolidating four GEO cohorts) (Spearman test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001). Sizes of circles represented relevant correlation coefficients. (H) Tumor burden (TMB) of low PS-score subgroup (n = 215) and high PS-score subgroup (n = 153) in TCGA-STAD (Wilcox test, * * * * P < 0.0001).

The Role of PS-Score in Anti-PD1/PD-L1 Immunotherapy
Previous studies have shown that PS-score may suggest the effect of immunotherapy. There is also evidence that patients with a high TMB status show long-lasting clinical responses to anti-PD1/PD-L1 immunotherapy (Zhang B. et al., 2020). Recent studies have shown that there was a clear connection between the expression of PD-L1 and pyroptosis. Therefore, we speculated that there may be a connection between PS-score and immunotherapy (Hou et al., 2020). We first checked the expression changes of immune checkpoints. We compared the differences in the expression levels of immune checkpoint genes in the ACRG cohort and the TCGA-STAD cohort ( Figure 7A and Supplementary Figure 5A). We found that the low-score group showed higher expression levels of immune checkpoint genes, thus suggesting a better response to immunotherapy. Tumor Immune Dysfunction and Exclusion (TIDE) is a computational framework to model two primary mechanisms of tumor immune evasion which can provide predicted results about immunotherapy (Jiang P. et al., 2018;Chen et al., 2021). High TIDE may predict patients with suppressive cells inhibiting T cell infiltration as non-responders. To better illustrate the predictive power of the PS-score for immunotherapy, we applicated TIDE in ACRG cohort. We were pleasantly surprised to find a positive correlation between TIDE and PS-score ( Figure 7B). Furthermore, predicted responses suggested that PS-score may be a good predictor of immunotherapy in GC ( Figure 7C). As we all know, "hot tumor" is more sensitive to immunotherapy (Galon and Bruni, 2019). We had also verified the higher PS-score represented worse prognostics in "hot tumor" such as breast cancer and kidney cancer (Figures 7D,E). Due to the lack of published data on anti-PD1/PD-L1 immunotherapy in GC, we investigated published three datasets: IMvigor210, NCT01358721, and GSE78220 (Choueiri et al., 2016;Hugo et al., 2016;Mariathasan et al., 2018). Interestingly, when we analyzed the prognosis of patients with these three cancers, we found that PS-score could be a good predictor for melanoma patients and metastatic renal cell carcinoma patients with 1year survival, and metastatic urothelial cancer patients' 2-year survival rate ( Figure 7F and Supplementary Figures 5B,C). So we guessed that PS-score had a connection with immunotherapy in other cancers. Then we detected changes in the expression of immune checkpoints presented by different PS-scores and also found that the low-score group had higher expression ( Figure 7G), thus suggesting that PS-score may be an important indicative index that plays the same important role in other cancers. As we hypothesized, we found that patients who responded to immunotherapy also showed a lower PS-score ( Figure 7H and Supplementary Figures 5D,E). In metastatic urothelial cancer, different immunological subtypes may lead to a completely different therapeutic response. We also found that they represented different levels of PS-score ( Figure 7I). These results explained the potential value of PS-score in immunotherapy, and to a certain extent proved that PS-score may be used as a predictor of immunotherapy responses.

DISCUSSION
Chronic infection of the gastric mucosa leads to the gradual development of atrophic gastritis and intestinal metaplasia, thus promoting the progression of GC (Smyth et al., 2020). Molecular signatures associated with distinct clinical outcomes have been delineated in various solid tumors to improve clinical management through the development of personalized medicine (Sorlie et al., 2001;Cho et al., 2011;Higgins and Baselga, 2011;Li et al., 2013;Roepman et al., 2014). Therefore, we need to explore the changes in the status and mechanisms of GC cells associated with the immune environment to facilitate treatment. Pyroptosis occurs in cells infected by pathogens, as an embodiment of programmed cell death, thus inducing the body's inflammatory response (Bedoui et al., 2020). Under the stimulation of pathogens, apoptosis can thus be converted into pyroptosis. Pyroptosis plays various roles in many cancers. It has the effect of inhibiting tumor growth in colorectal cancer, liver cancer, and skin cancer (Zaki et al., 2010;Ellis et al., 2011;Ma et al., 2016), but a two-way effect in breast cancer (Chen et al., 2012). So we cannot judge the prognostic value of GC based on the expression of several gasdermins alone. Therefore, we explored all the pathways directly related to pyroptosis and explored a prognostic signature by analyzing the influence of the involved pathways on the tumor microenvironment. Gasdermins blocker is under development, but there is insufficient evidence to support it. Our signature provides potential targets for targeted therapy of pyroptosis, especially CASP1 and GZMB. At present, pyroptosis has been considered for use in anti-tumor therapy, and our research suggests that pyroptosis combined with immunotherapy to improve the prognosis of patients may be an effective treatment direction.
The classification of samples based on predefined gene expression characteristics is a proven method (Cristescu et al., 2015). Our subtyping strategy drew on this method and classified GC patients based on the expression of pyroptosis-related regulators. We showed that the expression of these regulators was completely different when compared between the two clusters due to various heterogeneities. These regulators were also significantly associated with different survival risks. Our analysis culminated in several consensuses: (1) most pyroptosis-related regulators showed high expression levels in cluster 1; (2) cluster 2 was a separate subtype which showed a worse prognosis; (3) cluster 1 was determined to be an immune-inflamed phenotype and cluster 2 was determined as an immune-excluded phenotype; and (4) collectively considering clinical information and RNA data was more likely to reflect the cellular phenotypes.
Clinical trials have tested anti-tumor molecular targeted drugs in all GC types regardless of the molecular subtypes involved. For example, the expression of immune checkpoint molecules differs across different subtypes, therefore immunotherapy should be distinguished. To create a better clinical application value, we developed a scoring model (PS-score) to quantify the prognostic risk based on the two clusters. Our study provided strong evidence for the clinical management of GC. First, the PS-score takes into account the heterogeneity of patients. Second, this score can link pyroptosis and prognosis. Specifically, the PS-score featured both tumor suppressor genes and tumor-promoting genes and allocated these with different weightings. The PS-score included but was not limited to pyroptosis-related regulators such as GZMB and CASP1. Third, the PS-score represented patients with different clinical traits and was related to immunotherapy. A high PS-score showed worse clinical traits and a lower predicted survival time. TME cell infiltration data demonstrated that the PS-score holds important value for immunotherapy. FIGURE 7 | A powerful role of the PS-score scoring model in PD1/PD-L1 immunotherapy. (A) Differential expression of immune checkpoint genes in low PS-score subgroup (n = 209) and high PS-score subgroup (n = 91) of GSE62254 cohort (Wilcox test, * * * * P < 0·0001). (B) The relationship between TIDE and PS-score in GSE62254 cohort (Spearman test, p < 0·0001). (C) Different PS-score in responder group (n = 112) and non-responder group (n = 188) in GSE62254 cohort (Wilcox test, * * * * P < 0.00001). (D,E) OS curves for the PS-score with the cut-off value 1,258 of samples in TCGA-ACC cohort (Log-rank test, p = 0.0017) and TCGA-BRCA cohort (Log-rank test, p = 0.023). (F) ROC curves about PS-score in IMvigor210. (G) Differential expression of immune checkpoint genes in low PS-score subgroup (n = 126) and high PS-score subgroup (n = 222) of IMvigor210 cohort (Wilcox test, * * P < 0.01; * * * * P < 0.0001). (H). Different PS-score in CR/PR group (n = 68) and SD/PD group (n = 230) in IMvigor210 cohort (Wilcox test, * P < 0.05). SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. (I) Different PS-score in IMvigor210 cohort's immune phenotypes (Kruskal-Wallis test, P = 1.8e-06).
More activated immune cell infiltration in patients with a low PSscore led to a better response to immunotherapy. Compared with TMB and PD-L1 expression data, the PS-score is affordable and provides more informative outcomes. Fourth, the PS-score may also apply to other gastrointestinal tumors and immune-related tumors. Finally, compared with other models, the PS-score is directly focused on the death mode of GC cells. Researchers have investigated prognostic models of GC in hypoxia or under modified conditions of m6A already or have considered the immunoscore (Jiang Y. et al., 2018;Liu et al., 2020b;Zhang B. et al., 2020). Our study was more focused on the factors that directly caused tumor cell death and changed the tumor microenvironment. In this manner, our model is more valuable for facilitating treatment.
Viral and bacterial infections in the stomach can trigger downstream signaling pathways through pattern recognition receptors (PRRs) such as NLRs and TLRs. The activated caspases cleave the pyroptosis-related genes, causing cell dilation and death, and releasing mature inflammatory factors, especially IL-18 and IL-1. This process not only responds to infection but also alters the tumor microenvironment. More immune cell infiltration means more favorable immunotherapy. At the same time, the occurrence of pyroptosis, such as the increase of GSDMD splice, can inhibit MAPK, mTOR, and Wnt signaling pathways to various degrees. We speculate that it plays an anti-tumor role and improves prognosis mainly by inhibiting cell proliferation. In the process of chemotherapy, GSDME cleavage can promote the transformation of cells from apoptosis to pyroptosis and promotes TIL function, thus creating conditions for immunotherapy. GSDMB cleaved by GZMA also can convert apoptosis into pyroptosis, and IFN-γ promoted GZMA-or NK cell-induced pyroptosis in several target cells. High-level expression of GSDMB in cancer cells enhanced tumor clearance in a mouse model. These also mean that pyroptosis-related genes may be able to predict the prognosis and prompt immunotherapy.
Our study aimed was to classify patients with GC into subtypes, identify DEGs and build a prognostic model, and link pyroptosis with patient prognosis. Although we had performed multi-angle and multi-database verifications, this study still had limitations that need to be considered. The model created did not have a good predictive value in terms of all-time survival stages in patients undergoing immunotherapy. This may have been caused by the specificity of different cancers and requires more extensive research. Tumor heterogeneity is indeed a problem that cannot be ignored, and more targeted improvements for different types of tumors may be proposed with the development of pyroptosisrelated researches. There is also a need for excavating or building more gastric cancer immunotherapy data and H. pylori infection data. The results of some single-cell sequencing should be able to explain the specific changes in the tumor microenvironment, which is also an aspect of our attention in the future. Moreover, our model should be validated further by performing both in vitro and in vivo experiments to better evaluate the relationship between the PS-score and pyroptosis of cells after infection. These have not only increased the challenges but also added hope for us to make us more motivated to continue digging.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Committee of the School of Basic Medical Science, Cheeloo College of Medicine, Shandong University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
JJ designed the study, obtained the funding, and supervised the study. WS, ZY, LZ, and FL organized the data. WS and ZY performed the analysis. YF checked the statistical method. JJ, WS, ZY, and YF prepared the figures. WS, ZY, YF, and LC wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The National Natural Science Foundation of China (Nos. 81971901, 81772151, and 82002107) and the Department of Science and Technology of Shandong Province (No. 2018CXGC1208) were used for purchasing codes and training.