Transcriptome profiling data (FPKM value) of 1217 BC samples with their corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA) data sets (https://portal.gdc.cancer.gov/). Through screening the “negative HER2, PR, and ER” status based on clinical data, we finally identified 119 patients with TNBC and comprehensive clinical information. Other microarray data sets of 819 patients with TNBC and prognostic information were also downloaded from the Gene Expression Omnibus (GEO) data sets (https://www.ncbi.nlm.nih.gov/geo/), including 120 ER-negative BC in GSE16446, 198 HER2-negative BC in GSE25065, 310 HER2-negative BC in GSE25055, 107 TNBC in GSE58812, and 84 TNBC in GSE157284. In addition, corresponding miRNA sequencing data sets and mutation files were also obtained from the TCGA-BRCA to investigate the miRNA regulatory mechanism, and transcriptome profiling of 179 normal breast tissues was obtained from the Genotype-Tissue Expression (GTEx) database as normal controls (Carithers and Moore, 2015). The “ComBat” algorithm of the “sva” package was further applied to remove the nonbiological technical biases due to batch effects between different data sets (Leek et al., 2012). To remove the false positives caused by batch effects, we selected several stable internal reference genes (HPRT1, PPIA, RPS13, TBP, GAPDH, and HMBS) to perform the PCA analysis, which are reported as valid reference genes for human BC cell lines by Liu et al. (2015). Moreover, the IMvigor210 data sets (Hoffman-Censits et al., 2016), including 316 metastatic urothelial carcinomas (mUCs) with immunotherapy, were applied to investigate the therapeutic reaction, and the scRNA-seq data of 1534 cells from six patients with TNBC (GSE 118389) were used to validate molecular subtypes for TNBC. Detailed information on these data sets is listed in Supplementary Table S1.

According to previous studies, the caspase family, especially caspase1/4/5/8 (CASP1/4/5/8) was reported to specifically cleave GSDMD to further activate pyroptosis (Shi et al., 2015; Orning et al., 2018). In addition, Zhang et al. (2020), also found that CASP3 and granzyme B (GZMB) could convert cell apoptosis into pyroptosis through cleaving gasdermin E (GSDME). Granzyme A (GZMA) was also considered to be essential in inducing cell pyroptosis by cleaving gasdermin B (GSDMB) (Zhou et al., 2020). Moreover, inflammasome-associated families, such as absent in melanoma 2 (AIM2) and nucleotide-binding domain and leucine-rich repeat receptor (NLR), are demonstrated to induce the pyroptosis process through activating CASP1 and the release of IL1β and IL18 (Man and Kanneganti, 2015; Karki and Kanneganti, 2019). Therefore, based on the published research, a total of 33 pyroptosis-related signatures were chosen, and 24 genes were retained for subsequent analysis after filtering out the signatures with low expression (sum FPKM value of all samples less than one).

Based on the expression of pyroptosis-related signatures, we performed hierarchical clustering analysis and applied the “ConsensusClusterPlus” R package (Wilkerson and Hayes, 2010) to conduct unsupervised clustering based on Euclidean distance and Ward’s linkage methods 1000 repeated times to ensure the classification stability. During the process, the clusters from 2 to 9 were performed, respectively, and the optimal clustering model was determined based on the consensus cumulative distribution function (CDF) plot. Moreover, we performed multiple comparisons among different pyroptosis-subtypes, including for the tumor microenvironment (TME), prognosis, and vital clinical-pathological phenotypes to explore their characteristics. The R packages “survival” (Therneau and Lumley, 2015) and “survminer” (Kassambara et al., 2017) were used to perform Kaplan–Meier survival analysis and draw survival curves between pyroptosis subtypes.

To identify the DEGs between pyroptosis subtypes, the empirical Bayesian algorithm was applied through the “Limma” R package (Smyth, 2005), and the significance cutoff was set as adjusted p < .05 and absolute fold-change >1. To clarify the biological function and characteristics of pyroptosis clusters, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed by using the “ClusterProfiler” R package (Yu et al., 2012), and the results were visualized using the “ClueGO” plugin in Cytoscape v3.7.1 (Bindea et al., 2009).

To evaluate the immune cell infiltration (ICI) characteristics of TNBC tissues, we used the “CIBERSORT” R package (Chen et al., 2018) to quantitatively analyze the infiltration levels of different immune cells with the LM22 signatures by 1000 random permutations. The tumor purity scores, ICI levels, and stromal contents in different samples were evaluated via the “ESTIMATE” algorithm (Yoshihara et al., 2013). Moreover, through the “c2. cp.kegg.v6.2. symbols” data sets downloaded from the MSigDB database, we performed GSVA using the “GSVA” R package and drew a heatmap to exhibit the different immunogenic pathways (Hanzelmann et al., 2013).

To further identify a novel index representing the characteristics of the pyroptosis subtypes, we conducted univariate Cox proportional hazards regression analysis for overall survival (OS) to preliminarily screen significant genes through using the “coxph” function in the “survival” R package. Subsequently, to remove the multicollinearity among these candidate genes, LASSO regression was applied to screen independent prognosis-related genes with the optimal penalty parameter and a minimum 10-fold cross-validation (Ranstam and Cook, 2018). After further adjustment, the multivariate Cox regression (stepwise model) was conducted to identify hub genes, and the coefficients obtained from the regression algorithm were used to acquire the Ps-score based on the following formula: Ps _ score = val ( Gene 1 ) × β 1 + val ( Gene 2 ) × β 2 + … + val ( Gene  n ) × βn . The val (Gene) represents the expression FPKM value of each gene and β the corresponding regression coefficient. Moreover, according to the above formula, the Ps-scores of patients with TNBC were separately calculated, and the patients were divided into high and low subgroups according to the median value as the cutoff value (Sullivan et al., 2004). We also made similar comparisons between high and low Ps-score groups, including the TME, ICI, clinical phenotypes, pyroptosis-related signatures, and correlation of GSVA pathways based on the other four GEO data sets described earlier.

The multivariate Cox regression (stepwise model) was applied to construct the prognostic model for TNBC-combined Ps-scores and other clinical features, including age, clustering subtypes, clinical stages, and TNM stages. Variables with p-values < .05 were included in the Cox regression model, and the nomogram was further constructed to predict the probability of 1-, 3-, and 5-year survival in patients with TNBC using the “survival” R package. To evaluate and validate the prediction capability of the nomogram, we calculated the concordance index and plotted the calibration curves for 3- and 5-year survival through a bootstrapping method with 1000 resamples. To further investigate the expression of the Ps-score-related genes at the protein level, the Human Protein Atlas (HPA) (Ponten et al., 2008) was used to display the results of the immunohistochemistry (IHC) technique. The detailed information of patients is included in Supplementary Table S17.

Based on the IMvigor210 data sets with atezolizumab treatment, we performed a comprehensive comparison between different Ps-score subgroups, including response to immunotherapy, immune phenotype, and clinical remission rate. Moreover, to evaluate the potential therapeutic value of Ps-scores in chemotherapy for TNBC, we calculated the half-maximal inhibitory concentration (IC₅₀) of common chemotherapeutic drugs based on the Genomics of Drug Sensitivity in Cancer (GDSC) databases (Yang et al., 2013). Antitumor drugs such as 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel are recommended for BC treatment by current clinical guidelines. Differences in IC₅₀ of these chemotherapeutic drugs between Ps-score subgroups were compared by Wilcoxon test with the results exhibited in box diagrams using the “ggpubr” R package (Whitehead et al., 2019).

To validate molecular clusters and further seek biomarkers of each cluster, the Seurat pipeline was selected for subsequent analysis. Using the Seurat package v3.0 (Butler et al., 2018), we transformed the data matrix into a “Seurat object” through the “CreateSeuratObject” R function and performed the necessary quality control. The violin diagram exhibited the number of sequencing reads per sample and the expression of mitochondrial genes (Supplementary Figure S2A). Further, to remove the influence from mitochondrial and extreme genes, we kept the number of sequenced genes at 200–10,000, directly including the majority genes, and removed the cells with average gene expression <10 and mitochondrial genes >5%. Then, we conducted data standardization through the “NormalizeData” function with the method of “LogNormalize” and used the top 1500 variable counts to perform PCA using the “RunPCA” function (Supplementary Figures S2B,C). Subsequently, t-distributed statistical neighbor embedding (tSNE) was applied to visualize the density clustering, and the “SingleR” package was applied for cell-subtype annotation based on the marker genes (Aran et al., 2019). Moreover, the pseudotime trajectory analysis was further performed using the “Monocle” package v2.0 to expound the potential inner relationship among these cell clusters (Trapnell et al., 2014).

The miRNAs targeting Ps-score-related genes were predicted based on the following databases: TargetScan (http://www.targetscan.org/), starBase (http://starbase.sysu.edu.cn/starbase2/index.php), miRTar (https://mirtarbase.cuhk.edu.cn/), and miRDB (http://www.mirdb.org/). The expression of TNBC miRNA was downloaded from the TCGA-BRCA data sets, and the DEmiRNAs were further identified by the “Limma” R package. Subsequently, we identified the intersection of predicted miRNA by four databases and DEmiRNAs as regulated miRNAs for each hub gene and further visualized the miRNA–mRNA interaction network using Cytoscape v3.7.1.

After a series of rigorous screening and quality control steps, a total of 24 pyroptosis-related signatures remained for subsequent analysis in our study (Supplementary Table S2). Combined with the normal tissues in GTEx data sets, we first compared the expression of pyroptosis-related signatures between patients with TNBC and normal controls. We found that most pyroptosis-related genes were significantly upregulated in TNBC groups, including CASP1/3/5/8, GSDMA/C, NLRC4/P3/P7, IL18, IL1β, and TNF (Figure 1A). PCA indicated that the expression of these pyroptosis-related signatures could be used to divide the TNBC samples and controls into two distinct clusters (Figure 1B). In addition, the PCA of internal reference genes revealed that nonsignature genes failed to discriminate TNBC and control cohorts, indicating that the separation created a true distinction based on the pyroptosis-related DEGs in patients with TNBC (Supplementary Figure S1A). Moreover, the KEGG functional enrichment analysis revealed that these pyroptosis-related genes were predominantly focused on infectious diseases, immune response, and cellular signal conditioning mechanisms, including the NLR signaling pathway, p53 signaling pathway, TNF signaling pathway, and apoptosis (Figure 1C, Supplementary Table S3). In terms of genetics, 148 of the 203 samples (72.91%) manifested pyroptosis-related signatures in mutations and the NLR families, especially NLRP3 and NLRP7, exhibiting the highest frequency of mutations (Figure 1D, Supplementary Table S4). Moreover, the top 10 pyroptosis-related genes with the most frequent mutations were located on the 24 human chromosomes (Figure 1E).

Based on the expression of pyroptosis-related signatures, we used an unsupervised clustering method to identify the subtypes of patients with TNBC and identified k = 2 as the optimum clustering model from k = 2 to k = 9 clustering with the least area under the consensus CDF curve for 69 patients in clusterA and 50 patients in clusterB (Figure 2A, Supplementary Table S4; Supplementary Figure S1B). To further clarify the intrapatient heterogeneity of patients with TNBC, we performed a comparison in the clinical differences between subtypes and found that patients in clusterB were negatively associated with severe clinical stages, including the pathological and TMN stages. Furthermore, the survival analysis showed that patients in clusterB had a longer median survival time than those in clusterA with more patients in clusterB also receiving radioactive treatments, indicating that patients with TNBC in clusterB might have a better prognosis (Figure 2B). There was no significant difference in other clinical indexes including age, sex, M stage, and the ratio of pharmaceutical and surgical therapies (Supplementary Figure S1C). Notably, the expression of pyroptosis-related gene signatures was significantly increased in patients in clusterB compared with that of the clusterA cohort (Figure 2C). In terms of the immune infiltration scores, adaptive immune response-related lymphocytes (including memory B cells, activated memory CD4⁺ T cells, plasma cells, CD8⁺ T cells, and gamma delta T cells) were significantly increased in patients in clusterB compared with the clusterA cohort. However, innate immunity and immunoregulation-related cells were significantly infiltrated in clusterA cohorts, including neutrophils, activated natural killer cells, resting memory CD4⁺ T cells, and regulatory T cells (Tregs) (Figure 2D, Supplementary Table S5). Higher stromal scores and immune scores with lower tumor purity were also detected in patients in clusterB compared with the clusterA groups (Figures 3A,B).

Considering the biological characteristics of immune subtypes in TNBC, we conducted a DEG analysis between the two subtypes. Through comparing clusterA with clusterB groups, a total of 262 DEGs (including 13 clusterA- and 249 clusterB-related genes) in TNBC were identified (Figure 3C, Supplementary Table S6). To further interpret biological processes and pathways of pyroptosis-related subtypes, these DEGs were chosen to perform KEGG functional analysis. The results showed that clusterA-related genes were not enriched in any significant pathways while the clusterB-related signatures were predominantly enriched in immune activation-associated pathways, including natural killer cell–mediated cytotoxicity, the toll-like receptor signaling pathway, chemokine signaling pathway, cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, TNF signaling pathway, Th17 cell differentiation, and Th1 and Th2 cell differentiation (Figure 3D, Supplementary Table S7).

After successively including the 262 DEGs in univariate Cox regression, LASSO regression, and multivariate Cox regression analysis as candidate prognosis-associated genes, we identified six hub genes (including CFB, IFITM1, EPSTI1, MARCO, CXCL13, and CCL5) from the Ps-score signatures based on their β coefficients (Figures 3E,F, Supplementary Table S8). In addition, the expression of these hub genes was higher in clusterB subgroups. Based on the IHC data from the HPA database, the expression of these hub genes at the protein level was further validated in BC, especially CFB, ESPIT1, and IFITM1 (Figure 3G, Figure 6C). Based on the expression of these genes and their corresponding β coefficients, the Ps-score was defined by the following formula: Ps _ score = − 0.365 × C F B − 0.45 × I F I T M 1 − 0.298 × E P S T I 1 − 0.461 × M A R C O − 0.26 × C X C L 13 + 0.439 × C C L 5 (Supplementary Table S9, the gene name represents the corresponding gene expression FPKM values). Subsequently, those patients with TNBC were divided into a high and low Ps-score subgroups with median value (–5.28) as the cutoff; the high Ps-score cohorts exhibited a worse prognosis than that of low Ps-score patients in the TCGA data sets (Figure 4A). To prove the universal value of the Ps-score in TNBC, we also performed survival analysis of this score in four extrinsic GEO cohorts and obtained the same results (Figure 4B).

To investigate the biological characteristics of the Ps-scores, we compared the expression of pyroptosis-related genes and ICI between different Ps-score groups and further explored the correlation between significant clinical phenotypes and the Ps-scores. The results reveal that clusterB possessed a lower level of Ps-scores associated with pharmaceutical and surgical therapy as well as lower pathological stages in patients with TNBC (Figure 4C&E). Interestingly, the expression of pyroptosis-related signatures was significantly increased in the low Ps-score groups, including the CASP, GSDM, and NOD families, as well as inflammatory factors, suggesting the potential role of pyroptosis activation in the low-score of TNBC cohorts with better prognosis (Figure 4D). Moreover, ROC analysis showed 1-, 3-, and 5-year AUC values of the Ps-scores for predicting the prognosis of patients with TNBC of 0.867, 0.867, and 0.906, respectively, in the TCGA sets (Figure 5A). Furthermore, immune infiltration analysis revealed that substantial immune cells were significantly inhibited in the high Ps-score groups, including that of CD8⁺ T cells, follicular helper T cells, activated CD4⁺ memory T cells, and plasma cells (Figure 5D). The correlation analysis also indicates that the Ps-scores were significantly positively associated with the levels of tumor purity (R = 0.57, p < .001) but negatively associated with stromal and immune scores (R = –0.48, p < .001 and R = –0.50, p < .001, respectively) (Figure 5E).

Based on the Ps-scores and some primary clinical characteristics, multivariate Cox regression analysis was conducted to construct a nomogram that could accurately predict the probability of the 1-, 3-, and 5-year survival for patients with TNBC. The Ps-scores, age, pyroptosis-related cluster, and clinical stages were considered as related predictors for the prognosis of patients with TNBC and incorporated into the nomogram with significant regression coefficients and p-values (Figure 5B, Supplementary Table S10). From the nomogram, we could observe that the Ps-score contributed the most to the total score with a 0.716 concordance index (Figure 5B, Supplementary Table S11). Calibration curves exhibited that the nomogram had a good prediction capacity in both 3- and 5-year OS for patients with TNBC (Figure 5C).

The alluvial diagram visualized the status changes in the different characteristics of patients (Figure 6A). We found that patients with a low Ps-score in clusterB subtypes had a higher ratio of survival status. When using the TCGA and other external GEO data sets, the results of GSVA demonstrated the coincident negative correlation between Ps-scores and immunoregulation-related pathways, such as the toll-like receptor signaling pathway, antigen processing and presentation, rig I-like receptor signaling pathway, T cell receptor signaling pathway, NOD-like receptor signaling pathway, and JAK-STAT signaling pathway (Figure 6B, Supplementary Table S12).

To further explore the role of Ps-scores in predicting the therapeutic benefit in TNBC, we first calculated the Ps-scores of patients who accepted anti-PD-L1 immunotherapy from the GSE157284 and IMvigor210 cohorts before assigning them into high or low Ps-score groups. From the GSE157284 data set, patients with an effective response rate to anti-PD-L1 therapy showed lower Ps-scores, whereas the relative expression of immune-checkpoint genes was significantly increased in the high Ps-scores groups (Figures 6D,E, Supplementary Table S13). Congruously, we also found that the low Ps-score group had higher expression of immune-checkpoint genes in the IMvigor210 cohort and effective anti-PD-L1 responders also exhibited lower Ps-scores (Figures 6F,G). In metastatic urothelial cancers of the IMvigor210 data sets, distinct immunological subtypes might result in opposite therapeutic responses. Therefore, we also compared the Ps-score levels among these subtypes and found that the lowest Ps-score was in the inflamed subtypes (Figure 6H). Moreover, the rate of complete remission (CR) after immunotherapy was also increased in the low Ps-score cohort compared with the high Ps-score cohort, with the low Ps-score group validated to have a better prognosis for metastatic urothelial cancers (Figures 6I,J). All these results suggested that the Ps-score might serve as a significant indicator in immunotherapy decision making for cancers.

Besides checkpoint blocker therapy, we also investigated the potential associations between the Ps-scores and the curative efficacy of common chemotherapeutics in treating BCs. The IC₅₀ value was calculated for five common anti-BC chemotherapeutic drugs obtained from the GDSC databases, including 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel (Supplementary Table S14). Notably, all the drugs exhibited lower IC₅₀ values in the low Ps-score groups, indicating patients with low Ps-scores might obtain a better curative efficacy from common chemotherapy (Figure 6K). Collectively, these outcomes indicate that Ps-scores could be associated with the response to immunotherapy and common chemotherapy.

To validate the potential subtypes of patients with TNBC, the GSE118389 data set along with 1534 cells from six TNBC tissues were used to identify concrete cellular subtypes and corresponding marker genes. A total of 12 distinct cellular clusters were identified through tSNE analysis (Figure 7A) with the marker-genes of each cluster listed in Supplementary Table S15. Moreover, the results of cell-type annotation using the “SingleR” package indicated these cell clusters fell into six cellular subtypes, including epithelial cells, erythrocytes, CD8⁺ T cells, fibroblasts, endothelial cells, and monocytes, of which the epithelial cells were the most common cell type with six subtypes (Figure 7B). Notably, of the six Ps-score-related genes, five genes were identified as marker genes, and their expression in each cellular subtype is shown in a bubble diagram (Figure 7C). The expression of CCL5 and IFITM1 were generally increased in nearly all cellular subtypes and significantly higher than the expression of other signatures. Although the expression of remanent genes was relatively lower in these cells, significant cellular specificity was found in these Ps-score-related genes. For example, MARCO was particularly expressed in the epithelial cell subtype 4 and monocytes while CFB was particularly expressed in the fibroblast subtype 2 as well as epithelial cell subtypes 3 and 4. Interestingly, EPSTI1 was significantly expressed in immune-related cellular subtypes including monocytes and CD4⁺ T cells, consistent with the results of IHC. In addition, the pseudotime trajectory analysis also revealed the arrangement of different cellular subtypes that formed a certain rule based on its spatial relationships (Supplementary Figure S2E). Concretely, the trajectory analysis of epithelial cells revealed that a small quantity of epithelial cell subtype 3 was distributed at the start of the trajectory while mixed cells from epithelial cell subtypes 5 and 6 were distributed at the end of the trajectory. Moreover, epithelial cell subtype 2 was uniformly located behind epithelial cell subtypes 1, whereas epithelial cell subtype 4 nearly existed throughout the trajectory (Figure 7D).

To further explore the potential regulatory role of miRNAs targeted to these Ps-score-related genes in TNBC, we found 223 probable miRNAs and successfully constructed the miRNA-hub gene interaction network based on the prediction of the TargetScan, starBase, miRTar, and miRDB databases (Figure 7E, Supplementary Table S16). In addition, to validate the regulatory role of miRNAs in TNBC, we also identified 37 DEmiRNAs between pyroptosis-related clusters, before screening the top 20 miRNAs serving as vital regulatory factors (Figures 7F,G). Finally, we successfully simplified the interaction network with 20 miRNAs and four targeted genes to verify that ESPIT1 and CXCL13 are the most active targets regulated by massive DEmiRNAs (Figure 7H).