The Cancer Genome Atlas (TCGA) genomic data for LUAD samples were obtained from the Genomic Data Commons. Gene Expression Omnibus (GEO) was used to download gene expression microarrays of LUAD samples (GSE31210) and non-small cell lung cancer (NSCLC) samples (GSE37745 and GSE50081) and lung cancer (GSE30219). The Robust Multichip Average (RMA) method and R package “affy” normalized GSE37745 gene expression data. Detailed information of the cohorts is presented in Supplementary Tables S1, S2.

IMvigor210 was a single-arm phase Ⅱ study that looked into an anti-PD-L1 agent (atezolizumab) in patients with metastatic urothelial carcinoma (mUCC) (NCT02108652 and NCT02951767) (Mariathasan et al., 2018). The R package “IMvigor210CoreBiologies” obtained all the expression and clinical data from the IMvigor210 trials. GEO provided RNA-seq data for a total of 27 advanced NSCLC patients who were treated with anti-PD-1/PD-L1 (GSE135222).

A total of 47 PRGs were obtained from the study of Song et al. (2021). The R package “maftools” was used to demonstrate PRG mutation. The R package “ggpubr” was used to visualize the copy number variation (CNV) information of PRGs. The R package “limma” was used to examine the differential expression of PRGs in tumor samples.

The Pathway Commons database was used to find PRG protein–protein interactions. Pearson correlation was used to examine the co-expression status of PRGs (Supplementary Table S3). Cytoscape software was utilized to visualize the correlation network.

Based on the pyroptosis genes and R package “pheatmap,” K-means clustering was used to determine the pyroptosis-related subtypes (subtypes 1 and 2). The Kaplan–Meier survival analysis was performed to analyze patient differences between the two subtypes in conjunction with the log-rank test. The difference between two subtypes based on the PRG expression was investigated using principal component analysis (PCA).

We obtained 25 cancer treatment-predicted signature sets from various publications (Sweis et al., 2016; Ayers et al., 2017; Mariathasan et al., 2018; Kamoun et al., 2020). The R package “GSVA” was used to calculate the therapeutic signature gene set enrichment score using gene set variation analysis (GSVA). Detailed information of 25 cancer treatment-predicted signature sets is listed in Supplementary Table S4. The one-sided Wilcoxon rank-sum test was used to analyze the differences in the therapeutic enrichment scores between subtypes.

The range of infiltration of 22 immune cells in TCGA LUAD samples was inferred by the CIBERSORT (Cell-type Identification by Estimating Relative Subsets of RNA Transcripts) method (Newman et al., 2015). CIBERSORT can compute the abundances of specific cell types in a mixed sample based on the bulk expression. In addition, the ESTIMATE (Estimation of STromal and Immune Cells in MAlignant Tumor Tissues Using Expression Data) method was used to calculate the abundances of immune cells by the R package “estimate.” We focused on the mRNA expression of five immune checkpoints: PD-1, PD-L1, CTLA4, CD47, and BTLA. The one-sided Wilcoxon rank-sum test was utilized to analyze the differences between subtypes.

The R package “limma” discovered 719 differentially expressed genes (DEGs) between two subtypes with |log2FC| > 0.5 and p < 0.001. A web-based program, Metascape, was used to perform the enrichment analysis on 719 DEGs using ontology sources such as KEGG Pathway, GO, Reactome, and other canonical pathways (Supplementary Table S5). Then, a selection of enriched terms with a similarity greater than 0.3 was chosen and shown as a network plot.

WGCNA (weighted gene co-expression network analysis) is a data reduction method and an unsupervised classification method (Langfelder and Horvath, 2008; Langfelder and Horvath, 2012). The co-expression network was built using the Sangerbox 3.0 tool and DEG expression profile. Module–trait association analysis was used to determine which co-expression module was the most relevant to the clinical features. The genes were clustered, and a heatmap was created to illustrate the relationship between modules and phenotype.

The least absolute shrinkage and selection operator (LASSO) approach and Cox regression model were employed to screen the prognostic genes in the key module. One standard error (SE) over the minimum threshold was chosen. The R package “glmnet” managed the entire process. Finally, a seven-gene risk score formula was developed, and multivariate Cox regression coefficients were computed using the R package “survival”: Pyroptosis subtype-related risk score (PSR_score) = (exp Gene1 * coef Gene1) + (exp Gene2 * coef Gene2) + … +(exp Gene7* coef Gene7).

Patients were classified based on the median of their PSR_score. The R package “survival” used the log-rank test to compare the survival times of patients with high PSR_score and patients with low PSR_score. Furthermore, stratified analysis was performed to determine the protective effect of PSR_score based on the T stage, N stage, M stage, and tumor stage. Chi-square tests were used to examine the connections between the PRG score and clinical factors such as age, gender, T stage, N stage, and M stage. The data were presented using Kaplan–Meier graphs (Supplementary Table S1, S2).

The one-sided Wilcoxon rank-sum test was used to determine the difference between the two subtypes or high- and low-PSR_score groups. R version 4.1.2 was used for all statistical studies. p < 0.05 was considered statistically significant.

Supplementary Figure S1 depicts the analytical process used in this study. We first explored the landscape of variation in PRGs in the genome and transcriptome. A relatively high mutation frequency of PRGs was observed in LUAD (Figure 1A). TP53 exhibited the highest mutation frequency (55%), followed by NLRP3, NLRP7, and NLRP2. Then, we looked at the link between TP53 mutation and PRG expression. CHMP7, IRF2, CASP4, ELANE, BAX, and TIRAP were all downregulated in TP53 mutation samples (Supplementary Figure S2, p < 0.1). Following that, we investigated the CNV landscape of PRGs in LUAD (Figure 1B). Copy number amplification was common in HMGB1, BAX, CASP3, IRF2, IL18, and GPX4, whereas copy number deletion was common in GSDMC, GSDMD, AIM2, and CHMP6.

Furthermore, we investigated the difference in PRG expression levels between tumor and normal tissues (Figure 1C). A total of 30 (63.83%) PRGs showed differential expression (p < 0.05), with 23 genes showing substantial upregulation and seven showing significant downregulation in tumor samples.

We built an interaction network to investigate the relationship between PRGs (Figure 2A). The color of the edges indicated the five types of protein–protein interactions, and the thickness of the edges indicated the level of co-expression between PRGs, as determined by Pearson correlation (Supplementary Table S3). The network showed a strong relationship between PRGs.

To investigate the heterogeneous features of LUAD further, a K-means clustering algorithm was used to categorize patients based on PRG expression profiles. Patients with LUAD were classified into two subtypes (Figure 2B). Survival analysis revealed that subtype 1 had a considerably greater overall survival than subtype 2 (Figure 2C, p = 0.039, log-rank test). According to principal component analysis (PCA), LUAD patients had unique PRG expression patterns between two subtypes (Figure 2D).

The therapeutic differentiation between the subtypes was investigated, and the GSVA approach was utilized to determine the score of 25 therapeutic signature sets in TCGA LUAD data (Figure 3A). A total of 23 (92%) therapeutic signatures differed significantly between the two subtypes, with 20 therapeutic signature scores in subtype 2 significantly higher than those in subtype 1 and three therapeutic signature ratings significantly lower (Figure 3B, p < 0.05). Patients in subtype 2 were found to be more amenable to treatment.

The differentiation of TME between two subtypes is then evaluated. According to the CIBERSORT algorithm, infiltration of “B cells naive,” “dendritic cells activated,” “mast cells resting,” “monocytes,” and “neutrophil plasma cells” were higher in subtype 1 than in subtype 2 (Figure 3C, p < 0.05). “B cells memory,” “macrophages M1,” “NK cells resting,”, “T cells CD4 memory activated,” and “T cells CD8” showed significantly lower infiltration in subtype 1 than in subtype 2 (Figure 3D, p < 0.05). Furthermore, we investigated the tumor purity differentiation across the subtypes, finding that the ESTIMATE score, stromal score, and immune score in subtype 1 were considerably lower than those in subtype 2 (Figure 3E, p < 0.05). Furthermore, we investigated the distinction between the subtypes in the ability to recognize tumor cells and execute immune responses. We looked at the differential expression of five immunological checkpoints and discovered that the expression of all the five immunological checkpoints was considerably greater in subtype 2 than that in subtype 1 (Figure 3F, p < 0.05). The result indicated that samples in subtype 2 had a higher level of immune infiltration.

To investigate the potential biological activity of the subtypes, we detected DEGs between the two subtypes, and Metascape performed enrichment analysis on 719 DEGs (Figures 4A,B). The DEGs were found to be significantly enriched in a variety of immune-related pathways and processes, including “leukocyte activation,” “inflammatory response,” “innate immune response,” and “positive regulation of immune response” (Supplementary Table S5).

Then, WGCNA was used to build co-expressed networks based on the expression of 719 DEGs and identify important modules linked with clinical traits. The power value for modules was screened to ensure an average connection and high independence. The power value in this study was set at 5 as the soft-thresholding parameter to ensure a scale-free network (Figures 4C,D). In total, four modules have been identified (Figure 4E). The module–trait association analysis was used to discover co-expression modules that were highly relevant to clinical traits. Figure 4F depicts the relationship between modules and phenotype. Correlation analysis revealed that the blue module, which comprises 91 genes, was identified as a correlation between the prognosis and tumor stage. The top five highly enriched phrases for blue module genes were “secretion,” “cellular-modified amino acid metabolic process,” “epidermis development,” “NABA MATRISOME ASSOCIATED,” and “malignant pleural mesothelioma” (Supplementary Figure S3; Supplementary Table S5).

A model was built with seven pyroptosis subtype-related co-expression prognostic genes, FOSL1, KRT6A, GPR133, TMPRSS2, PRDM16, SFTPB, and SFTA3, to investigate the prognostic value of the selected subtype-related co-expression blue module genes (Figures 5A,B). Then, using the expression of seven genes, we established a predictive model according to the multivariate Cox proportional hazard model: PSR_score = (0.1072 * FOSL1 exp) + (0.09327 * KRT6A exp) + (−0.1144 * GPR133 exp) + (0.04062 * TMPRSS2 exp) + (−0.1238 * PRDM16 exp) + (−0.02503 * SFTPB exp) + (−0.04079 * SFTA3 exp).

The PSR_score of each patient in TCGA was calculated using the seven-gene-involved formula. The patients were divided into two groups using the median as the cutoff value: those with a high PSR_score and those with a low PSR_score. Patients with a high PSR_score had a substantially shorter life expectancy (Figure 5C, p = 6. 8e-10, log-rank test). The area under the curve (AUC) of the receiver operating characteristic (ROC) curve revealed that PSR_score correctly predicted mortality (Figure 5D, AUC = 0.679). We then investigated PSR_score’s ability to predict patient prognosis within clinicopathological subgroups. In most cancer stages, high PSR_score patients had a substantially worse OS than low PSR_score patients (Supplementary Figure S4, p < 0.05, log-rank test).

Following that, we validated the prognosis power of PSR_score in independent datasets. Survival analysis was carried out in four GEO lung cancer cohorts (GSE30219, GSE31210, GSE37745, and GSE50081), and the results revealed that a high PSR_score indicated a poor prognosis in all GEO datasets (Supplementary Figure S5, p < 0.1, log-rank test). We combined four GEO lung cancer cohorts into a big dataset to confirm the robustness of PSR_score. Similarly, patients with a high PSR_score had a significantly poor OS (Figure 5E, p = 1.1e-16, log-rank test), with an AUC of 0.682 (Figure 5F).

Pearson correlation analysis was performed to assess the relationship between PRG_score and the number of immune cells to study the link between PRG_score and immunological infiltration. Infiltration of “macrophages M1,” “T cells CD4 memory activated,” “macrophages M0,” “NK cells resting,” “NK cells activated,” “T cells CD8,” and “dendritic cells activated” was significantly positively connected with PRG_score (Figure 6A-G, p < 0.05, Pearson correlation analysis). Furthermore, ESTIMATE score of high PRG_score samples was higher than that of low PRG_score samples (Figure 6H). We also investigated the relationship between the expression of seven genes in the model and immune cells. We discovered that the quantity of most immune cells was associated with the expression of these genes (Figure 6I). In TCGA LUAD cohorts, the expression of PD-1 and PD-L1 was significantly higher in high PRG_score samples than in low PRG_score samples (Figures 6J,K, p < 0.05).

To further explore if the risk score can predict patients' responses to immunotherapy, we compared OS of patients with a high PRG_score versus low PRG_score who were receiving immunotherapy. In IMvigor210 and GSE135222 cohorts, patients with a high PRG_score had a significantly worse prognosis (Figures 7A,B, p < 0.05, log-rank test). In addition, we looked at the differences in immune checkpoint gene expression between high and low PRG_score groups. PD-L1 and CD47 in the high PRG_score group of the IMvigor210 cohort were significantly greater than those in the low PRG_score group (Figures 7C,D; p < 0.05).