Exploring the Pyroptosis-related Genes of Breast Cancer and their Effects on Tumor Immune Microenvironment


 Background: Breast cancer （BC） is the most common cancer among women, with high rates of metastasis and recurrence. Some studies have confirmed that pyroptosis is an immune-related programmed cell death. However, the correlation between the expression of pyroptosis-related genes in BC and its prognosis remains unclear. Methods: In this study, we identified 38 pyroptosis-related genes that were differentially expressed between BC and normal tissues. The prognostic value of each pyroptosis-related gene was evaluated using patient data from The Cancer Genome Atlas (TCGA). The Cox regression method was performed to establish a prognostic model for 16-gene signature, classifying all BC patients in the TCGA database into a low-or high-risk group. Results: The survival rate of BC patients in the high-risk group was significantly lower than that in the low-risk group (P＜0.01). Prognostic model is independent prognostic factor for BC patients compared to clinical features. Single sample gene set enrichment analysis (ssGSEA) showed a decrease for immune cells and immune function in the high-risk group. Conclusions: Pyroptosis-related genes influence the tumor immune microenvironment and can predict the prognosis of BC.


Introduction
Breast cancer (BC) is the most common cancer among women aged 20 to 59, accounting for nearly a third of new cancers diagnosed in women [1]. It is a hereditary and multifactorial disease [2,3]. Worldwide, one million women are diagnosed with breast cancer, and about half of them die from the disease each year [4]. By 2020, approximately 26% of BC patients progressed to invasive BC, which may be associated with high mortality rate [5]. Despite advances in diagnosis and treatment, nearly 12% of breast cancer patients end up with metastasis or recurrence, with a poor prognosis [6]. At present, the main treatment of breast cancer remains surgery, combined with radiotherapy, chemotherapy and hormonal treatment when necessary. Immune treatment was also tried which has not bene ted most patients, due to the high metastatic potential and lacking of comprehensive understanding for immune microenvironment of the tumor. Therefore, new therapeutic targets need to be discovered, and a new prognostic model that can assess patient risk and guide targeted therapy is urgently needed.
Pyroptosis is an in ammatory form of programmed cell death (PCD) characterized by cellular swelling and rupture, lysis, nuclear condensation, DNA fragmentation and IL-1β and IL-18 leakage [7], exacerbating the in ammatory response in the extracellular space [8]. It has been con rmed to be associated with in ammatory diseases and malignancies in humans. Gasdermins are a family of pore-forming proteins. It was found that the pore forming characteristics of n-terminal segments of Gasdermin D (GSDMD) and Gasdermin E (GSDME) were the driving factors of pyroptosis [9][10][11], these two kinds of proteins are cleaved by activated Caspase-1 and Caspase-3, respectively, leading to pyroptosis [9,12]. Further studies showed that other members of the Gasdermin family, including GSDMA, GSDMB and GSDMC, also have membrane perforation ability and can activate pyroptosis [13]. Gasdermin protein family form pores of 10 20nm on the plasmolemma, subsequently cell contents, including IL-1β and IL-18, slowly leak from the pores, aggravating the in ammatory reaction in the intercellular space. Water molecules entering the cells lead to cell expansion, and rapid destruction of plasma membrane results in cell lysis characterized by nuclear condensation and DNA fragmentation of chromatin [14,15]. In previous studies, pyrophosis has been used primarily to ght infection by pathogens, which ght microbes by releasing damageassociated molecular patterns (DAMPs) that trigger in ammation [16][17][18]. Recent studies have shown that pyrophosis is closely related to antitumor immunity, which can trigger powerful antitumor immunity in vivo and in vitro, and immune checkpoint inhibitors (ICIs) can synergistically improve its e cacy [19,20]. Pyrophosis, though, is usually harmful to normal tissue [21], its exogenous activation induces powerful antitumor activity [22,23]. Cancers of the digestive [24,25], respiratory [26,27], reproductive [28,29] and hematopoietic [30,31] systems are sensitive to pyroptosis induction. Chemotherapeutic drugs, such as paclitaxel and cisplatin, effectively inhibit tumor proliferation and metastasis by inducing pyroptosis [32,33]. Zhang et al. showed that CD8+ T cells and NK cells induced pyroptosis of tumor cells through granase B (an enzyme that can lyse GSDME) in a pyroptosis activated microenvironment, thus establishing a positive feedback loop [34].
According to the existing reports, the role of pyroptosis in tumor development and anti-tumor process is undeniable. However, its role in BC is unclear. Therefore, we conducted a bioinformatics study to compare the expression level of pyroptosis-related genes in normal breast and BC tissues, construct a prognostic model according to the prognostic value of these genes, and study the effect of pyroptosis-related genes on tumor immune microenvironment.

Data preparation
The gene expression data of 1109 BC patients and 113 normal human breast tissue also with the corresponding clinical data were downloaded from TCGA database (https://portal.gdc.cancer.gov/repository).

Identi cation of differentially expressed pyroptosis-related genes
Based on 52 pyroptosis-related genes collected from previous literature, the R software was runned using the "limma" package to identify differentially expressed genes (DEGs) between tumor and normal samples with a P value <0.05. We established a PPI network for the DEGs by Search Tool for the Retrieval of Interacting Genes (STRING) (https://string-db.org/). A correlation graph for these DEGs was constructed applying R software. Construction of prognostic model according to pyroptosisrelated genes To further explore the prognostic value of pyroptosis-related genes in breast cancer, Cox regression analysis was used to analyze the relationship between each gene and the patient's survival status, and screened 7 genes associated with prognosis. We then performed the multivariable Cox regression analysis to narrow down the candidate genes, ultimately a 5-gene signature prognostic model was constructed. The risk score for TCGA expression data was calculated using the "scale" function in R. BC patients in TCGA database were classi ed into high-or low-risk group based on the medium risk score, then we compared the overall survival (OS) time between two groups applying Kaplan-Meier analysis.
PCA based on the 5-gene signature was performed by prcomp function of R. We employed the "survival", "survminer" and "timeROC" R packages to develop the 1, 3 and 5-year ROC curve.
Independent prognostic analysis of the model Clinical features including age, TNM stage and gender from TCGA database combined with risk score were analyzed employing Univariate and multivariable Cox regression analyses.
Functional enrichment analysis of the DEGs between the low-and high-risk groups BC patients in the TCGA data were classi ed into two groups based on the median risk score. We screened out the DEGs between the low-and high-risk groups according to speci c criteria (|log 2 FC| ≥1 and FDR<0.05). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on these DEGs were performed using an online tool called the Database for Annotation Visualization and Integrated Discovery (DAVID)( http://david.abcc.ncifcrf.gov/).

Statistical analysis
We use Single-factor analysis of variance to compare the differences of gene expression between BC tissue and normal breast in TCGA database. The Pearson chi-square test was applied to compare the categorical variables. The OS of patients between high-and low-risk groups was compared employing Kaplan-Meier analysis. When evaluating the independent prognostic value of prognostic model, univariate and multivariate Cox regression was performed. The immune cell in ltration and immune pathway action between two groups were compared by the Mann-Whitney test. All statistical analyses were performed using R software (v4.1.1).

Tumor clustering based on DEGs
A consensus clustering analysis with all 1109 BC patients in the TCGA database was carried out, showing the relations between the expression of the 38 pyroptosis-related DEGs and BC subsets. We increased the clustering variable K from 2 to 10, and the results showed that when k=2, the intragroup correlations were the highest and intergroup correlations were low, suggesting that 1109 BC patients could be well divided into two groups according to 38 DEGs ( Fig. 2A). The survival curve showed signi cant differences in OS between the two clusters (P 0.05) (Fig. 2B). Gene expression levels and clinical features including age (≤65 or 65 years), TNM and stages between the two clusters are shown on the heatmap (Fig. 2C), but no signi cant differences presented in clinical features between the two groups.
Identi cation of prognostic genes and construction of risk models based on TCGA database We performed univariate Cox regression analysis on 1089 BC samples who have complete survival information in TCGA database, and 7 prognosis-related genes were identi ed. Among them, 6 genes (CHMP6, IL18, IRF1, IRF2, TP63, GZMA) were protective genes with HRs 1 while TIRAP HR 1 was associated with high risk (Fig. 3A). By applying the Cox regression analysis, 5 characteristic genes (CHMP6, IL18, IRF2, TP63, TIRAP) were screened out to constructed prognostic model and the risk coef for each gene was calculated. 1089 BC patients were evenly assigned to the high-and low-risk groups based on the median risk score (Fig. 3B). The principal component analysis (PCA) demonstrated that patients could be well classi ed into high-or low-risk group (Fig. 3C). The survival chart showed that the number of deaths increased as the risk increased (Fig. 3D). We performed time-dependent receiver operating characteristic (ROC) analysis to evaluate the sensitivity and speci city of the prognostic model, the result showed that the area under the ROC curve (AUC) was 0.707 for 1-year, 0.669 for 3-year, and 0.663 for 5-year survival (Fig. 3E). Signi cant difference in OS between the high and low risk groups was presented (P 0.001) (Fig. 3F).

Independent prognostic value assessment of prognostic model
To evaluate whether the prognostic model based on gene signatures could be used as an independent prognostic factor, univariate and multivariable Cox regression analyses were applied. Both analyses demonstrated risk score as an independent prognostic factor (P 0.01) univariate HR=1.926, 95% CI: 1.516-2.447; multivariable: HR=2.024, 95% CI: 1.553-2.639 ( Fig. 4A and 4B). In addition, a heatmap showed the distribution of clinical features and signature genes in the high-and low-risk groups according to TCGA (Fig. 4C), in which we found that patient age and T stage distributed differently between the two groups (P 0.01).
Functional enrichment analysis of prognostic models By using R software, FDR 0.05, |log 2 FC| 1 as screening conditions, we identi ed 19 DEGs between highand low-risk groups. All the DEGs were downregulated in high-risk group. GO enrichment analysis and KEGG pathway analysis based on these DEGs were carried out, suggesting that the DEGs mainly enriched in immune response, chemokine-mediated signaling pathway, in ammatory response, cellular response to tumor necrosis factor, B cell receptor signaling pathway, positive regulation of B cell activation, Cytokine-cytokine receptor interaction and Chemokine signaling pathway ( Fig. 5A and 5B).

Comparison of immune function between high and low risk groups
The single-sample gene set enrichment analysis (ssGSEA) was used to assess the immune cell in ltration and immune pathway action of each sample between the high -and low-risk groups in TCGA data. The distribution of immune cells showed that most immune cells were reduced in the high-risk group,

Discussion
In this study, we rst identi ed 38 differentially expressed pyroptosis-related genes between BC and normal tissue based on TCGA database. Signi cant differences in OS of two clusters classi ed according to 38 DEGs were presented, but not in clinical features. Subsequently, to further evaluate the prognostic value of these pyroptosis-related genes, a prognostic model based on 5 gene signatures was constructed, applying the Cox regression analysis. The AUC of model for 1, 3 and 5 years were all greater than 0.6, indicating high accuracy of model predictions. Enrichment analysis suggested that DEGs between highand low-risk group were closely related to immune function. Analysis of immune cell in ltration and immune function between the two groups showed that the in ltration of immune cells and the level of immune function were signi cantly reduced in the high-risk group.
Pyroptosis, an in ammatory -related programmed cell death, was demonstrated to affect the development and treatment of tumor from two aspects. On the one hand, the activation of pyroptosis results in the release of in ammatory mediators, such as IL-1 and IL-18, which might promote cancer development and progression [35,36]. some researchers even regard pyroptosis as another protumorigenic mechanism [37]. On the other hand, induction of pyroptosis has been shown to eradicate tumor cells from a variety of cancers [38]. In BC, whether pyroptosis-related genes affect the OS of patients and how to work on the tumor remain unclear. Our study built a prognostic model containing 5 pyroptosis-related characteristic genes (CHMP6, IL18, IRF2, TP63, TIRAP), which could access the risk and predict the prognosis of BC patients. Charged multivesicular body protein 6 (CHMP6) is a human ortholog of yeast vacuolar protein sorting (vps) 20, an acceptor for ESCRT-III on endosomal membranes and observed to regulate cargo sorting [39]. Fu et al. con rmed that CHMP6 induces oncosis, a form of cell death [40]. Oncosis is remarkable in ischemia, trauma and forms of neurodegeneration [41,42]. However, no relationship between CHMP6 and tumors has been found, and the mechanism by which it induces oncosis is not fully understood. In this study, CHMP6 is a pyroptosis-related gene that reduces the risk of BC survival, and we speculate that it may eradicate cancer cells by inducing oncosis and pyroptosis. Although the speci c mechanism remains to be further studied, it is likely to be a regulatory factor and therapeutic target for BC. Interleukin 18 (IL-18) is a pleiotropic cytokine that produced by a variety of immune cells, but dendritic cells and macrophages are the primary sources of active IL-18 [43]. The main function of IL-18 is mediated through induction of interferon γ(IFN-γ) secretion from T helper (Th1) cells [44]. Studies have shown that IL18 has a dual effect on tumor cells which may due to the complexity of the human immune system [45]. On the one hand, some preclinical models proved that IL18 has antitumor activity [45]. It acts as a proin ammatory factor, recruiting immune cells to the tumor microenvironment to ght tumors, which matched our model that patients with increased immune cell in ltration in BC tissue have a reduced risk of survival. Existing evidence indicates that IL18 has a protective effect on colon cancer [46]. On the other hand, IL18 could promote metastasis by upregulating VEGF and CD44 [47]. However, whether the speci c mechanism is just to induce pyroptosis of cancer cells needs further study. IRF2, a member of the IRF family,is a transcription factor [48]. Studies suggested that IRF2 positively regulated MHC-I pathway and negatively regulates PD-L1 expression [49]. The MHC-I pathway is essential for immune recognition and tumor elimination by CD8 + T cells. Progressing tumors are frequently associated with the MHC-1 de ciency [50] and downregulating the MHC-1 pathway is one way for cancers evading the immune system [51]. Another way cancer can avoid immune elimination is upregulating PD-L1, thus inhibiting Ag-speci c CD8 + T cell function [52]. Therefore, IRF2 not only promotes pyroptosis, but also regulates immune microenvironment to ght tumors. Tumor protein 63 (TP63) is a member of the p53 family of transcription factors [53]. Previous studies showed that tumor cells acquired super-enhancers by recruiting master transcription factors (TFs) to activate the oncogene expression [54,55]. TP63 was con rmed as master TFs of squamous cell carcinoma (SCC) [56,57], thus regulating super-enhancers to promote SCC growth. Interestingly, others have reported that loss of TP63 expression occurs in a large number of humans early and advanced cutaneous squamous cell carcinomas [58]. In addition, low expression of TP63 has also been found associated with shorter survival in other SCCs, which is consistent with our research that TP63 is a protective gene for BC. This phenomenon may be related to the fact that different subtypes of TP63 have different or even opposite functions on cells. To apply TP63 as a predictive target for BC, more research would be needed. Toll-Interleukin 1 Receptor (TIR) Domain Containing Adaptor Protein (TIRAP) is an important adaptor protein which belongs to the TLR/IL-1R superfamily [59], containing two domains, a phosphatidylinositol 4,5-bisphosphate (PIP2) binding domain (PBD) and a TIR domain. TIR Domain is mainly involved in a variety of protein-protein interactions associated with in ammation [60]. Our study also demonstrated that TIRAP is a cell in ammatory death-related gene in BC. It has been proved that the abnormal expression of TIRAP leads to a variety of tumors, such as lymphocytic leukemia [61] and gastric cancer [62]. Meanwhile, its expression is associated with a higher risk and lower survival of BC in our study. However, Whether and how TIRAP plays a role in the tumor immune microenvironment remains unclear.
So far, the research on pyroptosis is not comprehensive enough. Three of the ve pyroptosis-related genes (IL18, IRF2, TIRAP) identi ed by our model could regulate in ammatory response, and we hypothesized that they affect the pyroptosis of BC cells through this way. Pyroptosis-related DEGs between the high-and low-risk groups in BC patients were analyzed for functional enrichment, and we found that they are mainly involved in immune response, chemotaxis of in ammatory cytokines and activation of in ammatory cells, indicating that pyroptosis could regulate the antitumor immunity of BC. Our study found that patients in the high-risk group of BC showed a decrease in immune cell in ltration and immune function activity, suggesting that impaired immune function in the high-risk group was associated with reduced anti-tumor ability. Interestingly, Treg cells have been reported to inhibit anti-tumor immunity and are associated with poor prognosis [63], while it reduced in high-risk group in our study, which may due to the different roles of its two subtypes. Therefore, it is necessary to identify the subtype of Treg cells in BC.

Conclusion
To sum up, pyroptosis-related genes are closely related to BC, due to its differentially expressed between BC and normal tissue. In addition, a prognostic model of BC according to 5 signature genes picked out based on TCGA database, and this model was proved to be an independent prognostic factor of BC. There were signi cant differences in immune cell in ltration and immune function activity between the high-and low-risk groups of the model, suggesting that pyroptosis-related genes have an effect on BC antitumor immunity. Our study provides novel predictive marker genes for the prognosis for BC and supplies new clues for immunotherapy of BC.

Consent for publication
All the authors agreed to publish the manuscript

Con icts of Interest
The authors report no con icts of interest in this work. the interactions of the pyroptosis-related genes (interaction score=0.4). C The correlation network of the pyroptosis-related genes (red line: positive correlation; blue line: negative correlation. The depth of the colours re ects the strength of the relevance).