A novel pyroptosis-related indicator of immune infiltration features and prognosis in breast cancer.

Breast cancer is the most common malignancy in women, and there is evidence for the dual role of cell pyroptosis in tumor development. However, little is known about the relationship between cell pyroptosis and breast cancer and its prognostic value. We aimed to construct a prognostic model using cell-pyroptosis-related genes to provide innovative insights into the prognosis and treatment of breast cancer. We screened candidate genes for pyroptosis using public databases and identified 10 cell pyroptosis signature genes with the random forest method. Finally, a nomogram for predicting 1-, 3-, and 5-year survival probabilities was constructed. The differences in immune cell distributions between survival periods were similar across the breast cancer datasets. The 10 identified key pyroptosis factors showed a significant correlation with Her2, tumor-node-metastasis (TNM) stage, and survival of breast cancer. The risk scores correlated positively with the infiltration features of naive B cells, CD8+ T cells, atpdelnd mast cells, while they correlated negatively with those of M0 macrophages and dendritic cells. In conclusion, our findings confirm that cell pyroptosis is closely associated with breast cancer. Importantly, the prognostic complex values generated from the 10 cell-pyroptosis-related genes based on various clinical features may provide an important basis for future studies on the prognosis of breast cancer.

Introduction

According to the latest data on the global cancer burden, in 2020, breast cancer had the highest morbidity worldwide and was the leading cause of cancer-related deaths, with approximately 2.26 million newly diagnosed cases and 685,000 reported deaths (1). Among women, breast cancer accounts for 24.5% of all cases and 15.5% of cancer-related deaths; it ranks first among cancers in terms of incidence and cancer-related mortality in most countries (2). Owing to the continuous progress in medical treatment, the survival rate of patients with breast cancer has improved greatly, and the current 1-, 3-, and 5-year survival rates are 0.92, 0.75, and 0.73, respectively. Although the survival rate associated with breast cancer is better than that for most other tumors, the 10-year survival rate remains low (0.61), and some patients experience distant recurrence at different times, with an average survival time of 2 years after recurrence (3). Given this situation, a reliable prognostic model may help improve the survival rate of patients.

Cell pyroptosis, also known as cellular inflammatory necrosis, is an inflammatory form of programmed cell death (PCD), but it differs considerably from other types of PCD (4). It was first identified in infected macrophages in 1992 (5), and the term was coined by Cookson et al. (6). PCD induced under various stimuli has been studied extensively in many disease models, and cell pyroptosis was found to be mediated by gasdermin (GSDM) (7). The GSDM family includes gasdermin A (GSDMA), gasdermin B (GSDMB), gasdermin C (GSDMC), gasdermin D (GSDMD), gasdermin E (GSDME), and PJVK (8). Except for PJVK, the other members of the GSDM family have both C- and N-terminal structural domains. GSDM is cleaved, and the GSDM-N structural domain is released (9). The released N-terminal domain perforates the cell membrane and leads to characteristic morphological changes associated with pyroptosis, including cytoplasmic swelling, membrane rupture, and release of inflammatory factors into the extracellular environment, thus directly amplifying the systemic immune responses (9). The main manifestation of pyroptosis in cells is swelling that results in cell membrane rupture, which leads to the efflux of contents and activation of an intense inflammatory response (10).

The role of cell pyroptosis in cancers is gradually becoming evident. In esophageal squamous carcinoma, head and neck squamous carcinoma, and hepatocellular carcinoma, pyroptosis can be induced by different drugs to achieve tumor suppression (11–13). On the other hand, pyroptosis inhibition protects microglia and neurons, rescues dopaminergic neurons, inhibits neuroinflammation, and alleviates neurodegeneration (14). Therefore, pyroptosis is a double-edged sword that plays a key role in antitumor immunity in certain tumors and may provide an effective treatment strategy for cancer; however, its induction in normal tissues and immune cells can lead to severe damage (15). Thus, in the case of tumors, on the one hand, pyroptosis causes inflammation that contributes to the generation and maintenance of an inflammatory microenvironment around cancer cells, thus promoting tumor progression; on the other hand, acute activation of pyroptosis leads to the infiltration of multiple immune cell types that inhibit tumor progression (16). In non-small cell lung cancer, GSDME-mediated pyroptosis can make immune factors MIP-1α, MIP-1β, MIP-2, and IP-10 increase, thereby recruiting T cells to achieve anti-tumor effect (17). Some chemotherapeutic drugs can induce the expression of GSDMC and activate caspase-8 and induce the pyroptosis of breast cancer cells (18). Pyroptosis also has different anti-tumor mechanisms in the study of hepatocellular carcinoma, gastric cancer, and ovarian cancer (19). To examine the prognostic impact of pyroptosis, a pyroptosis score was constructed for melanomas. The results indicated that cell pyroptosis is an independent prognostic factor that can improve survival by strong immune clearance in tumors through immune cells, including T cells, B cells, and natural killer (NK) cells (20). Furthermore, a recent study found that N-oxide induced pyroptosis in tumor cells by activating the endoplasmic reticulum stress kinase PERK and, thus, enhanced CD8+ T-cell-mediated antitumor immunity in triple-negative breast cancer (TNBC) in the in vivo setting (21). All these findings reveal the close relationship between pyroptosis and breast cancer. Specifically, cellular pyroptosis plays an important role in tumorigenesis and antitumor processes, but its specific function in breast cancer is unclear.

This systematic study was conducted to determine the expression of genes associated with cell pyroptosis in normal breast tissues and breast cancer tissues in order to evaluate their prognostic value and further investigate the correlation between cell pyroptosis and the tumor immune microenvironment. We aimed to comprehensively assess the relationship between cell pyroptosis and the prognosis of breast cancer patients and also constructed a nomogram based on breast cancer pyroptosis-related genes. The efficacy of this nomogram was evaluated for prognostic prediction, molecular characterization, clinical significance, and assessment of the regulatory function of the immune microenvironment.

Material and methods

Data collection

The Cancer Genome Atlas Breast Cancer (TCGA-BRCA) expression profile, variants, clinical information, and follow-up data were downloaded from XENA, and samples with complete phenotypic and survival data were retained for further analyses. The GSE96058 expression data and corresponding sample information were downloaded from the GEO database. The METABRIC expression data and corresponding sample information were downloaded from the ciBoportal database. The data were preprocessed as follows: first, probes corresponding to genes according to the annotation file were included, while unannotated probes were removed. If there were multiple probes that corresponded to the same gene, the probe with the maximum value was considered to represent the level of expression of that gene. All genes with low expression were filtered out at a cutoff of >1 in at least 10% of the samples.

Differential gene analysis

Differential analysis of data normalized to fragments per kilobase of transcript per million (FPKM) in TCGA-BRCA was performed using the limma package in R. For screening, the criteria were set as |FC| > 1.5 and p < 0.05.

Cox regression and survival analyses

For Cox univariate analysis, regression modeling of individual genes or clinical features [age, tumor–node–metastasis (TNM) stage, LumA, LumB, TNBC, and Her2 status] was performed using the coxph function in the survival package. Prognosis-related genes or clinical features were screened at a cutoff p-value of <0.05. After the corresponding modeling parameters were extracted, forest plots were drawn using the forestplot package.

For survival analysis, related genes were screened and grouped according to their median expression level. Overall survival (OS) and grouping information were fitted using the survfit function in the surv package and finally analyzed and visualized using the ggsurvplot function in the survminer package.

Screening of pyroptosis factors and gene set variation analysis by the random forest method

Potential pyroptosis-related factors were screened using the random forest (RF) algorithm. A three-time and 10-fold crossover model was constructed using the carat package. The optimized parameters were input into the RF package, and finally, genes with the top 10 MeanDecreaseGini scores were selected. The pROC package was used for prognosis prediction of the tumor samples, and receiver operating characteristic (ROC) curves were drawn. Finally, the genes with the top 10 MeanDecreaseGini scores were selected as the key pyroptosis factors.

The pyroptosis signature gene expression values were determined from normalized expression data in TCGA-BRCA, and gene set variation analysis (GSVA) was performed on the normalized expression data using the GSVA package to obtain the complex values.

Prediction of survival using the prognostic complex

Regression modeling of the complex values was performed using the coxph function in the survival package, and the samples were classified into high- and low-risk groups according to the median of the complex values. Finally, OS was analyzed.

Validation of the prognostic complex and nomogram construction

In order to verify the independent prognostic efficacy of the high and low complex values, univariate Cox analysis was first performed on the TCGA-BRCA dataset in combination with other clinicopathological characteristics, including age, TNM stage, LumA, LumB, TNBC, and Her2 status. Next, overall prognosis based on the above factors was analyzed by multi-factor Cox regression analysis to verify the independent prognostic efficacy of the risk score. The Cox proportional hazards regression model was constructed using the cph function in the R package rms. Finally, the survival package was used to calculate survival probability, and the nomogram function was used to construct the nomogram and plot the calibration curves to evaluate its accuracy and predict its utility.

Immune infiltration analysis

The immune infiltration score files for TCGA were downloaded from the TIMER2 (http://timer.comp-genomics.org/) database, and TCGA-BRCA sample data were screened. The corresponding CIBERSORT score data were used for comparison of immune cell differences among samples with different survival values, using the t_test function in the rstatix package in R. Difference analysis and correlation analysis between the complex value and immune cell proportions were conducted by the Pearson method. Immune infiltration scoring was performed using TIMER2 along with normalized expression data from GSE96058 and METABRIC. Finally, immune cell differential analysis and mapping were performed.

Enrichment analysis and correlation analysis of key genes associated with the clinical features

The enrichment of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was analyzed using the clusterProfiler package. The enriched GO terms were further analyzed by GSEA. For GO terms that were enriched in both steps, similarity calculation and hclust clustering were performed using the GOSemSim package.

For screening tumor samples for TCGA molecular typing and stage, mosaic correlation analysis was performed on survival-related clinical information, and the expression levels of pyroptosis-related key genes were analyzed with the vcd package.

Mutation analysis and molecular docking

The Mutect2 mutation files of TCGA-BRCA were downloaded using the R package, TCGAbiolinks, and this was followed by visualization of mutation types using the R package maftools. The corresponding compound structures were downloaded from the DrugBank database (https://go.drugbank.com/) and screened according to Lipinski’s rules (hydrogen bond acceptor ≤ 10, hydrogen bond donor ≤ 5, rotatable bond ≤ 10, the logarithmic value of lipid–water partition coefficient ≤ 5, molecular weight 180–480, and polar surface area ≤ 140). Finally, 5,464 small-molecule compounds were obtained. The 3D structural information of the proteins encoded by the characteristic pyroptosis-related genes in breast cancer was queried in PDB (https://www1.rcsb.org/). Relevant structural information was obtained for PAK7, and its corresponding PDB file, 2F57, was downloaded. The approximate docking box range was identified according to the ligand information therein, and other relevant parameters of AutoDock Vina were set. Docking with small-molecule compounds was performed with AutoDock Vina, and demonstration and analysis of interaction forces were performed using Pymol and Ligplus.

Drug sensitivity analysis

We analyzed the drug resistance of genes based on the GDSC drug database. The mRNA expression data and drug sensitivity data were merged. Pearson correlation analysis was performed to determine the correlation between mRNA expression of the screened pyroptosis-related genes and the IC50 values of drugs. The p-value was adjusted by false discovery rate (FDR).

Codes

The codes used to analyze the data can be found here: https://github.com/wzhlc1206/ A-Novel-Pyroptosis-Related-Signature-Indicates-Immune-infiltration-Features-and-Predicts-Prognosis-.git

Results

Identification of pyroptosis-related genes

First, we combined the pyroptosis-related genes from the PMID:338828074, MsigDB, and GeneCards databases ( ) and obtained a pyroptosis gene set containing 57 genes ( ). Next, the correlation of these genes with the normalized expression data in the TCGA-BRCA dataset was analyzed. Genes with Pearson correlation coefficients >0.6 and p-values <1e−10 were selected as candidate pyroptosis-related genes. Finally, 4,939 genes and 4,996 pyroptosis-related genes were obtained. The heatmap of partial genes is shown in , and the RNA expression profile for breast cancer is shown in .

Figure 1

Screening of candidate genes for pyrosis. (A) Heatmap of the correlation between partial differences in the TCGA-BRCA dataset and pyroptosis genes. (B) Volcano plot depicting the results of differential analysis between tumor and normal datasets and (C) survival time. (D) Venn diagram of intersecting pyroptosis genes. *P < 0.05, **P < 0.01, ***P < 0.001.

FPKM expression data from TCGA-BRCA were processed using limma and filtered to exclude genes with low expression levels. These values were normalized for principal component analysis (PCA). Short-term survival was set as <5 years, while long-term survival was set as more than 5 years. By comparing the expression of differentially expressed genes, we found that there were significant differences between normal and tumor samples ( ). In contrast, differences in short- and long-term survival were less significant ( ).

Differential analysis of TCGA-BRCA data was performed using the limma package according to the presence or absence of disease. Differential genes were screened between breast cancer and normal tissue samples ( ). The genes that were associated with survival are listed in . The screening criteria were as follows: |FC| > 1.5 and p < 0.05. A total of 10,913 genes were upregulated, while 10,577 were downregulated in breast cancer samples as compared to the normal samples ( ). A total of 2,255 intersecting genes related to pyroptosis were obtained ( ). The breast cancer samples were grouped according to survival time for differential analysis: 994 genes were upregulated, while 906 were downregulated in patients with long-term survival ( ). The last two sets of intersecting genes were overlapped, and 29 candidate genes for pyroptosis were identified ( ).

Screening of core genes associated with pyroptosis

Tumor samples in the TCGA-BRCA dataset were screened, and univariate Cox regression analysis was performed for 29 pyroptosis-related genes using the survival package. Based on the median gene expression values, the genes were classified into high- or low-expression groups. The results showed that a total of 19 genes were associated with breast cancer prognosis ( , p < 0.05). Survival analysis of these prognosis-related pyroptosis genes based on the median threshold showed that 14 potential pyroptosis-related genes were significantly associated with survival in breast cancer: IGLV1-44, IGLV1-51, IGLV1-66, CALML3, IGHA2, PAK7, TNN, TP53AIP1, PCP2, TRBV7-4, EDAR, IGLV3-19, FREM1, and ZMYND10 ( ). High expression of these genes for up to 6,000 days was associated with better survival, and these results were highly significant for TP53AIP1.

Figure 2

Genes related to the prognosis of breast cancer. (A) Results of univariate Cox analysis of the TCGA-BRCA dataset. (B–O) Prognosis-related pyroptosis gene expression and breast cancer survival curve.

To identify the core genes associated with pyroptosis, we further screened these 14 potential pyroptosis genes. A three-time 10-fold crossover model was constructed using the carat package, and finally, the normalized dataset was analyzed using the RF method with optimized parameters. Gini scoring was performed, and finally, the top 10 genes according to their MeanDecreaseGini scores were selected as the core genes associated with pyroptosis. These genes included PAK7, TP53AIP1, PCP2, ZMYND10, IGHA2, IGLV1-44, TNN, FREM1, IGLV1-51, and TRBV7-4. The potential pyroptosis-related genes according to their MeanDecreaseGini scores are shown in . The heat map for their expression is shown in , and show the expression of these genes, along with survival, in breast cancer patients.

Figure 3

Identification of core pyroptosis-related genes and relationship between the expression of these genes and survival. (A) MeanDecreaseGini score after RF modeling of potential scores (light blue indicates the top 10 pyroptosis factor scores). (B) Heat map of the association between the level of expression of the top 10 pyroptosis factors and breast cancer survival. (C–L) Survival analysis of the pyroptosis factors.

Prediction of OS based on prognostic complex value and validation of independent prognostic factors

The pyroptosis-related genes identified after screening were used as the defined set for GSVA. GSVA was performed for breast cancer samples using the GSVA package, and the complex prognostic value was obtained. The samples were then divided into high- or low-risk groups according to the median complex prognostic values. The correlation between survival and pyroptosis-related gene expression was determined. Significant differences in survival and LumB, LumA, and Her2 subtypes were found between the high- and low-risk groups. Association analysis of survival and expression was performed for the overall sample and five subtypes according to LumB, Basal, LumA, Normal, and Her2 status ( ), based on the complex prognostic values calculated using GSVA for breast cancer samples.

To verify the independent prognostic efficacy of the high and low complex prognostic values, Cox univariate analysis was first performed on the TCGA-BRCA dataset in combination with other prognostic factors, including stage, age, Her2, LumA, LumB, TNBC status, and TNM stage ( ). Overall prognosis based on the above six factors (including the complex prognostic values) was then analyzed by multifactor Cox regression ( ). The results confirmed that the complex values were independent prognostic factors and could be used to develop prognostic column line plots that could predict the probability of survival at 1, 3, and 5 years. Nomograms of 1-, 3-, and 5-year survival probabilities were verified by the calibration graph method, and the standard curve of the calibration graph was in good agreement with the calibration prediction curve ( ). The predicted and observed values for 1-, 3-, and 5-year survival were in good agreement.

Figure 4

Independent prognostic efficacy of complex prognostic values. Results of (A) univariate and (B) multivariate Cox regression analyses, (C) nomograms, and calibration curves for survival at (D) one (E) 3 and (F) 5 years.

Correlation analysis of complex prognostic values and immune-infiltrating cells

The immune infiltration score files for TCGA were downloaded from the TIMER2 database, and data related to the breast cancer samples in TCGA-BRCA were selected and grouped according to survival (with 5 years set as the threshold). The CIBERSORT score data were used for immune cell differential analysis and mapping ( ). Immune infiltration scoring was performed using TIMER2 for standardized expression data from GSE96085 and METABRIC, followed by immune cell differential analysis and mapping. The results of immune infiltration analysis showed that the complex values were positively correlated with infiltration of naive B cells, CD8+ T cells, and mast cells, while they were negatively correlated with M0 macrophages and dendritic cells ( ).

Figure 5

Correlation of complex prognostic values and immune infiltrating cells. (A) Differences in the distribution and proportion of immune cells in the TCGA-BRCA, GSE96058, and METABRIC datasets. (B–F) Correlation analysis of complex values with the proportion of different immune cell types.

Enrichment analyses of key functional genes and pathways

The clusterProfiler package was used to perform GO-BP ( ) and KEGG pathway enrichment analyses ( ) ( ) for the top 10 pyroptosis factors identified with the RF method and enriched pathways from GSEA ( ) ( ). For BPs enriched in the first two steps, similarity calculation and hclust clustering were performed using the GOSemSim package. Immune responses, defense responses, and cell recognition were enriched ( ); this indicates that these genes were involved in regulating similar BPs.

Figure 6

Enrichment analyses of key functional genes and pathways. (A) Bubble plot for enriched GO-BP and KEGG pathways related to key pyroptosis factors. (B) GSEA for enriched pathways in the TCGA-BRCA samples. (C) The enriched pathways in panels (A, B) are classified into three major categories of biological functions based on similarity.

Mutation analysis of key genes and their correlation with clinical features

The TCGA-BRCA point mutation maf files of Mutect analysis were downloaded from TCGA. The mutation data of tumor samples and pyroptosis core genes were screened, and this was followed by visualization of mutation types using the R package mafTools. Mutation analysis showed that FREM1, TNN, and PAK7 were the top 3 key pyroptosis genes with the highest mutation frequency ( ).

Figure 7

Mutation analysis of key genes and their correlation with clinical features. (A) Waterfall plot for point mutation analysis of key pyroptosis factors in samples from the TCGA-BRCA dataset. (B) Statistical plot of point mutation types. (C, D) Point mutation distribution for FREM1 and PAK7 and (E–N) the relationship between key pyroptosis factors and TCGA molecular typing, Her2 status, and stage of breast cancer.

Screening of clinical information related to TCGA molecular typing and mosaic correlation analysis of Her2 status, stage, survival, and levels of expression of the 10 key genes of pyroptosis in breast cancer samples using the vcd package showed that there was a significant correlation between clinical features and expression of the 10 pyroptosis genes ( ).

Screening of potential therapeutic compounds based on molecular docking experiments

Protein structure information was downloaded from the PDB database. AutoDock Vina was used for molecular docking experiments. Compounds were filtered based on their binding energies, and finally, Pymol was used for drawing protein–compound binding maps. The corresponding compound structures were downloaded from the DrugBank database (https://go.drugbank.com/) and filtered according to Lipinski’s rule. Finally, 5,464 small-molecule compounds were obtained. 3D structural information about the proteins encoded by characteristic pyroptosis-related genes in breast cancer was obtained from the PDB database (https://www1.rcsb.org/), and only relevant structural information for PAK7 was obtained (the other genes did not have corresponding structural files containing ligands based on which the docking box range could be inferred). The corresponding PDB file, 2F57, was downloaded. The approximate docking box range was calculated based on the ligand information therein, and other relevant parameters were set for AutoDock Vina. AutoDock Vina was used for docking small-molecule compounds, and Pymol and Ligplus were used for demonstration and analysis of the interaction forces. The top 3 small-molecule compounds with the strongest binding abilities were DB08676, DB08674, and DB06888 ( ). Docking conformation and interaction force analysis of pak7 and these three small-molecule compounds are shown in . The complete docking scores are provided in .

Table 1

Top 10 compounds with the best docking scores with PAK7.

DrugBank ID	Hydrogen	Hydrogen	Rotatable	LogP	Molecular	TPSA	Affinity
	Acceptors	Donors	Bonds		Weight		(kcal/mol)
DB08676	5	1	0	2.5	453.5	83.2	-11.1
DB08674	5	1	0	2.8	435.5	83.2	-10.9
DB06888	6	2	0	2.2	403.4	86.3	-10.7
DB15054	5	2	3	2.5	424.5	86.2	-10.7
DB12012	8	1	4	3.3	455.4	80.2	-10.5
DB12200	2	2	2	2.9	369.4	66.9	-10.4
DB12523	5	1	5	3.2	432.5	65.8	-10.2
DB06469	4	3	1	2.2	439.5	88.6	-10.1
DB08683	3	1	0	3.8	393.4	65.3	-10.1
DB12611	5	4	4	1.5	419.5	130	-10.1

Figure 8

Docking conformation and interaction force analysis for PAK7 and three small-molecule compounds. (A) DB08676, (B) DB08674, and (C) DB06888. Top panel: Pymol shows the docking conformation and hydrogen bonding, wherein cyan is the small molecule, the yellow dashed line is the hydrogen bond, and blue is the amino acid residue forming the hydrogen bond with the small molecule. Bottom panel: Ligplus force analysis, wherein the small molecule is in the middle, surrounded by related protein amino acid residues; the green dashed line is the hydrogen bond; and the amino acid in green is the amino acid residue forming a hydrogen bond.

To explore the relationship between pyroptosis genes and drug sensitivity, we analyzed the drug sensitivity of genes based on the GDSC drug database. The results showed that ZMYND10, PCP2, TNN, and FREM1 were associated with drug response. Patients with cancer with high expression of the PCP2 and TNN genes were probably sensitive to AS605204 and FK886 and resistant to XAV939. Furthermore, cancers with high expression of the FREM1 gene were sensitive to XAV939 and dasatinib and resistant to SB590885 and gemcitabine, and cancers with high expression of ZMYND10 were resistant to XAV939, gemcitabine, and SB590885 ( ).

Figure 9

Drug sensitivity analysis of the pyroptosis genes. Data from the GSCA database represent the correlation between the expression of proptosis genes and drug sensitivity. A negative correlation indicates that patients with cancer with overexpression of the proptosis genes are sensitive to the drug, and vice versa.

Discussion

According to the latest statistics, breast cancer ranks first among tumors in women in terms of incidence (22), and the number of cases is expected to reach 4.4 million in 2070 (23). Increasing evidence shows that cell pyroptosis is closely associated with the development of lung, gastric, cervical, and breast cancers (24). Some scholars have used long-chain non-coding RNAs or lncRNAs associated with pyroptosis to predict the prognosis and tumor immune microenvironment of breast cancer patients (25); however, the role of pyroptosis-related genes in breast cancer remains unclear. Therefore, there is a need to examine the prognostic profiles of pyroptosis-related genes by combining independent prognostic factors associated with breast cancer. Herein, we describe the expression, mutation, immune infiltration, and clinical significance of pyroptosis-related genes in breast cancer.

In our study, most pyroptosis-related genes were differentially expressed in breast cancer, and this indicated their potential roles in tumorigenesis (26). We constructed a prognostic gene signature based on 10 cell pyroptosis genes (PAK7, TP53AIP1, PCP2, ZMYND10, IGHA2, IGLV1-44, TNN, FREM1, IGLV1-51, and TRBV7-4) and generated a nomogram for predicting the probability of survival at 1 and 3 years. Survival analysis showed significant differences in OS between high- and low-risk groups classified according to the complex prognostic values and between groups with high and low complex prognostic values for three breast cancer subtypes based on LumB, LumA, and Her2 status. The results of univariate and multifactorial Cox analyses revealed the independent prognostic effects of the complex prognostic values. Furthermore, the differences in immune cell distributions according to survival were similar across the different breast cancer datasets.

In the present study, the 10 key pyroptosis-related factors identified were significantly correlated with TCGA molecular typing based on the Her2 status and stage of breast cancer samples. Some scholars have pointed out that GSDMB overexpression/amplification occurs in 60% of HER2-positive breast cancers (27). In addition, higher expression of DRD2 was positively correlated with longer survival times, especially in HER2-positive patients, and DRD2 also triggered GSDME-mediated pyroptosis (28). This implies that HER2 may be more closely related to pyroptosis. Mutation analysis showed that FREM1, TNN, and PAK7 were the top 3 key pyroptosis genes with the highest mutation frequencies. Only structural data for PAK7 were available in the PDB database, and the results of molecular docking revealed that PAK7 had the strongest binding ability for the three small-molecule compounds examined, namely, DB08676, DB08674, and DB06888. Some studies have shown that small molecular drugs that can regulate immune cells will greatly improve the efficacy of cancer immunotherapy (29). Small molecular compounds can change the molecular pathway to deal with immune tolerance and overcome tumor-related immunosuppression, so as to produce effective anti-tumor response and improve the efficacy of cancer immunotherapy (30). Small molecular drugs are compatible with systemic administration and are suitable for extracellular and intracellular targets (31). In future research, we will pay more attention on how to improve the efficacy of breast cancer by combining immunotherapy with DB08676, DB08674, and DB06888 small molecular compounds. PAK7, an oncogene, is significantly upregulated in both gastric and esophageal squamous carcinoma (32, 33). Its expression is significantly enhanced in breast cancer tissues and correlates positively with its pathological differentiation and TNM stage. PAK7 is involved in breast cancer progression through the activation of the Wnt/β-linked protein signaling pathway; this indicates the potential applicability of PAK7 as a therapeutic target for the treatment of breast cancer (34). In fact, one study has reported that inhibition of osteosarcoma can be achieved by targeting PAK7 (35). Studies on non-small cell lung cancer suggest that the PAK7 mutant phenotype is also associated with the tumor immune microenvironment (36). Similarly, in our study, it was confirmed that PAK7 was correlated significantly with the TCGA molecular type, Her2 stage, TNM stage, and the immune microenvironment in breast cancer. Notably, TP53AIP1, ZMYND10, and FREM1 are associated with the apoptotic pathway (37). Moreover, apoptosis and pyroptosis are tightly linked to each other through multiple mechanisms (38). ZMYND10 suppresses breast cancer oncogenicity by inhibiting the miR145-5p/NEDD9 signaling pathway (39). As an oncogene, ZMYND10 promotes apoptosis in tumor cells by regulating the activity of sMEK1 (40) and inhibiting angiogenesis. In addition, ZMYND10 can upregulate BAX expression, which facilitates the pro-apoptotic pathway, thereby promoting paclitaxel-induced apoptosis in ovarian cancer cells. Reduced FREM1 expression is usually associated with hormone-receptor-negative and TNBC status, and correlates significantly with poor OS and recurrence-free survival. Furthermore, FREM1 expression correlates positively with the level of immune infiltration of CD4+, CD8+ T cells, and CD86+ M1 macrophages, while it correlates negatively with CD68+ pan− and CD163+ M2 macrophages. These findings suggest that FREM1 is a potential biomarker for assessing immune infiltration status and BC prognosis (41, 42). Interestingly, TILRR enhances IL-1-induced anti-apoptotic responses (43) and is a splice variant of the FRAS1 family (44). It has also been shown that high IGHA2 mRNA expression is associated with a more favorable prognosis (45); this is consistent with our findings. IGHA2 is significantly associated with the OS of patients with esophageal squamous cell carcinoma, as evaluated using a prognostic risk score model of immune-related genes (46). The close relationship between apoptosis and cell pyroptosis may explain the dual roles of TP53AIP1, ZMYND10, and FREM1. IGLV1-44 and IGLV1-51 are immune-related genes (47–49) that can influence immune infiltration. Our results suggest that several pathways associated with the pyroptosis-related genes were enriched. Thus far, the relationship of TNN, TRBV7-4, and PCP2 with cancer has not been reported, and our results suggest that these genes, as core pyroptosis genes, have important associations with breast cancer development and the immune microenvironment.

Breast cancer is often considered a cold tumor with reduced immune cell infiltration, suppressed immune microenvironment, and low mutation frequency (50). Importantly, immunotherapy may produce lasting therapeutic effects. Therefore, many scholars have studied and explored immunotherapy strategies for breast cancer. For example, atezozumab is a monoclonal antibody against PD-L1. In a phase III study, atezozumab combined with paclitaxel was found to affect the progress of TNBC and improve its prognosis. Accordingly, treatment with atezozumab has been approved by the Food and Drug Administration (FDA) for PD-L1-positive advanced or metastatic TNBC (51). Other immunotherapies for breast cancer, including cytotoxic T-lymphocyte-associated antigen 4 or CTLA-4, tumor-infiltrating lymphocytes, and tumor vaccines, have also achieved relatively good results in the research setting and may pave the way for more breast cancer treatments that can be applied in the clinical setting (52). The combination of cell pyroptosis induction and immunotherapy may synergistically increase anticancer activity (53). Accordingly, it has been indicated that cell-pyroptosis-related gene signatures correlate with almost all steps of the cancer immune cycle (54). The results of immune infiltration analysis in the present study showed that the complex prognostic values correlated positively with the infiltration of naive B cells, CD8+ T cells, and mast cells, while they correlated negatively with M0 macrophages and dendritic cells. Thus, cell pyroptosis pathways may have an impact on the immune-oncology landscape. Finally, the data obtained from the GDSC database indicated that patients with increased expression of the ZMYND10, PCP2, TNN, and FREM1 genes may not respond well to treatment with XAV939 and dastatinib; however, they may respond to treatment with AS605240, FK886, etc.

Some limitations of this study need to be mentioned. As the data were obtained from public databases, they need to be validated with laboratory/real clinical data. Yu et al. have identified different pyroptosis genes from us in humans through the TCGA and GEO databases (55), but the good thing is that those pyroptosis genes found by both of us can predict the prognosis of breast cancer. Therefore, in the next study, we will analyze the reported data and verify the relationship between pyroptosis genes and breast cancer through in vivo and in vitro experiments. However, the immune microenvironment has a very important impact on tumors, so we will also focus on the regulation of pyroptosis genes in the immune microenvironment and the effects on related immune cells. In conclusion, our results suggest that induction of cell pyroptosis may be a novel strategy for breast cancer immunotherapy and has broad clinical applicability. The findings may provide an important basis for future studies in terms of developing efficient treatment strategies for breast cancer.

Data availability statement

Publicly available datasets were analyzed in this study. The data can be found here: GSE96058 dataset at https://www.ncbi.nlm.nih.gov/geo/, TCGA dataset at https://xenabrowser.net/datapages/, METABRIC dataset at http://www.cbioportal.org/, and GDSC drug data at https://www.cancerrxgene.org.

Ethics statement

Ethical review and approval were not required for this study, in accordance with the local legislation and institutional requirements. Written informed consent for participation was also not required in accordance with the national legislation and the institutional requirements. Therefore, written informed consent was not obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

CW, LZ, LR, NS, YZ, and XQ conceived and designed the analysis and implemented the experimental studies. CW performed statistical analysis, interpreted the results, graphed the data, and wrote the paper. GZ, AW, SX , HT, ZP, TZ, and PG modified the draft. NS, YZ, and XQ approved the draft of the paper. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by grants from Army Medical University (No. XZ-2019-505-042), Chongqing Basic and Frontier Research Exploration Project (No. cstc2018jcyjA0137), and Chongqing Talent Program Project (No. CQYC20210303406).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.