Comprehensive characterization of pyroptosis phenotypes with distinct tumor immune profiles in gastric cancer to aid immunotherapy

Results

Genetic and transcriptomic alterations of pyroptosis-related genes across gastric tumors

Here, we observed the functions of 33 pyroptosis-relevant genes in gastric tumors. We firstly identified their chromosomal locations (Figure 1A). The mRNA expression levels of above genes were investigated in 32 controls and 375 gastric tumors from TCGA-STAD cohort. The results showed that most of pyroptosis-related genes had higher expression levels in tumors (Figure 1B). The incidence of somatic mutations of pyroptosis-related genes in gastric cancer was then summarized. Among 433 specimens, 117 occurred somatic mutations of pyroptosis-associated signatures, with frequency of 27.02% (Figure 1C). We investigated that PLCG1, CASP5, CASP8 and NLRP3 displayed the highest mutation frequencies in gastric cancer specimens. Further analysis exhibited the widespread CNV frequencies in 33 pyroptosis-related genes and most were focused on the amplification in CNV (Figure 1D). These findings indicated that such genetic variations affected their instability in gastric tumors. The PPI network showed the extensive interplay of such pyroptosis-relevant genes (Figure 1E), indicative of the potential implication of their cross-talk in gastric cancer pathogenesis and progression.

Figure 1. Comprehensive characterization of genetic and expression alterations of pyroptosis-related genes in gastric cancer. (A) Circos 2D track plots showing the location of 33 pyroptosis-related genes on chromosomes. (B) The mRNA expression of 33 pyroptosis-relevant genes in 32 normal and 375 gastric cancer tissue specimens from TCGA-STAD cohort. The asterisks indicated the statistical p-values (^*p < 0.05; ^**p < 0.01; ^***p < 0.001) and ns indicated no statistical significance. (C) The somatic mutation frequency of 33 pyroptosis-relevant genes in 433 gastric cancer samples in TCGA-STAD dataset. The upper barplot displayed somatic mutation. The stacked barplot exhibited the fractions of conversions of specimens. The number on the right showed the mutation frequency in each pyroptosis-related gene. The right barplot showed the proportion of each mutation type. (D) Summary of the CNV frequency of 33 pyroptosis-related genes. The column height indicated the mutation frequency. Amplification, red; deletion, blue. (E) The PPI network visualizing the interactions between 33 pyroptosis-relevant genes.

Identification of three distinct pyroptosis phenotypes in gastric cancer

Four gastric cancer cohorts (GSE15459, GSE62254, GSE84437 and TCGA-STAD) were integrated and batch effects were removed (Figure 2A, 2B). Finally, 1275 gastric cancer patients were integrated as the meta-cohort. Univariate-cox regression analysis of 33 pyroptosis-related genes showed that GPX4 (HR: 1.24 (1.06–1.46); p = 0.008), NLRP7 (HR: 0.86 (0.76–0.96; p = 0.009), CASP6 (HR: 0.87: (0.76–0.99); p = 0.032), NLRP6 (HR: 0.90 (0.82–0.99); p = 0.04), IL6 (HR: 1.07 (1.00–1.13); p = 0.047), CASP5 (HR: 0.90 (0.82–0.99); p = 0.023), CASP4 (HR: 0.84 (0.71–0.999); p = 0.048) and TIRAP (HR: 0.73 (0.59–0.91); p = 0.004) had significant correlations to prognosis across 1275 gastric cancer patients. Through the ConsensusClusterPlus package, we stratified gastric cancer specimens to three different pyroptosis phenotypes according to the expressions of above 8 prognostic pyroptosis-related genes with unsupervised clustering method (Figure 2C–2E; Supplementary Table 1). We named these phenotypes as pyroptosis cluster A (n = 625), B (n = 352) and C (n = 298). There were remarkable differences in the transcriptional profiles of the prognostic pyroptosis-relevant gene set among three diverse pyroptosis phenotypes (Figure 2F). Prognostic analysis for the three pyroptosis phenotypes showed the particularly prominent survival advantage in pyroptosis cluster-C, while pyroptosis cluster-B displayed the poorest survival outcomes (Figure 2G). The accuracy of this pyroptosis-related classification was externally proven in the GSE26253 (Supplementary Figure 1A–1D). The differences in prognosis among three pyroptosis phenotypes were also confirmed in the GSE26253 cohort (Supplementary Figure 1E).

Figure 2. Identification of three distinct pyroptosis phenotypes in gastric cancer. (A, B) Integration of four gastric cancer cohorts (GSE15459, GSE62254, GSE84437 and TCGA-STAD) into one meta-cohort and removal of batch effects. (C) Heat map showing the optimal clustering by k = 3 based on prognostic pyroptosis-related genes in the meta-cohort. (D) The consensus cumulative distribution function (CDF) for the cumulative distributive function of the consensus matrix for k = 2–9 identified by different colors. (E) The CDF plots of consensus clustering matrix for k = 2–9. (F) The t-SNE plots of the mRNA expression profiles of prognostic pyroptosis-related genes showing the three distinct subtypes indicated by different colors. (G) Kaplan–Meier survival curves for the three pyroptosis subtypes based on the 1275 gastric cancer specimens in the meta-cohort (log-rank test). (H) Heatmap showing the GSVA score of the 50 hallmark pathways in three pyroptosis phenotypes.

Pyroptosis phenotypes with different immune profiles

For uncovering potential biomolecular alterations in three different pyroptosis phenotypes, this study carried out GSVA on these 50 hallmark gene sets. We investigated that pyroptosis cluster-C owned distinct enrichment by stromal pathways (e.g., TGF-β signaling, epithelial mesenchymal transition and angiogenesis; Figure 2H). Pyroptosis cluster-B was highly enriched in carcinogenic activation pathways such as PI3K-AKT-mTOR signaling, Wnt-β-catenin signaling, glycolysis and KRAS signaling. Moreover, pyroptosis cluster-A markedly displayed the activation of immune-relevant processes like allograft rejection, interferon gamma response, inflammatory response, complement, IL2-STAT5 signaling and IL6-JAK-STAT3 signaling. These data were indicative of the survival outcomes of the three pyroptosis phenotypes. Moreover, ssGSEA was adopted for estimating the abundant levels of 28 immune cell populations within the TME among distinct pyroptosis phenotypes. Immune-excluded phenotype is rich in immune cells that are reserved in the matrix ambient the tumor cell nest instead of penetrating its parenchyma [20]. Pyroptosis cluster-C could be categorized as immune-excluded tumors with innate immune cell infiltrations as well as stromal activation (Figure 3A). The cluster-B could be categorized as immune-desert tumors with low infiltration levels of nearly all immune cells, indicative of immunity suppression. Moreover, the cluster-A could be categorized as immune-inflamed tumors with rich adaptive immune cells as well as immune activation. The immune infiltration traits across three pyroptosis phenotypes were confirmed via CIBERSORT algorithm (Supplementary Figure 2). Figure 3B showed the prominent heterogeneity in immune checkpoints among pyroptosis phenotypes. Pyroptosis cluster-C exhibited the highest expression of LAG3, CTLA4, TNFRSF9, ICOS, CD80, TNFSF9, ICOSLG, KIR3DL1, PDCD-1, TNFRSF8, TNFRSF15, TNFRSF14, HHLA2, CD244, CD27, BTLA, LGALS9, TMIGD2, CD28, TNFRSF25, CD40LG, CD160, and CD200R1. Through ESTIMATE algorithm, we quantified the entire infiltrations of immune/stromal cells or tumor purity among the pyroptosis phenotypes. Our results showed that pyroptosis cluster-B had the highest stromal score, followed by pyroptosis cluster-A and pyroptosis cluster-C (Figure 3C). Inversely, pyroptosis cluster-B displayed the lower tumor purity, indicating that tumors in the cluster-B were potentially neighbored by rich non-cancer cells such as stromal/immune cells (Figure 3D). In Figure 3E, pyroptosis cluster-B had the highest immune score, followed by pyroptosis cluster-A and pyroptosis cluster-C. Collectively, our data demonstrated that the different pyroptosis phenotypes had the features of distinct immune landscapes.

Figure 3. Pyroptosis phenotypes with different immune landscapes in the gastric cancer meta-cohort. (A) The relative abundance of 28 immune cells in the TME among distinct pyroptosis phenotypes through ssGSEA. (B) The mRNA expression of the main immune checkpoints among three pyroptosis phenotypes. The asterisks indicated the statistical p-values (^*p < 0.05; ^**p < 0.01; ^***p < 0.001). (C–E) The stromal score, tumor purity and immune score in three pyroptosis phenotypes.

Identification of pyroptosis phenotype-related gene signatures

The pyroptosis phenotype-related transcriptional expression alterations across three pyroptosis phenotypes were detected in the gastric cancer meta-cohort. Through limma package, we identified 1629 overlapped pyroptosis phenotype-relevant DEGs among the three phenotypes (Figure 4A; Supplementary Table 2). These DEGs represented the key distinguishing indicators for three pyroptosis phenotypes. GO enrichment analyses of above DEGs suggested that immunity-relevant mechanisms presented the marked enrichment (Figure 4B, 4C). This indicated that such DEGs were pyroptosis phenotype-relevant genes. According to the 1629 pyroptosis phenotype-relevant gene signatures, this study established three reliable genomic phenotypes (Figure 4D–4G). Gastric cancer patients were clustered into three different pyroptosis genomic phenotypes, named as pyroptosis gene cluster-A (n = 511), cluster-B (n = 510) and cluster-C (n = 254). Moreover, Figure 4H shows the prominent heterogeneity in pyroptosis phenotype-associated genes and clinicopathological characteristics among three pyroptosis genomic phenotypes.

Figure 4. Identification of pyroptosis phenotype-related gene signatures in the gastric cancer meta-cohort. (A) Venn diagram showing the overlapped pyroptosis phenotype-related DEGs among three pyroptosis phenotypes. (B, C) GO and KEGG enrichment results of the overlapped pyroptosis phenotype-related DEGs. (D) Heat map showing the optimal clustering by k = 3 based on the overlapped pyroptosis phenotype-related DEGs. (E) The CDF for the cumulative distributive function of the consensus matrix for k = 2–9 identified by different colors. (F) The CDF plots of consensus clustering matrix for k = 2–9. (G) The t-SNE plots of the mRNA expression of pyroptosis phenotype-related gene signatures visualizing three distinct pyroptosis genomic phenotypes. (H) Heatmap showing the differences in the expression of pyroptosis phenotype-related gene signatures and clinicopathological characteristics among three pyroptosis genomic phenotypes.

Pyroptosis genomic phenotypes with distinct prognosis and immune landscapes

Survival analysis showed the prominent prognostic differences among three pyroptosis genomic phenotypes in the gastric cancer meta-cohort. Pyroptosis gene cluster-C owned the poorest prognosis, while pyroptosis gene cluster-B exhibited a significant survival advantage (Figure 5A). Consistent with pyroptosis phenotypes, pyroptosis genomic cluster-C could be classified into immune-excluded tumors with innate immune cell infiltrations as well as stromal cells, with the genomic cluster-B categorized into immune-desert tumors with immunity suppression, and the genomic cluster-A categorized into immune-inflamed tumors with rich adaptive immune cells as well as immune activation (Figure 5B–5E). The heterogeneity in transcriptomic levels of immune checkpoints was also investigated across three pyroptosis genomic phenotypes (Figure 5F).

Figure 5. Pyroptosis genomic phenotypes with distinct prognosis and immune landscapes in the gastric cancer meta-cohort. (A) Kaplan-Meier survival curves for three pyroptosis genomic phenotypes (log-rank test). (B) The relative abundance of 28 immune cells among three pyroptosis genomic phenotypes through ssGSEA. (C–E) The immune score, stromal score, and tumor purity across three genomic phenotypes. (F) The mRNA expression of immune checkpoints across three genomic phenotypes. The asterisks indicated the statistical p-values (^*p < 0.05; ^***p < 0.001) and ns indicated not significant.

Generation of the pyroptosis score and evaluation of its clinical implication

To accurately quantify pyroptosis phenotypes in individual gastric cancer, this study proposed a scoring system named the pyroptosis score according to the pyroptosis phenotype-relevant gene signatures with PCA algorithm. Alluvial diagram showed the heterogeneity of the quantification of pyroptosis score among different pyroptosis phenotypes, pyroptosis genomic phenotypes and survival status (Figure 6A). Notably, pyroptosis cluster-B owned the strongest pyroptosis score, with subsequent cluster-A or cluster-C (Figure 6B). Moreover, pyroptosis genomic cluster-C displayed the strongest pyroptosis score, with subsequent cluster-A or cluster-B (Figure 6C). Next, gastric tumors in the meta-cohort were stratified as low or high pyroptosis score group. The efficiency of the scoring system in inferring survival outcomes was then observed. In Figure 6D, high pyroptosis score displayed poorer survival. The pyroptosis score distribution across different clinical characteristics was also analyzed. Patients with advanced grade and stage had higher pyroptosis score (Figure 6E, 6F). Further analysis of time-independent ROC curves confirmed the well performance of the scoring system (Figure 6G). Thus, the scoring system could be utilized for evaluating certain clinicopathological characteristics (grade and stage) of gastric cancer patients.

Figure 6. Generation of the pyroptosis score and evaluation of its clinical implication in the gastric cancer meta-cohort. (A) Alluvial diagram of pyroptosis phenotypes with different pyroptosis genomic phenotypes, pyroptosis score and survival status. (B) Distribution of pyroptosis score in three pyroptosis phenotypes. (C) Distribution of pyroptosis score in three pyroptosis genomic phenotypes. (D) Kaplan-Meier survival curves for high and low pyroptosis score groups (log-rank test). (E) Distribution of pyroptosis score in different grades (G1-3). (F) Distribution of pyroptosis score in different stages (stage I–IV). (G) Time-independent ROC curves showing the predictive performance of the pyroptosis score, age, gender, grade, and stage in gastric cancer prognosis.

Characteristics of the pyroptosis score in molecular subtypes and tumor genome somatic mutations

GSEA results uncovered that high pyroptosis score was significantly connected to oncogenic or immune-related pathways (Figure 7A), whereas low pyroptosis score exhibited significant associations with genetic mutations (Figure 7B). Significantly mutated genes were compared in two groups. As shown in mutational landscape, higher somatic mutational frequency was found in low pyroptosis score group (33.49% vs. 36.26%; Figure 7C, 7D). TTN (17% vs. 20%), TP53 (15% vs. 20%), MUC16 (11% vs. 12%), etc. exhibited increased genetic mutation frequencies for patients with low pyroptosis score. Somatic copy number alterations (SCNA) exhibit positive correlation to immune evasion as well as proliferative capacity of tumor cells. We observed an increased SCNA for patients with low pyroptosis score (Figure 7E). The high microsatellite instability (MSI-H) subtype, characterized by favorable survival outcomes exhibited distinct correlation to lower pyroptosis score, while low MSI (MSI-L) and microsatellite stable (MSS) status displayed higher pyroptosis score (Figure 7F, 7G).

Figure 7. Characteristics of the pyroptosis score in molecular subtypes and tumor genome somatic mutations. (A, B) GSEA showing the activated signaling pathways in samples with high and low pyroptosis scores in the gastric cancer meta-cohort. (C, D) Mutational landscape in TCGA-STAD cohort clustered by high and low pyroptosis score subgroups. (E) Relative distribution of SCNA in high and low pyroptosis score groups from TCGA-STAD cohort. (F) Relative distribution of pyroptosis score in different molecular subtypes (MSI-H, MSI-L and MSS) in TCGA-STAD cohort. (G) Relative distribution of MSI status in patients with high and low pyroptosis scores from TCGA-STAD cohort.

The pyroptosis score predicts immunotherapeutic responses and chemosensitivity

We evaluated the interaction of the pyroptosis score with several signaling pathways by Spearman analyses. Consequently, the pyroptosis score presented negative correlations to DNA damage repair processes but was positively correlated to stromal-related signatures in the meta-cohort (Figure 8A). Immunotherapies like CTLA-4/PD-1 blockage have strictly made progress for anti-cancer treatment. Tumor mutation burden (TMB) has become an effective biomarker for predicting immunotherapeutic response. Low pyroptosis score presented notable connection to increased TMB in TCGA-STAD cohort (Figure 8B). IPS is also recommended as an estimator of immunotherapeutic efficacy. The analysis revealed that IPS score of CTLA-4/PD-1 inhibitor treatment was distinctly elevated in low pyroptosis score in TCGA-STAD cohort (Figure 8C–8F). The response to common chemotherapeutic agents was estimated in each sample from the meta-cohort (Figure 8G–8J). We observed that high pyroptosis score exhibited markedly reduced IC50 values of cisplatin and docetaxel as well as low pyroptosis score had significantly increased IC50 value of paclitaxel. Thus, patients with high pyroptosis score could be more sensitive to cisplatin and docetaxel as well as low pyroptosis score was more likely to benefit from paclitaxel. Due to the close association between the pyroptosis score and immune response, this study further observed the predictive performance of the pyroptosis score in immunotherapeutic responses in two anti-PD-1 cohorts. The prominent therapeutic benefits and immune responses to PD-1 inhibitor therapy were confirmed in low pyroptosis score group (79%) than another group (21%) in the GSE78220 cohort (Figure 9A). Patients with low pyroptosis score displayed prominent survival advantages (Figure 9B). The consistent results also confirmed that low pyroptosis score was distinctly correlated to clinical benefits and prolonged survival time in another anti-PD-1 therapy (Liu et al. study; Figure 9C, 9D). Collectively, above data prominently demonstrated that the pyroptosis score could estimate immunotherapeutic efficiency.

Figure 8. The roles of the pyroptosis score in predicting immunotherapeutic responses and chemosensitivity. (A) Interactions of the pyroptosis score with several signaling pathways through Spearman analyses. Red, positive correlation and blue, negative correlation. (B) Distribution of TMB score in high and low pyroptosis score groups in TCGA-STAD cohort. (C–F) Distribution of IPS score of CTLA-4/PD-1 in high and low pyroptosis score groups in TCGA-STAD cohort. (G–J) The estimated IC50- values of cisplatin, docetaxel, gemcitabine, and paclitaxel in patients with high and low pyroptosis scores in the meta-cohort.

Figure 9. Validation of the predictive efficacy of the pyroptosis score for immunotherapeutic responses in two anti-PD-1 cohorts. (A) The fraction of patients with clinical response to anti-PD-1 therapy in high and low pyroptosis score groups from the GSE78220 cohort. (B) Kaplan-Meier curves for patients with high and low pyroptosis scores in the GSE78220 cohort (log-rank test). (C) The fractions of patients who presented therapeutic responses to anti-PD-1 therapy in high and low pyroptosis score groups from the Liu cohort. (D) Kaplan-Meier curves of patients with high and low pyroptosis scores in the Liu cohort (log-rank test). Abbreviations: CR: complete response; PR: partial response; SD: stable disease; PD: progressive disease.

Validating prognostic pyroptosis-relevant genes in gastric tumors

Through western blot, we validated the transcriptomic levels of such prognostic pyroptosis-relevant genes in three paired gastric cancer tissues and normal tissues (Figure 10A). Our results confirmed that the expression of GPX4 (Figure 10B), CASP6 (Figure 10C), IL-6 (Figure 10D), CASP5 (Figure 10E), CASP4 (Figure 10F), and TIRAP (Figure 10G) was greatly up-regulated in gastric tumors relative to normal tissues. Meanwhile, the expression of NLRP7 (Figure 10H) and NLRP6 (Figure 10I) was greatly attenuated in gastric cancer tissues than controls. Above evidence confirmed the abnormal expression of the prognostic pyroptosis-related genes in gastric cancer.

Figure 10. Validation of the expression of prognostic pyroptosis-related genes. (A–I) Western blot for the expression of (B) GPX4, (C) CASP6, (D) IL-6, (E) CASP5, (F) CASP4, (G) TIRAP, (H) NLRP7 and (I) NLRP6 in three paired gastric cancer tissues and normal tissues. The asterisks indicated the statistical p-values (^*p < 0.05; ^**p < 0.01; ^***p < 0.001).

Materials and Methods

Sample curation and preprocessing

Transcriptome data or complete patient information of gastric cancer were collected from the Cancer Genome Atlas (TCGA) or Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). In total, five gastric cancer cohorts (GSE15459, GSE62254, GSE84437, GSE26253 and TCGA-stomach adenocarcinoma (STAD)) were ultimately included for further analysis. RNA-seq profiles (FPKM value) of 32 controls and 375 gastric tumors were retrieved from TCGA-STAD through the Genomic Data Commons (https://portal.gdc.cancer.gov/) with TCGAbiolinks package, with subsequent conversion of FPKM to TPM value [29]. The original “CEL” files of microarrays on the Affymetrix platform were obtained, with subsequent background adjustment and quantile standardization. Furthermore, the standardized matrix of microarrays on the Illumina platform was acquired. Utilizing “ComBat” function from sva tool, log2 (TPM + 1) was constructed [30], and the negative values that represented abnormal values were removed, and thus batch effects of integrated GSE15459, GSE62254, GSE84437 and TCGA-STAD datasets were removed. Somatic mutations or copy number variations (CNVs) were curated from TCGA. The maftools package was employed for analyzing and visualizing the somatic mutation data. Figure 11 shows the workflow diagram of this study.

Figure 11. The workflow diagram in this study.

Unsupervised clustering for pyroptosis-related genes

Thirty-three pyroptosis-related genes were curated from the published studies. Their chromosomal locations were drawn via Rcircos tool [31]. Protein-protein interactions (PPIs) of 33 pyroptosis-related genes were conducted with the STRING (https://string-db.org/) [32]. Univariate cox regression analysis was applied to screen prognostic pyroptosis-related genes. By unsupervised clustering analysis, different pyroptosis phenotypes based upon prognostic pyroptosis-associated gene set were established and gastric cancer cases were classified for subsequent exploration. The amount and reliability of clustering were evaluated with consensus clustering algorithm. Through ConsensuClusterPlus tool, such steps as well as 1000 times repetitions were carried out for insuring the clustering reliability [33]. The number of clusters k was set as 2~9. Moreover, 80% of the samples were sampled utilizing a re-sampling approach. T-distributed stochastic neighbor embedding (t-SNE) was run for proving the subtype assignment utilizing above ferroptosis genes.

Gene set variation analysis (GSVA)

The 50 hallmark gene sets were curated from the Molecular Signatures Database (MSigDB) [34]. GSVA enrichment analysis was then presented to estimate the activity of these pathways and biological processes in gastric cancer samples with GSVA package [35].

Estimation of immune cell infiltration

Through single sample gene set enrichment analysis (ssGSEA), the relative abundance of 28 immune cells was quantified in the TME. Immune cell markers were curated from previous research [36, 37]. The relative abundance of each immune cell was normalized to ssGSEA score, ranging from 0 to 1. CIBERSORT, a deconvolution approach based upon linear support vector regression, was also applied for estimating the abundance of 22 immune cell populations [38].

Quantifying immune response predictors

The mRNA expression levels of the main immune checkpoints (IDO1, LAG3, CTLA4, TNFRSF9, ICOS, CD80, PDCD1LG2, TIGIT, CD70, TNFSF9, ICOSLG, KIR3DL1, CD86, PDCD1, LAIR1, TNFRSF8, TNFSF15, TNFRSF14, IDO2, CD276, CD40, TNFRSF4, TNFSF14, HHLA2, CD244, CD274, HAVCR2, CD27, BTLA, LGALS9, TMIGD2, CD28, CD48, TNFRSF25, CD40LG, ADORA2A, VTCN1, CD160, CD44, TNFSF18, TNFRSF18, BTNL2, C10orf54, CD200R1, TNFSF4, CD200 and NRP1) were quantified in each gastric cancer sample. The ESTIMATE algorithm was employed for calculating immune/stromal score, thus inferring immune/stromal cell infiltration levels as well as tumor purity based upon the transcriptional profiling [39]. Tumor mutation burden (TMB) score was calculated as the amount of mutations/length of exons (30 Mb). Immunophenoscore (IPS), an indicator in predicting response of CTLA-4 or anti-PD-1 blockade, quantifies determining factors of cancer immunogenic characteristics as well as characterizes the immune profiles within the tumor and cancer antigenome [37]. IPS was calculated based on the MHC-associated signatures, checkpoint or immunomodulatory molecules, effector or suppressor cell populations.

Identification of pyroptosis phenotype-related differentially expressed genes (DEGs)

Through limma package, DEGs (adjusted p < 0.05) between distinct pyroptosis phenotypes were determined [40]. Then, overlapped pyroptosis phenotype-related DEGs were extracted for further analysis.

Construction of the pyroptosis score

The consensus clustering algorithm was adopted to define genomic clustering number and reliability according to the expression of overlapped pyroptosis phenotype-related DEGs. Prognostic analysis of each pyroptosis phenotype-related DEG was carried out utilizing univariate cox regression models. The DEGs with p < 0.05 were chosen for subsequent exploration. Principal component analysis (PCA) was adopted for constructing the pyroptosis score. Principal component (PC) 1 and 2 were extracted to define the pyroptosis score. Such method possessed the advantages of focusing the score on the set with the most block of well correlated (or anticorrelated) genes in the set, and down-weighting contributions from genes that don’t track with other set members. The pyroptosis score was calculated utilizing an approach like GGI [41, 42], as follows: pyroptosis score = ∑(PC1i+PC2i), of which i indicated prognostic pyroptosis phenotype-associated DEG expression.

Functional enrichment analyses

Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment on pyroptosis phenotype-related DEGs was run with clusterProfiler package [43]. Gene set enrichment analysis (GSEA) was presented to observe the biological pathways involved in pyroptosis score based upon the “c2.cp.kegg.v6.2.symbols” set.

Relationship of pyroptosis score with well-established biological signatures

The gene sets of several well-established biological signatures were collected according to the published literature [44–46]. Spearman correlation was carried out for assessing interactions between pyroptosis score and above biological pathways.

Assessment of the sensitivity of chemotherapeutics

Cisplatin, docetaxel, gemcitabine, and paclitaxel were selected as candidate chemotherapeutic agents. To predict the sensitivity of these chemotherapeutic agents, IC50 of each sample was measured utilizing pRRophetic tool [47].

Transcriptome and follow-up data of immunotherapy cohorts

Two immunotherapy cohorts were collected for our analysis: metastatic melanoma patients receiving PD-1 blockade treatment (GSE78220 cohort) [48] and Liu et al. [49]. Gene expression profiles were curated and converted to TPM values for further analysis.

Patients and specimens

Three gastric cancer cases who had surgically proven primary gastric cancer and received D2 radical gastrectomy were recruited from Xiamen Haicang Hospital. Tumors and adjacent normal tissue specimens of three cases were collected during surgery, which were immediately placed in liquid nitrogen and stored at −80°C. This study gained the approval of the Medical Ethics Committee of Xiamen Haicang Hospital (KY-2020014). All patients provided written informed consent.

Western blot

Radioimmunoprecipitation assay (RIPA; Beyotime, China) buffer was applied for protein extraction from tissues. The supernatants were run on 8%–12% acrylamide gels via SDS-PAGE and then transferred onto polyvinylidene difluoride membranes (Millipore, USA). Antibodies against GPX4 (#P36969; RayBiotech, USA; 1:1000), NLRP7 (#AB117732; Abcam, USA; 1:1000), CASP6 (#AB108335; Abcam; 1:1000), NLRP6 (#P59044; RayBiotech; 1:1000), IL-6 (#AB271042; Abcam; 1:2000), CASP5 (#AB40887; Abcam; 1:3000), CASP4 (#AB238124; Abcam; 1:1500), TIRAP (#P58753; RayBiotech; 1:1000), and GAPDH (#60004-1-Ig; Proteintech, China; 1:20000) and HRP-conjugated secondary antibodies (#ab97080 or ab47827; Abcam; 1:2000) were employed for western blot. The chemiluminescence western blot detection system (Bio-Rad, USA) was utilized for protein detection.

Statistical analysis

Statistical analysis was generated by R 3.6.1. Kaplan-Meier survival analysis was run through Survminer and survival packages. Statistical significance between groups was estimated through student’s t or Wilcoxon test. Meanwhile, three groups were compared with Kruskal-Wallis test or one-way analysis of variance. The receiver operating characteristic (ROC) curve was conducted to investigate the predictive efficacy of pyroptosis score and common clinical features (age, gender, grade and stage) in gastric cancer prognosis, with subsequent estimation of area under the curve (AUC) utilizing timeROC package. P values were two-side, with p < 0.05 as statistical difference.

Availability of data and materials

The data used to support the findings of this study are included within the supplementary information files.

Abbreviations

STAD: stomach adenocarcinoma; TME: tumor microenvironment; TCGA: the Cancer Genome Atlas; GEO: Gene Expression Omnibus; CNV: copy number variation; PPI: protein-protein interaction; STRING: Search Tool for the Retrieval of Interacting Genes; t-SNE: t-distributed stochastic neighbor embedding; GSVA: Gene set variation analysis; MSigDB: Molecular Signatures Database; ssGSEA: single sample gene set enrichment analysis; ESTIMATE: Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data; TMB: tumor mutation burden; IPS: Immunophenoscore; DEGs: differentially expressed genes; PCA: principal component analysis; PC: principal component; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GSEA: Gene set enrichment analysis; pan-F-TBRS: pan-fibroblast TGF-β response signature; EMT: epithelial-mesenchymal transition; IC50: half-maximal inhibitory concentration; ROC: receiver operating characteristic; AUC: area under the curve.

Author Contributions

Shuitu Feng, Fanghong Luo, Jianmin Zhuang and Wenhui Zheng conceived and designed the study. Kaida Huang, Yubiao Lin, Guoqin Qiu and Shengyu Wang conducted most of the experiments and data analysis, and wrote the manuscript. Lihua Feng, Zhigao Zheng, Yingqin Gao and Xin Fan participated in collecting data and helped to draft the manuscript. All authors reviewed and approved the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Ethical Statement and Consent

The study was approved by the Ethics Committee of Xiamen Haicang Hospital (KY-2020014). All patients provided written informed consent.

Funding

This research was sponsored by Natural Science Foundation of Xiamen, China (3502Z20227423, 3502Z20227426, 3502Z20227307); Fujian Natural Science Foundation Project (No. 2023D007); The Fujian provincial health technology project (2022QNB031); The science and technology program of Xiamen city (3502Z20209252 and 3502Z20209164); The science and technology planning project of Haicang District Bureau of industry and information technology (350205Z20202006, 350205Z20212001, 350205z2021 2003, 350205Z20222004, 350205Z20222005).

Editorial Note

This corresponding author has a verified history of publications using a personal email address for correspondence.

References

• 1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
• 2. Qiu XT, Song YC, Liu J, Wang ZM, Niu X, He J. Identification of an immune-related gene-based signature to predict prognosis of patients with gastric cancer. World J Gastrointest Oncol. 2020; 12:857–76. https://doi.org/10.4251/wjgo.v12.i8.857 [PubMed]
• 3. Wang X, Cheng G, Miao Y, Qiu F, Bai L, Gao Z, Huang Y, Dong L, Niu X, Wang X, Li Y, Tang H, Xu Y, Song X. Piezo type mechanosensitive ion channel component 1 facilitates gastric cancer omentum metastasis. J Cell Mol Med. 2021; 25:2238–53. https://doi.org/10.1111/jcmm.16217 [PubMed]
• 4. Gambardella V, Castillo J, Tarazona N, Gimeno-Valiente F, Martínez-Ciarpaglini C, Cabeza-Segura M, Roselló S, Roda D, Huerta M, Cervantes A, Fleitas T. The role of tumor-associated macrophages in gastric cancer development and their potential as a therapeutic target. Cancer Treat Rev. 2020; 86:102015. https://doi.org/10.1016/j.ctrv.2020.102015 [PubMed]
• 5. Zhao J, Zhong S, Niu X, Jiang J, Zhang R, Li Q. The MHC class I-LILRB1 signalling axis as a promising target in cancer therapy. Scand J Immunol. 2019; 90:e12804. https://doi.org/10.1111/sji.12804 [PubMed]
• 6. Kawazoe A, Fukuoka S, Nakamura Y, Kuboki Y, Wakabayashi M, Nomura S, Mikamoto Y, Shima H, Fujishiro N, Higuchi T, Sato A, Kuwata T, Shitara K. Lenvatinib plus pembrolizumab in patients with advanced gastric cancer in the first-line or second-line setting (EPOC1706): an open-label, single-arm, phase 2 trial. Lancet Oncol. 2020; 21:1057–65. https://doi.org/10.1016/S1470-2045(20)30271-0 [PubMed]
• 7. Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, Liu XQ, Sher X, Jung H, Lee M, Lee S, Park SH, Park JO, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018; 24:1449–58. https://doi.org/10.1038/s41591-018-0101-z [PubMed]
• 8. Wang F, Wei XL, Wang FH, Xu N, Shen L, Dai GH, Yuan XL, Chen Y, Yang SJ, Shi JH, Hu XC, Lin XY, Zhang QY, et al. Safety, efficacy and tumor mutational burden as a biomarker of overall survival benefit in chemo-refractory gastric cancer treated with toripalimab, a PD-1 antibody in phase Ib/II clinical trial NCT02915432. Ann Oncol. 2019; 30:1479–86. https://doi.org/10.1093/annonc/mdz197 [PubMed]
• 9. Alptekin A, Parvin M, Chowdhury HI, Rashid MH, Arbab AS. Engineered exosomes for studies in tumor immunology. Immunol Rev. 2022; 312:76–102. https://doi.org/10.1111/imr.13107 [PubMed]
• 10. Kong X, Fu M, Niu X, Jiang H. Comprehensive Analysis of the Expression, Relationship to Immune Infiltration and Prognosis of TIM-1 in Cancer. Front Oncol. 2020; 10:1086. https://doi.org/10.3389/fonc.2020.01086 [PubMed]
• 11. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m⁶A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020; 19:53. https://doi.org/10.1186/s12943-020-01170-0 [PubMed]
• 12. Segal BH, Giridharan T, Suzuki S, Khan ANH, Zsiros E, Emmons TR, Yaffe MB, Gankema AAF, Hoogeboom M, Goetschalckx I, Matlung HL, Kuijpers TW. Neutrophil interactions with T cells, platelets, endothelial cells, and of course tumor cells. Immunol Rev. 2023; 314:13–35. https://doi.org/10.1111/imr.13178 [PubMed]
• 13. Bousso P. Diving into the mechanism of action of tumor immunotherapies with intravital imaging. Immunol Rev. 2022; 306:218–23. https://doi.org/10.1111/imr.13032 [PubMed]
• 14. Neel DV, Basu H, Gunner G, Chiu IM. Catching a killer: Mechanisms of programmed cell death and immune activation in Amyotrophic Lateral Sclerosis. Immunol Rev. 2022; 311:130–50. https://doi.org/10.1111/imr.13083 [PubMed]
• 15. Tang R, Xu J, Zhang B, Liu J, Liang C, Hua J, Meng Q, Yu X, Shi S. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol. 2020; 13:110. https://doi.org/10.1186/s13045-020-00946-7 [PubMed]
• 16. Wu D, Wang S, Yu G, Chen X. Cell Death Mediated by the Pyroptosis Pathway with the Aid of Nanotechnology: Prospects for Cancer Therapy. Angew Chem Int Ed Engl. 2021; 60:8018–34. https://doi.org/10.1002/anie.202010281 [PubMed]
• 17. Ju X, Yang Z, Zhang H, Wang Q. Role of pyroptosis in cancer cells and clinical applications. Biochimie. 2021; 185:78–86. https://doi.org/10.1016/j.biochi.2021.03.007 [PubMed]
• 18. Siska PJ, Singer K, Evert K, Renner K, Kreutz M. The immunological Warburg effect: Can a metabolic-tumor-stroma score (MeTS) guide cancer immunotherapy? Immunol Rev. 2020; 295:187–202. https://doi.org/10.1111/imr.12846 [PubMed]
• 19. Rosenbaum SR, Wilski NA, Aplin AE. Fueling the Fire: Inflammatory Forms of Cell Death and Implications for Cancer Immunotherapy. Cancer Discov. 2021; 11:266–81. https://doi.org/10.1158/2159-8290.CD-20-0805 [PubMed]
• 20. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017; 541:321–30. https://doi.org/10.1038/nature21349 [PubMed]
• 21. Dubyak GR, Miller BA, Pearlman E. Pyroptosis in neutrophils: Multimodal integration of inflammasome and regulated cell death signaling pathways. Immunol Rev. 2023; 314:229–49. https://doi.org/10.1111/imr.13186 [PubMed]
• 22. Zhu X, Chen L, Liu L, Niu X. EMT-Mediated Acquired EGFR-TKI Resistance in NSCLC: Mechanisms and Strategies. Front Oncol. 2019; 9:1044. https://doi.org/10.3389/fonc.2019.01044 [PubMed]
• 23. Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia-Ramentol J, Iglesias M, Sevillano M, Ibiza S, Cañellas A, Hernando-Momblona X, Byrom D, Matarin JA, Calon A, et al. TGFβ drives immune evasion in genetically reconstituted colon cancer metastasis. Nature. 2018; 554:538–43. https://doi.org/10.1038/nature25492 [PubMed]
• 24. Desbois M, Wang Y. Cancer-associated fibroblasts: Key players in shaping the tumor immune microenvironment. Immunol Rev. 2021; 302:241–58. https://doi.org/10.1111/imr.12982 [PubMed]
• 25. Zeng QZ, Yang F, Li CG, Xu LH, He XH, Mai FY, Zeng CY, Zhang CC, Zha QB, Ouyang DY. Paclitaxel Enhances the Innate Immunity by Promoting NLRP3 Inflammasome Activation in Macrophages. Front Immunol. 2019; 10:72. https://doi.org/10.3389/fimmu.2019.00072 [PubMed]
• 26. Wu M, Wang Y, Yang D, Gong Y, Rao F, Liu R, Danna Y, Li J, Fan J, Chen J, Zhang W, Zhan Q. A PLK1 kinase inhibitor enhances the chemosensitivity of cisplatin by inducing pyroptosis in oesophageal squamous cell carcinoma. EBioMedicine. 2019; 41:244–55. https://doi.org/10.1016/j.ebiom.2019.02.012 [PubMed]
• 27. Fan JX, Deng RH, Wang H, Liu XH, Wang XN, Qin R, Jin X, Lei TR, Zheng D, Zhou PH, Sun Y, Zhang XZ. Epigenetics-Based Tumor Cells Pyroptosis for Enhancing the Immunological Effect of Chemotherapeutic Nanocarriers. Nano Lett. 2019; 19:8049–58. https://doi.org/10.1021/acs.nanolett.9b03245 [PubMed]
• 28. Li C, Qiu J, Xue Y. Low-dose Diosbulbin-B (DB) activates tumor-intrinsic PD-L1/NLRP3 signaling pathway mediated pyroptotic cell death to increase cisplatin-sensitivity in gastric cancer (GC). Cell Biosci. 2021; 11:38. https://doi.org/10.1186/s13578-021-00548-x [PubMed]
• 29. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016; 44:e71. https://doi.org/10.1093/nar/gkv1507 [PubMed]
• 30. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012; 28:882–3. https://doi.org/10.1093/bioinformatics/bts034 [PubMed]
• 31. Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics. 2013; 14:244. https://doi.org/10.1186/1471-2105-14-244 [PubMed]
• 32. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, Jensen LJ, von Mering C. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017; 45:D362–8. https://doi.org/10.1093/nar/gkw937 [PubMed]
• 33. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
• 34. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015; 1:417–25. https://doi.org/10.1016/j.cels.2015.12.004 [PubMed]
• 35. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
• 36. Jia Q, Wu W, Wang Y, Alexander PB, Sun C, Gong Z, Cheng JN, Sun H, Guan Y, Xia X, Yang L, Yi X, Wan YY, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun. 2018; 9:5361. https://doi.org/10.1038/s41467-018-07767-w [PubMed]
• 37. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
• 38. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
• 39. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
• 40. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
• 41. Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, Bin J, Liao Y, Rao J, Liao W. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res. 2019; 7:737–50. https://doi.org/10.1158/2326-6066.CIR-18-0436 [PubMed]
• 42. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, Nordgren H, Farmer P, Praz V, Haibe-Kains B, Desmedt C, Larsimont D, Cardoso F, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006; 98:262–72. https://doi.org/10.1093/jnci/djj052 [PubMed]
• 43. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
• 44. Şenbabaoğlu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, Miao D, Ostrovnaya I, Drill E, Luna A, Weinhold N, Lee W, Manley BJ, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016; 17:231. https://doi.org/10.1186/s13059-016-1092-z [PubMed]
• 45. Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, Dawson N, O'Donnell PH, Balmanoukian A, Loriot Y, Srinivas S, Retz MM, Grivas P, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet. 2016; 387:1909–20. https://doi.org/10.1016/S0140-6736(16)00561-4 [PubMed]
• 46. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, Jhunjhunwala S, Banchereau R, Yang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018; 554:544–8. https://doi.org/10.1038/nature25501 [PubMed]
• 47. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014; 9:e107468. https://doi.org/10.1371/journal.pone.0107468 [PubMed]
• 48. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, Seja E, Lomeli S, Kong X, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2016; 165:35–44. https://doi.org/10.1016/j.cell.2016.02.065 [PubMed]
• 49. Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, Zimmer L, Gutzmer R, Satzger I, Loquai C, Grabbe S, Vokes N, Margolis CA, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019; 25:1916–27. https://doi.org/10.1038/s41591-019-0654-5 [PubMed]