Gene expression profiles and clinical information were downloaded from TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases. The TCGA–BRCA cohort contained 1,109 BRCA samples and 113 normal tissue samples, and the GSE42568 cohort contained 104 BRCA samples and 17 normal tissue samples. For TCGA–BRCA cohort, fragments per kilobase million (FPKM) values were transformed into transcripts per million (TPM). The “sva” package of R software was used to address the batch effect. Patients whose survival information was unavailable were excluded from the study. A total of 1,200 patients were included in the study. The patients’ clinical information is given in Supplementary Table 1 .

The pyroptosis-related literature revealed 33 pyroptosis-related genes ( Supplementary Table 2 ) (20–22). Unsupervised clustering analysis was applied to identify distinct pyroptosis patterns based on the expression of 33 pyroptosis-related genes and classify the patients (TCGA–BRCA cohort and GSE42568 cohort) for further analysis. We used the R package “ConsensuClusterPlus” to perform the above analysis and 1000 times repetitions for guaranteeing the stability of clustering. The optimal number of clusters was determined according to the consensus clustering algorithm.

To study the differences in biological processes responsible for the characteristic patterns of pyroptosis, the “GSVA” R package was used to perform gene set variation analysis (GSVA) (23, 24). The “clusterProfiler” package was used for functional annotation and the gene set file (c2.cp.kegg.v7.2.symbols.gmt) was obtained from the MSigDB database (https://www.gsea-msigdb.org) (25–27).

We used the single-sample gene-set enrichment analysis (ssGSEA) program to measure the relative abundance of each cell infiltration in the BRCA TME. The ImmuneScore, StromalScore, and ESTIMATEScore were calculated using the “ESTIMATE” package (28–31).

A scoring system was established to quantitatively evaluate pyroptosis in individual BRCA patients. The process of establishing the scoring system is as follows: The DEGs identified from different pyroptosisclusters were firstly normalized among all samples and the overlap genes were extracted. Differential analysis and Venn diagram showed that there are 8 common differential genes among the three PyroptosisCluster. Then, we performed univariate Cox regression analysis for each gene. These genes with a significant prognosis were extracted for the next step of the analysis. We performed principal component analysis (PCA) to calculate pyroptosis scores using the following formula:

where i is the expression of pyroptosis-related genes (32, 33).

Immunophenoscore (IPS) is a good predictor of CTLA-4 and PD-1 responsiveness and thus the response to immunotherapy. The immune checkpoint inhibitor (ICI) immunophenoscore file was downloaded from the Cancer Immunome Database (TCIA, https://tcia.at/home) (34, 35).

All statistical analyses were conducted using the R (R 4.1.0) software. The student’s t-test (unpaired, two-tailed) was used to evaluate the differences between the two independent groups. One-way analysis of variance (ANOVA) and Kruskal–Wallis test were used as parametric and non-parametric methods, respectively, for data from more than two groups. The best cut-off score between the two groups with high and low PyroptosisScore was derived using the surv-cutpoint function. We applied the “limma” R package to identify the differentially expressed genes (DEGs) (36, 37). The enrichment analysis of gene ontology (GO), functional annotation, and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed using the “clusterProfiler” package. R packages “survival” and “survminer” were used for survival analysis (38, 39). The mutation landscape in patients was shown using the waterfall function of the “maftools” package (40, 41). If not specified above, P  < 0.05 was considered significant.

Figure 1 summarizes the common pathways of pyroptosis. We analyzed the incidence of copy number variation and somatic mutations of pyroptosis-related genes in BRCA. Out of 986 samples, only 95 had mutations in pyroptosis-related genes, with a mutation frequency of 9.63% ( Figure 2A ). The mutation frequency of CASP8 was the highest (2%), and there were significant differences in the expression of certain genes (ELANE, CASP5, CASP1, TNF, N0D2, IL18, NLRP7, NLRP3, and IL1B) between CASP8 wild and CASP8 mutation groups. These results indicated that CASP8 has a mutational co-occurrence relationship with these genes ( Supplementary Figure 1 ).

The expression analysis results showed that compared with the control group, the expression of CASP3, CASP6, CASP8, GSDMD, IL18, NOD2, and PYCARD was significantly higher in the tumor group, whereas that of CASP1, CASP4, ELANE, GPX4, GSDMB, GSDMC, GSEMD, IL1B, IL6, NLRP1, NLRP3, NOD1, PJVK, and SCAF11 was significantly lower ( Figure 2B ). The results of survival analysis showed that the dysregulation of most pyroptosis-related genes was significantly related to prognosis ( Supplementary Figure 2 ). The high expression of CASP6, CASP5, TIRAP, SCAF11, NLRP7, PLCG1, GSDMC, GSDMD, and NLRC4 was associated with poor prognosis. The high expression of ELANE, CASP9, CASP8, GSDMB, CASP4, CASP1, TNF, NOD1, PYCARD, NLRP6, NLRP3, NLRP2, IL6, NLRP1, IL18, and IL1B was associated with a better prognosis.

The location of CNV (copy number variations) alteration of pyroptosis-related genes on chromosomes is shown in Figure 2C . The CNV alteration in pyroptosis-related genes was mostly related to amplification in the copy number, whereas CASP3, CASP1, CASP4, NLRP6, CASP5, TIRAP, IL18, GPX4, ELANE, and CASP9 had a widespread frequency of CNV deletion ( Figure 2D ).

Our analysis cohort consisted of two BRCA datasets (TCGA–BRCA, GSE42568), OS data, and clinical information ( Supplementary Table 1 ). A univariate Cox regression analysis used to screen pyroptosis-related genes associated with prognosis in BRCA showed that IL18, CASP1, CASP4, and NLRP3 significantly correlated with prognosis ( Supplementary Table 3 ). Based on the expression of pyroptosis-related genes, the “ConsensusClusterPlus” R package was used to classify patients with qualitatively different pyroptosis modification patterns. The three different modification patterns were determined using an unsupervised clustering analysis ( Supplementary Figure 3 ) as PyroptosisCluster A, PyroptosisCluster B, and PyroptosisCluster C ( Figure 3A ).

We used the ssGSEA algorithm to analyze the immune cell infiltration of three clusters. We found that PyroptosisCluster C was remarkably rich in innate immune cell infiltration including B cells, CD4 T cells, CD8 T cells, macrophages, eosinophils, mast cells, MDSCs, mast cells, and T helper cells ( Figure 3B ). Moreover, we found that patients in the PyroptosisCluster C had a survival advantage ( Figure 3G ). The results of GSVA enrichment analysis showed that PyroptosisCluster A was significantly enriched in matrix pathways, such as ECM receptor interaction, cell adhesion molecules (CAMs) ( Figure 3C ). PyroptosisCluster C presented enrichment pathways associated with immune activation including the T cell receptor signaling pathway, B cell receptor signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, chemokine signaling pathway, cytokine receptor interactions, and JAK/STAT signaling pathway ( Figure 3E ), whereas PyroptosisCluster B was prominently related to immune suppression ( Figure 3D ). In addition, we found significant differences in the transcription profiles of pyroptosis-related genes between the three different PyroptosisClusters ( Figure 3F ). Based on the above analysis, we classified the PyroptosisCluster A as an immune rejection phenotype, characterized by innate immune cell infiltration and matrix activation. The PyroptosisCluster B was classified as an immune desert phenotype, characterized by immunosuppression. The PyroptosisCluster C was classified as an immunoinflammatory phenotype, characterized by adaptive immune cell infiltration and immune activation.

To further study the potential biological characteristics of pyroptosis-related genes, we identified eight overlapping genes, namely NLRP2, GPR132, HLA-E, FXYD5, CDH3, HLA-F, PML, and MSN in the three subtypes (PyroptosisCluster A, PyroptosisCluster B, and PyroptosisCluster C) and performed KEGG and GO enrichment analyses ( Figure 4A ). These genes showed enrichment of biological processes significantly related to pyroptosis and immunity, confirming that pyroptosis plays a crucial role in the immune regulation of TME ( Figures 4B, C ). To further verify this regulatory mechanism, we performed an unsupervised cluster analysis using the obtained genes. The results were consistent with the clustering grouping of PyroptosisCluster ( Supplementary Figure 4 ). The unsupervised clustering algorithm revealed three different pyroptosis genome phenotypes, termed gene clusters A–C. We observed significant differences in the expression of pyroptosis-related genes in three gene clusters, which was consistent with the expected results ( Figure 4D ). The results of survival analysis showed that gene cluster A showed a significant survival advantage, and gene cluster B was significantly related to poor prognosis ( Figure 4E ). The heat map shows clinical characteristics of PyroptosisCluster and gene cluster ( Figure 4F ). The opposite characteristics were observed in gene cluster A and gene cluster B.

To further analyze the functions of pyroptosis in BRCA, we constructed a scoring system, termed as PyroptosisScore, based on these pyroptosis genes to quantify the pyroptosis pattern of individual BRCA patients. We found that the majority of immune cells were negatively correlated with PyroptosisScore ( Figure 5A ). The alluvial diagram was used to visualize the changes in the attributes of individual patients ( Figure 5B ). To study the relationship between PyroptosisScore and patients’ prognosis, the “survminer” program was used to find the best cut-off value. The patients were divided into high and low PyroptosisScore groups. The high PyroptosisScore group was associated with a poor prognosis ( Figure 5C ). In addition, both gene cluster B ( Figure 5D ) and PyroptosisCluster B ( Figure 5E ) had a high PyroptosisScore, whereas both gene cluster B ( Figure 4E ) and PyroptosisCluster B ( Figure 3G ) had a poor prognosis. This finding also validated our analysis results. The results of the nomogram plot showed that PyroptosisScore may have a good advantage in long-term survival prediction ( Figure 5F ). The calibration chart showed that PyroptosisScore had a good performance, with a harmonious consistency (C index = 0.69) between the predicted and observed survival rates ( Figure 5G ).

The clinical correlation analysis revealed that patients older than 55 years had a higher PyroptosisScore ( Figure 6A ). In the patient subgroup, the high PyroptosisScore of types <55 age, T1-2, N1-3, M0, and stage I–II patients was significantly correlated with poor prognosis ( Figures 6A–E ). The R package “maftools” was used to visualize the differences in the distribution of somatic mutations in the high ( Figure 7A ) and low ( Figure 7B ) PyroptosisScore groups. We found that the mutation rate of the two was similar (83.03% and 87.59%). The results of survival analysis showed that the higher tumor mutation burden (TMB) had a poor prognosis ( Figure 7C ). Among these, H-TMB+H-PyroptosisScore had the worst prognosis ( Figure 7D ). In addition, no significant differences in TMB between the high and low PyroptosisScore groups were observed ( Figure 7E ). Furthermore, no significant correlation between PyroptosisScore and DNAss ( Figure 7F ) and RNAss ( Figure 7G ) was observed.

We further analyzed the characteristics of TME of patients with BRCA and found that a high ImmuneScore was significantly associated with a better prognosis ( Figure 8A ). No significant correlation between ImmuneScore ( Figure 8B ), ImmuneScore ( Figure 8C ), TumorPurity, ( Figure 8D ), and the prognosis was present. Except for TumorPurity, ImmuneScore, ImmuneScore, and ImmuneScore, all were significantly expressed in the low PyroptosisScore group ( Figure 8E ). In addition, we found that the low PyroptosisScore group had a more abundant immune cell infiltration pattern ( Figure 8F ). These findings emphasized the impact of pyroptosis on the BRCA microenvironment and validated our conclusion that the low PyroptosisScore group had more survival advantages than the high PyroptosisScore group ( Figure 5C ). The results of correlation analysis showed that PyroptosisScore was significantly and positively correlated with macrophages M2, macrophages M0, and resting mast cells. A significant negative correlation was observed with macrophages M1, regulatory T cells, gamma delta T cells, follicular helper T cells, CD8⁺ T cells, and CD4⁺ memory-activated T cells ( Figure 8G ).

Immunotherapy, especially immune checkpoint blocking therapy, has recently emerged as an important player in the treatment of tumors. We analyzed certain key immune checkpoint genes (PDCD1, PDCD1LG2, HAVCR2, IDO1, CD274, and CTLA4) in BRCA. The expression of these six immune checkpoint genes was significantly higher in the low PyroptosisScore group than in the high PyroptosisScore group ( Figure 9A ). In addition, PyroptosisScore had a significant negative correlation with these six immune checkpoint genes ( Figure 9B ). These six immune checkpoint genes were positively correlated with each other ( Figure 9C ). The results of GSEA enrichment analysis showed that the low PyroptosisScore group was enriched in cancer and immune-related pathways, such as antigen processing and presentation, apoptosis, B cell receptor signaling pathway, chemokine signaling pathway, cytokine receptor interactions, JAK/STAT signaling pathway, MAPK signaling pathway, natural killer cell-mediated cytotoxicity, cancer pathways, primary immunodeficiency, T cell receptor signaling pathway, and VEGF signaling pathway ( Figure 9D ). The ICI therapy represented by PD-L1 and PD-1 blockade is an effective method to treat certain tumors. To further evaluate the application of PyroptosisScore in BRCA, we obtained the immunotherapy profile of BRCA patients from the TCIA database and found that the low PyroptosisScore group had a higher ICI score and was more sensitive to immunotherapy than the high PyroptosisScore group ( Figure 9E ). In summary, tumors patients with low PyroptosisScore were characterized by inflammation, abundant immune infiltration, high expression of immune checkpoints, and a better prognosis and response to immunotherapy.