Pyroptosis‐related molecular classification and immune microenvironment infiltration in breast cancer: A novel therapeutic target

Abstract The underlying role of pyroptosis in breast cancer (BC) remains unknown. Herein, we investigated the correlations of 33 pyroptosis‐related genes (PRGs) with immune checkpoints and immune cell infiltrations in BC patients based on The Cancer Genome Atlas cohort (n = 996) and Gene Expression Omnibus cohort (n = 3,262). Enrichment analysis revealed that these PRGs mainly functioned in pyroptosis, inflammasomes and regulation of autophagy pathway. Four prognostic independent PRGs (CASP9, TIRAP, GSDMC and IL18) were identified. Then, cluster 1/2 was recognized using consensus clustering for these four PRGs. Patients from cluster 1 had a favourable prognosis and diverse immune cell infiltrations. A nomogram was developed based on age, TNM stage, tumour subtype and pyroptosis score. Patients with the high‐risk group exhibited worse 5‐year OS, and the result was consistent in the external cohort. Additionally, high‐risk group patients were associated with downregulated immune checkpoint expression. Further analysis suggested that the high‐risk group patients were associated with a higher IC50 of paclitaxel, doxorubicin, cisplatin, methotrexate and vinorelbine. In summarizing, the pyroptosis score‐based nomogram might serve as an independent prognostic predictor and could guide medication for chemotherapy. Additionally, it may bring novel insight into the regulation of tumour immune microenvironment in BC and help to achieve precision immunotherapy.

identifying those patients with high risk and predicting the overall survival and immunotherapeutic response in advanced BC patients.
Pyroptosis, a form of programmed cell death apart from apoptosis and autophagy, is featured by cytoplasmic swelling, cells lysis and release of inflammatory factors, which triggered inflammation and immune response. 10 The pyroptosis can be activated by the caspase-1/3/4/5/11 inflammasome pathway and is induced by the gasdermin family. [11][12][13][14] Recent studies on pyroptosis had found that it may be closely related to malignant tumours and play a vital role in either antitumour effect or carcinogenesis, progression and drug resistance. 15 For instance, studies had revealed that the overexpression of AIM2 inflammasome exhibited antitumour effects in hepatocellular cancer, and decreased expression or absence of AIM2 may cause colorectal carcinoma. 16,17 On the contrary, overexpression of GSDMB was reported to be associated with a worse prognosis and less sensitivity to HER2-targeted therapy in BC. 18 Besides, DFNA5/ GSDME acted as an antioncogene in hepatocellular cancer, and loss of DFNA5/GSDME contributed to drug resistance in melanoma and lung cancer. [19][20][21] Signatures based on pyroptosis showed substantial prognostic predictive values in ovarian cancer, 22 glioblastoma, 23 melanoma, 24 and gastric cancer, 25 but pyroptosis-related signatures on BC are lacking. Besides, pyroptosis was identified in many types of immune cells, but the potential function of pyroptosis in the tumour immune microenvironment (TIME) remains undefined. 26 A study demonstrated that pyroptosis-produced cytokines were released into the TIME, exhibiting tumour-promoting effects. 27 However, pyroptosisinduced inflammation in the TIME could enhance the sensitivity of immunotherapy through activating certain immune cells. 28 Overall, existing studies had confirmed that pyroptosis had an essential role in both oncogenic and antioncogenic aspects, and we may utilize it as a prognosis predictor to help in making decisions in the treatment selected.
This study aimed to comprehensively evaluate the association of pyroptosis with the prognosis, TIME and chemotherapy responsiveness in BC patients. We constructed a PRGs-based nomogram and estimated its efficacy in prognostic prediction, chemotherapeutic responsiveness and function in TIME regulation.

| Data extraction and identification of pyroptosis-related genes
Gene expression files of breast normal and tumour tissues, clinicopathological characteristics and the corresponding follow-up information of BC patients were extracted from The Cancer Genome Atlas (TCGA) data set (https://portal.gdc.cancer.gov/). Moreover, mutation information and copy number variation information for BC were also downloaded from the TCGA data set. The gene expression files and clinical information mentioned above of the external validation cohort were retrieved from the Gene Expression Omnibus (GEO) databank (https://www.ncbi.nlm.nih.gov/geo/). Analyses of these data were conducted by R (version 4.0.5) and R Bioconductor packages. Pyroptosis-related genes (PRGs) were acquired from published literatures. [29][30][31] Then, the Search Tool for the Retrieval of Interacting Genes (STRING), version 11.0 (https://strin g-db.org/), was applied to construct a protein-protein interaction (PPI) network for the PRGs.

| Functional enrichment analysis and mutation analysis of PRGs
We performed GO analysis in the heatmap and network of enriched terms on the Metascape website (http://metas cape.org) to identify the biological process function of these selected PRGs. The 'ggplot2' package in R was applied to perform KEGG analysis. GSEA analysis was achieved in the Hallmark gene set with GSEA v4.1.0 software.
Pathways were defined as enriched when p < 0.05.

| Consensus clustering for the PRGs and immune cell infiltration (ICI) analysis
The unsupervised clustering analysis based on the expression levels of PRGs was conducted by the 'Consensus Cluster Plus' R package.
The 'CIBERSORT' R package was utilized to clarify the distinct infiltration levels for immune cells and to quantify the percentage of 22 common immune cell types in the BC samples in different clusters and different risk groups. Spearman analysis was applied to investigate the correlation coefficient between any two types of immune cells.

| Development and validation of pyroptosis risk score-based nomogram
According to prior reviews, a total of 33 PRGs were identified. 15,32,33 Next, univariate and multivariate Cox regression analyses (UMCRA) were implemented to search the correlation between the PRGs expression level and 5-year OS of BC patients and to explore independent prognostic PRGs in the TCGA cohort. Specifically, univariate Cox regression analysis was utilized to screen out PRGs that were significantly related to the OS, and PRGs that meet the requirements of p < 0.05 were extracted for multivariate Cox regression analysis. Additionally, the pyroptosis risk score (PS) was established using the screening independent prognostic PRGs. The PS was calculated after centralization and standardization by applying the 'scale' function in R. Formulae are as follows: PS = expression level of PRGs × sum of coefficients. What's more, Cox analysis was utilized to screen independent clinical variables for models. Finally, a novel PS-based nomogram containing the PS, age, tumour subtype and TNM stage was developed. The AUC of the time-dependent ROC curve and the calibration curve was conducted to weight the 5-year OS prediction capability and the calibration capability of the PS-based nomogram, separately. Data from the external cohort were analysed to further validate this model. Stratification analysis in subgroups was conducted to estimate the discriminative abilities for risk stratification of this model as well.

| Assessment of the therapeutic response
To evaluate the sensitivity of drugs commonly treated for BC patients, we measured the IC50 of five chemotherapeutic agents and one small-molecule inhibitor agent (paclitaxel, doxorubicin, cisplatin, lapatinib, methotrexate and vinorelbine) acquired from the GDSC website and the R package pRRophetic was applied to calculate the IC50 of these agents in the BC patients. Wilcoxon signed-rank test was performed to compare the difference in the IC50 between the high-and low-risk groups.

| Statistical analysis
Statistical analyses were carried out utilizing Stata/MP, version 14.0 (Stata Corp LP, College Station, TX), and R software (www.r-proje ct.org, version 3.6.1). The χ2 test and Fisher's exact test were applied to compare the categorical variables in the primary and external validation cohorts. The Kruskal-Wallis test was utilized to calculate the differences of nomogram score across different subgroups which were stratified by following clinical characteristics such as tumour size, lymph node status, TNM stage and subtype. The Wilcoxon test was performed to compare the stromal score, estimated score, immune score and immune checkpoint expression levels between two risk groups. Pearson correlation tests were employed to finish the correlation analysis of immune infiltration levels. Kaplan-Meier survival analysis and a log-rank test were utilized to draw the survival curve and to compare the difference between groups. The independent prognostic predictors were determined by UMCRA.
Statistically significant was defined as the p value of analyses < 0.05. respectively. The RNA expression levels of these 33 PRGs between normal and tumour tissue were illustrated as a heatmap ( Figure 1A).

| Identification of PRGs between breast normal and tumour tissues
Then, we used the spearman method to analyse the correlation coefficient between any two PRGs ( Figure 1B). We identified that the expression levels of TNF, IL18, NLRP3, NLRP4, NOD1, NLRP1, CASP4, CASP1, CASP5 AIM2, CASP8 and GSDME were positively correlated with at least one of them. Next, we performed a proteinprotein interaction (PPI) analysis to investigate the interactions of the 33 PRGs with the minimum required interaction score setting at 0.9 ( Figure 1C). In addition, the correlation network (red: positive correlations; blue: negative correlations) of part of these PRGs that could link to each other was also presented ( Figure 1D).

| Functional enrichment analysis and mutation information of PRGs
To investigate the biological function and behaviour of the 33 PRGs, we performed GO analysis. The results revealed that these PRGs were chiefly involved in these terms: nucleotide-binding oligomerization domain (NOD) pathway, pyroptosis, inflammasomes, protein processing, apoptotic signalling pathway, necrotic cell death, inflammatory response to an antigenic stimulus, regulation of autophagy and negative regulation of NF-kappa B transcription factor activity (Supplementary Figure 1A

| Construction of prognostic-related PRGs risk score in BC patients
UMCRA were utilized to screen independent prognostic-related PRGs at the base of 33 PRGs. After the analysis, four PRGs were recognized to be significantly correlated with the survival of BC patients, including CASP9, TIRAP, GSDMC and IL18. Then, four PRGs were applied to calculate the pyroptosis score (PS).  Specifically, 138 cases were grouped to cluster 1 and 858 cases were grouped to cluster 2. The 5-year overall survival rate of cluster 2 was significantly poorer than that of cluster 1 (Supplementary Figure 2D).

| Differences in clinical features and ICI between the two clusters
The clinicopathological features consisted of the tumour size (T1-T4), lymph node status (N0-N3), TNM stages (I-IV) and survival status (alive or dead) were presented in a heatmap (Supplementary Figure 2E).
Then, we analysed the distribution differences of ICI between two clusters and found that M1 macrophages, activated CD4 memory T cells, CD8 T cells and follicular helper T cells were enriched in cluster 1 while M0 macrophages, resting T4 memory cells, resting mast cells and M2 macrophages were abundant in cluster 2 (Supplementary Figure 2F). The result indicated that activated innate ICI was distributed much more in cluster 1, thus conferring a better survival.

| Establishment and validation of the PSbased nomogram
A total of 996 BC patients were extracted from the TCGA cohort as the primary cohort, and 3,262 BC patients were acquired from the GEO cohort (GSE96058) as the external validation cohort. No significant differences were detectable in the baseline characteristics between these two cohorts (all p > 0.05; Table 1). Next, UMCRA of clinical variables and PS were conducted and the results are summarized in Table 2. It demonstrated that age, TNM stage, tumour subtype and PS were independent prognostic predictors. Therefore, to further create a clinical tool as applied to predict the 5-year OS, we constructed a PRGs-related nomogram based on the patient's age, TNM stage, tumour subtype and pyroptosis score in the primary cohort and it was validated in the external validation cohort (Figure 2A).

| Clinical evaluation by PSbased nomogram score
By using the optimal nomogram score of each cohort as cut-off values, we separated the BC patients into a low-risk group (N = 640) and a high-risk group (N = 356). The distribution of risk score and the survival status of each case are shown in Figure 3A (primary cohort) and Figure 3B (external validation cohort). As the nomogram score

| The nomogram retains its ability to predict OS in stratification analysis
To better understand the relationship between clinical features and the nomogram risk model in the primary cohort, a stratification analysis was conducted to determine whether the nomogram retained its prognostic prediction ability in diverse subgroups. As shown in and T3), N0 and N1-3 subgroups, respectively. Likewise, the nomogram also retained its prognostic prediction capacity for different molecular subtyping of BC, such as luminal subtype, triple-positive subtype, HER2-positive (hormonal receptor-negative) subtype and triple-negative subtype. As expected, in these molecular subtypes, all patients in the low-risk group had a survival advantage, compared to that in the high-risk group. These data suggested that the nomogram had the potential to predict OS for BC patients.

| Associations between clinical features and the levels of nomogram score
A series of Wilcoxon rank tests were performed to explore the cor-  Figure 5D).

| Gene set enrichment analysis of PRGs in different risk groups
To further investigate the potential molecular mechanisms leading to the differences between low-risk and high-risk groups, the

| The landscape of immune cell infiltration in risk subtypes
CIBERSORT algorithm, including the gene expression levels, was utilized to investigate the differences of immune cell infiltration (ICI) between the two risk groups of BC patients. We found that the significant differences of the 22 immune cell proportion existed between the two groups ( Figure 6A). Then, to investigate the correlation coefficient values between any two types of immune cells, we used the spearman method and found that activated CD4 memory cells were positively related to CD8 T cells (p < 0.001, r = 0.37), while macrophage M0 was negatively correlated with resting CD4 memory T cells (p < 0.001, r = −0.42) ( Figure 6B). In the next step, three indicators related to BC immunogenicity including immune score, estimate score and stromal score were computed by the expression data (ESTIMATE) algorithm in both groups. All three indicators were evidently higher in the low-risk group ( Figure 6C), indicating that low-risk group patients may have a better response to immune therapy. Additionally, we investigated the distinct distributions of ICI between the two risk groups. As displayed in Figure 6D, the low-risk group had an abundance of naive B cells, resting CD4 memory T cells and resting dendritic cells. In contrast, the high-risk group had a notably higher proportion of macrophage M2. These results indicated that the heterogeneity of ICI between high-and low-risk groups may encompass targets for immunotherapy.

| The association between immune checkpoint genes and risk groups
To assess whether different risk groups had distinct sensitivity to immunotherapy, we also explored the expression levels alteration of immune checkpoints. Specifically, the differential expression levels of immune checkpoints were compared between different risk groups in the primary cohort ( Figure 7A-L). The analysis displayed that the low-risk group had higher expression levels of PD-L1, PD-1, CTLA4, HAVCR2, IDO1, CD8, CXCL9, GZMA and GZMB, thus revealing a better response to immunotherapy. However, no significant difference was found in the expression levels of LAG3, CXCL10 and TNF.

| Assessment of the sensitivity of therapeutic agents
Most early BC patients need to accept chemotherapy to reduce the risks of recurrence and metastasis. However, not all of these patients were sensitive to these drugs. Thus, we selected five chemotherapeutic agents and one small-molecule inhibitor agent (paclitaxel, doxorubicin, cisplatin, lapatinib, methotrexate and vinorelbine) that might benefit patients in the high-risk group. Then, the estimated half inhibitory centration (IC50) value of these drugs was compared and the result revealed that patients in the high-risk group were associated with a higher IC50 of paclitaxel, doxorubicin, cisplatin, methotrexate and vinorelbine, whereas it was associated with lower IC50 for lapatinib ( Figure 7M-R). The results implied that this model might serve as a predictor for chemotherapy responsiveness.

| DISCUSS ION
Pyroptosis has a two-way effect in cancers. Stimulation by pyroptosis releasing inflammatory factors can transform normal cells into tumour cells. 33 In contrast, pyroptosis can facilitate cancer cell death, thus play an antitumour role. 34 In BC, different expression levels of pyroptosis-related cytokines, including inflammasome, gastermin F I G U R E 2 Nomogram for predicting the 5-year OS of patients in the primary cohort (A). ROC curve and calibration plot for validating the accuracy of the nomogram in the primary cohort (B, D) and the external validation cohort (C, E) family and caspase may lead to exactly opposite effects. However, the expression levels of several gastermins did not appear to act as reliable prognostic factors. The role of PRG signature in BC has not yet been elucidated. Therefore, a prognostic signature by combining independent prognostic PRGs in BC needed to explore.
In our study, the expression patterns, biological functions and genetic alterations of 33 PRGs were illustrated. As expected, functional analyses showed that these PRGs were chiefly enriched in cell death-related pathways such as pyroptosis, necrotic cell death, apoptotic signalling pathway and regulation of autophagy.
Additionally, they were also enriched in the cancer-related NF-κB pathway. However, only four (CASP9, TIRAP, GSDMC and IL18) independent prognostic predictors were finally confirmed for further analysis. Among these four PRGs, CASP9 might act as antioncogenes while TIRAP, IL-18 and GSDMC might be oncogenes.
Upregulation of CASP9 by miR-182-5p inhibition can induce apoptosis and antiproliferative effect in breast cancer. 35 Similarly, upregulation of CASP9 through NF-κB contributed to apoptosis activation in cancer cells, thus exerting the antitumour effect. 36 On the contrary, downregulation of CASP9 by microRNA-224 37 or inactivation of CASP9 by lncRNA POU3F3 38 promoted proliferation and progression in TNBC. TIRAP was reported to be an upstream regulator of the NF-κB pathway and its expression enhanced cell proliferation in lymphoma. 39 Another study revealed that downregulation of TIRAP contributed to the anti-proliferation effect in non-small-cell lung cancer cells. 40 Studies about GSDMC suggested that it functioned as an oncogene, facilitating cell proliferation in colorectal cancer 41 and serving as an independent prognostic predictor in lung adenocarcinoma. 42 Besides, GSDMC was found to be highly expressed in BC and correlated with poorer survival. 43 A previous study reported that tumour-derived IL18 could induce PD-1 expression on NK cells, which resulted in a poor prognosis in patients with TNBC. 44 Under the stimulation of leptin, IL18 expression increased and promoted BC cell migration and invasion. 45 Another study reported a similar conclusion that high serum levels of IL18 resulted in worse outcomes in patients with BC. 46 These results showed that different expres-

| CON CLUS ION
In summary, our study comprehensively evaluated the prognostic value, functions in TIME regulation, and potential therapeutic

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets for this study can be found in The Cancer Genome Atlas (TCGA) data set (https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO) databank (https://www.ncbi.nlm.nih. gov/geo/).