Molecular Characteristics, Clinical Implication, and Cancer Immunity Interactions of Pyroptosis-Related Genes in Breast Cancer

Objective: Pyroptosis represents an emerging inflammatory form of programmed cell death. Herein, specific functions and clinical implications of pyroptosis-related genes were systematically characterized in breast cancer. Methods: Expression, somatic mutation and copy number variation of 33 pyroptosis-related genes were assessed in breast cancer from TCGA dataset. Their interactions, biological functions and prognostic values were then observed. By stepwise Cox regression analysis, a pyroptosis-related gene signature was generated. The predictive efficacy in survival was examined by survival analyses, ROCs, univariate and multivariate analyses and subgroup analyses. Associations between risk score (RS) and cancer immunity cycle, HLA, immune cell infiltrations, and immune checkpoints were analyzed. Results: Most of pyroptosis-related genes were abnormally expressed in breast cancer. CASP8, NLRC4, NLRP3, NLRP2, PLCG1, NLRP1, NLRP7, SCAF11, GSDMC, and NOD1 occurred somatic mutations as well as most of them had high frequency of CNV. There were closely interactions between them. These genes were distinctly enriched in immune-related processes. A three-gene signature was generated, containing IL-18, GSDMC, and TIRAP. High RS predicted poorer overall survival, progression, and recurrence. After verification, this RS was an independent and sensitive predictive index. This RS was negatively correlated to cancer immunity cycle. Also, low RS was characterized by high HLA, immune cell infiltrations and immune checkpoints. A nomogram including age and RS was generated for accurately predicting 5-, 8-, and 10-year survival probabilities. Conclusion: Pyroptosis-related genes exert key roles in cancer immunity and might be applied as a prognostic factor of breast cancer.


INTRODUCTION
Breast cancer represents a frequently diagnosed malignancy among women globally, with a high mortality (1). This malignancy affects 1/20 globally and 1/8 in high-income countries (2). Females with high risk of developing breast cancer are a heterogeneous population (3). Further research requires to improve prognostic models to stratify high-risk patients. The biology in breast cancer progress is complex in which genetic and environmental elements are involved (4). Conventional breast cancer classification primarily replies on clinicopathologic characteristics and routine markers, not capturing various clinical courses of individual patient (5). Indepth understanding of the molecular mechanisms could lead to improvement in patients' prognoses.
Varied factors are in relation to carcinogenesis, such as activations of proto-and antioncogenes, immune microenvironment as well as chronic inflammation (6,7). Pyroptosis represents a form of programmed cell death, which may induce the cleavage of gasdermin D along with activation of immune and inflammatory response (8). Activated pyroptosis induces the release of the inflammatory factors IL-1 and IL-18, thereby promoting breast cancer initiation. It is featured by cell swelling as well as bubble-like protrusions. Interactions between pyroptosis and cancers are complex due to varied influence of pyroptosis on cancers in distinct tissue specimens as well as genetic background (9). Increasing evidences highlight the roles of pyroptosis in carcinogenesis (10). Previously, increased GSDMB expression is related to poor survival and high metastases of breast cancer (11). Furthermore, recent research has demonstrated the crosstalk between pyroptosis and anti-cancer immunity (12). Research has displayed that chemotherapy drugs, miRNAs, etc., may trigger cancer pyroptosis, which inhibits malignant development of cancers (13). Based on pyroptosis regulators, a seven-gene signature has been generated for predicting ovarian cancer prognoses (14). However, no studies have reported prognostic implications of pyroptosis-related gene signature in breast cancer.
Herein, we analyzed molecular characteristics and clinical implication of pyroptosis-related genes as well as their interactions with cancer immunity in breast cancer.

Data Acquisition
RNA-seq profiling (FPKM values) of breast cancer was retrieved from the Cancer Genome Atlas (TCGA; https://cancergenome. nih.gov/) database. After removing samples without complete clinical information, 1,082 breast cancer samples were included in our study (Supplementary Table 1). Meanwhile, RNA-seq profiles of 113 adjacent normal tissues were also obtained from TCGA database. Then, FPKM values were converted to TPM values according to the following formula: TPMi = FPKMi * 1000000/(FPKM0 + . . . . + FPKMm), where i represents gene i and m represents the total number of all genes. Somatic mutation and copy number variation (CNV) data were also retrieved from TCGA database.

Protein-Protein Interaction (PPI) and Gene Ontology (GO) Annotation Analysis
Thirty pyroptosis-related genes were uploaded onto the STRING database (https://string-db.org) and their interaction pairs were obtained (20). A PPI network of pyroptosis-related genes was visualized with Cytoscape (https://cytoscape.org/) (21). The gene set of 33 pyroptosis-related genes was analyzed via clusterProfiler package (22). Biological process, cellular component and molecular function of these pyroptosis-related genes were analyzed. Terms with adjusted p < 0.05 were significantly enriched.

Development of a Prognostic Model
To screen which pyroptosis-related genes were related to breast cancer prognoses, univariate analyses were carried out. Genes with p < 0.05 were included for conducting a stepwise Cox regression model. Risk score (RS) was calculated by linearly combining regression coefficients multiplied with expression values. Based on the median of RS, patients were stratified into two groups. Overall survival was compared between groups by Kaplan-Meier curves and log-rank tests. Expression patterns of genes in this model were visualized into heatmap. Area under the curve (AUC) of the receiver operating characteristic (ROC) curve was conducted for assessing the predictive efficacy of this model.

Univariate and Multivariate Cox Regression Analyses
Univariate analyses were applied for estimating the associations between prognoses and age, T, N, M, stage, and RS. Furthermore, multivariate analyses were presented for observing whether these factors were independently predictive of prognoses of breast cancer. Hazard ratio (HR), 95% confidence interval (CI), and p-values were separately determined.

Subgroup Analysis
To investigate the sensitivity of RS in predicting prognoses, patients were stratified into subgroups according to clinical characteristics, including age ≥65 and age <65, M0 and M1, N0 and N1-3, stage I-II and stage III-IV, T1-2 and T3-4 subgroups. OS was compared in high and low RS patients in each subgroup. P-values were determined with log-rank tests.
Assessment of Activated Pathways, HLA, Immune Cell Infiltration, and Immune Checkpoints Enrichment scores of several cancer-related pathways were estimated by single sample gene set enrichment analysis (ssGSEA) algorithm (23), including IFN-Gamma signature, APM signal, base excision repair, cell cycle, DNA replication, Fanconi anemia pathway, homologous recombination, microRNAs in cancer, mismatch repair, nucleotide excision repair, oocyte meiosis, p53 signaling pathway, progesterone-mediated oocyte maturation, proteasome, pyrimidine metabolism, spliceosome, systemic lupus erythematosus, viral carcinogenesis. The gene sets of above pathways were listed in Supplementary Table 2.
Associations between RS and enrichment scores of pathways were then analyzed with Spearman correlation test. Wilcoxon test was applied for estimating the differences in expression of HLA signatures and immune checkpoints between highand low-risk groups. Also, infiltration levels of 28 immune cells were inferred by ssGSEA algorithm and were compared between groups. Adjusted p-value < 0.05 indicated the significant difference between high-and low-risk groups.

Gene Set Enrichment Analysis (GSEA)
GSEA method was employed for identifying enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) in high and low RS groups based on transcriptomic data (25). Terms with nominal enrichment score >2 and false discovery rates <0.05 were significantly enriched.

Nomogram Construction
Independent prognostic factors were incorporated for constructing a prognostic nomogram in predicting 5-, 8-, and 10-year survival duration with stepwise Cox regression analyses. Calibration plots were plotted for comparing nomogram-predicted and observed 5-, 8-, and 10-year survival.

Cell Culture
Human normal breast cells (MCF-10A) and breast cancer cells (MDA-MB-231 and HCC70) were retrieved from Shanghai Cell Bank of Chinese Academy of Sciences (China), which were maintained in DMEM plus 10% FBS (Gibco, USA) and 1% penicillin-streptomycin. All cells were incubated at 37 • C under the condition of 5% CO 2 .

Statistical Analysis
Data were analyzed with the R (version 3.6.1) and R Bioconductor packages. Comparisons between two groups were evaluated with student' t-test or Wilcoxon test. Differences in disease-free interval (DFI), disease-free survival (DFS), diseasespecific survival (DSS) and progression-free interval (PFI) were compared in high and low RS groups with Kaplan-Meier curves and log-rank tests. P < 0.05 indicated statistical significance.

Interactions Between Pyroptosis-Related Genes and Their Biological Implications
Spearman test was utilized for evaluating correlations between 33 pyroptosis-related genes in breast cancer from TCGA dataset. In Figure 2A, there were tight interactions between them, such as CASP4 and CASP1. A PPI network was than constructed based on 33 pyroptosis-related genes ( Figure 2B). Among all nodes, CASP1 had the highest degree. Furthermore, biological implications of 33 pyroptosis-related genes were evaluated by clusterProfiler package. Our data showed that these pyroptosis-related genes were mainly involved in regulating IL-1β production and secretion biological processes ( Figure 2C). Furthermore, they participated in cytosolic, inflammasome and membrane cellular components as well as apoptotic process.  (Figure 3A). According to their coefficients and expression levels, we calculated risk score (RS) for each patient (Table 1). Then, we evaluated association between RS and survival outcomes. In Figure 3B, patients with high RS exhibited poorer OS than those with low RS (p = 7.124e-05). Furthermore, we observed expression patterns of IL-18, GSDMC, and TIRAP in high and low RS samples. As a result, IL-18 was down-regulated as well as GSDMC and TIRAP were up-regulated in high RS group compared with low RS group ( Figure 3C). We also showed their expression in different clinical features including stage, age, T, N, and M. Time-dependent ROC curves were plotted for verifying the predictive performance. In Figure 3D, AUC was 0.652, indicating that RS possessed high accuracy and sensitivity in predicting prognoses of breast cancer.

Association Between the Pyroptosis-Related Gene Signature and Immunogenicity
We further assessed whether the pyroptosis-related gene signature affected cancer immunity cycle. As a result, we found that the risk score almost negatively participated in each step of cancer immunity cycle, including release of cancer cell antigens, cancer antigen presentation, priming and activation, T cell recruiting, Th1 cell recruiting, dendritic cell recruiting, Th22 cell recruiting, macrophage recruiting, monocyte recruiting, neutrophil recruiting, NK cell recruiting, eosinophil cell recruiting, basophil recruiting, Th17 cell recruiting, B cell recruiting, Th2 cell recruiting, Treg cell recruiting, infiltration of immune cells into tumors, recognition of cancer cells by T cells, and killing of cancer cells ( Figure 6A). Also, RS was distinctly related to the activation of IFN-gamma signature, APM signal, cell cycle, Fanconi anemia pathway, homologous recombination, microRNAs in cancer, oocyte meiosis, P53 signaling pathway, proteasome, spliceosome, and viral carcinogenesis ( Figure 6A). The low RS breast cancer samples were characterized by increased expression of human lymphocyte antigen (HLA), as  Figure 6B). Also, higher infiltration levels of nearly all immune cells were detected in the low RS group than the high RS group, including activated B cell, activated CD4 T cell, activated CD8 T cell, central memory CD4 T cell, effector memory CD4 T cell, effector memory CD8 T cell, gamma delta T cell, immature B cell, memory B cell, regulatory T cell, T follicular helper cell, type 1 T helper cell, type 17 T helper cell, type 2 T helper cell, activated dendritic cell, CD56bright natural killer cell, CD56dim natural killer cell, eosinophil, immature dendritic cell, macrophage, mast cell, MDSC, monocyte, natural killer cell, natural killer T cell, and neutrophil ( Figure 6C). Moreover, we compared the expression of immune checkpoints between high and low RS groups. In Figure 6D, higher expression of BTLA , CD200R1, CD244, CD27,  CD28, CD40, CD40LG, CD48, CD70, CD86, CTLA4, HAVCR2,  HHLA2, ICOS, IDO1, IDO2, KIR3DL1, LAG3, LAIR1, LGALS9,  PDCD1, PDCD1LG2, TIGIT, TMIGD2, TNFRSF14, TNFRSF18,  TNFRSF25, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF14, TNFSF9, and VSIR was found in the low RS group. Meanwhile, the high RS group was characterized by increased expression of CD160, CD276, CD44, ICOSLG, and NRP1. Above data suggested that RS was related to suppressive immunity in breast cancer.
Pathways Involved in the Pyroptosis-Related Gene Signature GSEA was employed for exploring pathways involved in the pyroptosis-related gene signature. In Figures 7A-F, high RS was distinctly associated with cell cycle (NES = 1.85 and FDR =

Development of a Prognostic Nomogram for Breast Cancer
Our multivariate analyses demonstrated that age and the pyroptosis-related gene signature were independent risk factors of breast cancer, which were used for developing a prognostic nomogram. In Figure 8A, this nomogram could predict 5-, 8and 10-year survival probabilities. The predictive performance was evaluated by calibration plots. There were high consistencies in nomogram-predicted and actual 5-, 8-and 10-year survival (Figures 8B-D). Our data suggested that the nomogram exhibited the well predictive efficacy in 5-, 8-and 10-year survival probabilities.

Validation of the Expression of Genes in the Pyroptosis-Related Gene Signature
We further validated the expression of IL-18, GSDMC, and TIRAP in the pyroptosis-related gene signature in vitro. As depicted in western blotting, IL-18, GSDMC, and TIRAP expression was all markedly up-regulated in breast cancer cells MDA-MB-231 and HCC70 compared with normal breast cells MCF-10A (Figures 9A-D). Immunofluorescence also confirmed the significant up-regulation of IL-18, GSDMC and TIRAP expression in MDA-MB-231 and HCC70 cells than MCF-10A cells (Figures 9E-G).

DISCUSSION
Programmed cell fate nearly focuses on apoptosis and necroptosis. It is of importance to find an alternative option when these cell deaths are compromised (26). Pyroptotic death represents a form of programmed cell death, which is induced by inflammasomes (27). Inducing apoptosis of cancer cells has been applied for eliminating malignant cells (28). Nevertheless, due to escaping apoptosis, induction of pyroptosis may be especially critical in treating antiapoptotic cancers. Immunotherapies show remarkable efficacy in treating breast cancer. Nevertheless, therapeutic effects are still limited (29). Pyroptosis offers an opportunity for alleviating immunosuppression as well as promoting an immune response in treating breast cancer (30).
Here, we characterized expression, genetic mutations, and clinical implications of pyroptosis-related genes in breast cancer. Among pyroptosis-related genes, GPX4, NLRP7, CASP6, CASP3, IL1B, IL18, CASP8, NLRP6, IL6, GSDMC, PYCARD, AIM2, NOD2, NLRP3, and CASP4 were aberrantly expressed in breast cancer. Furthermore, CASP8 (2%), NLRC4 (1%), NLRP3 (1%), NLRP2 (1%), PLCG1 (1%), NLRP1 (1%), NLRP7 (1%), SCAF11 (1%), GSDMC (1%), and NOD1 (1%) occurred somatic mutations as well as most of them had high frequencies of CNV in breast cancer. Our spearman test and PPI network both revealed the tight interactions between pyroptosis-related genes. In the PPI network, CASP1 had the highest degree. Consistently, CASP1 is a prognostic factor as well as therapeutic target in breast cancer (31). As shown in GO enrichment analyses, pyroptosis-related genes were primarily involved in mediating IL-1β production and secretion biological processes as well as cytosolic, inflammasome and membrane cellular components and apoptotic process, indicating the key biological implications of pyroptosis in tumorigenesis. We generated a pyroptosisrelated gene signature, containing IL-18, GSDMC, and TIRAP. High RS was indicative of undesirable OS, recurrence, and progression of breast cancer. AUC = 0.652 demonstrated the well predictive efficacy. Multivariate analyses and subgroup analyses suggested that this RS was an independent risk factor of breast cancer prognoses. Previously, a pyroptosisrelated gene signature was generated for prediction of ovarian cancer prognoses (14). Our univariate cox regression analyses showed that IL-18 was a protective factor as well as GSDMC and TIRAP were risk factors for breast cancer prognoses. IL-18, a proinflammatory cytokine, modulates inflammation, and immune response. As confirmed by previous studies, mesenchymal stem cells expressing IL-18 suppresses breast cancer proliferation and metastasis, suggesting the antitumor activities of IL-18 (32,33). GSDMC is specifically cleaved by caspase-8 through TNFα to generate GSDMC N-terminal domain, thereby forming pores on the cell membrane as well as inducing pyroptosis. High expression of GSDMC is in relation to undesirable survival of cancer patients (34). Furthermore, there is evidence shows that TIRAP is a risk factor of cancer prognosis (35).
Immunotherapies especially immune checkpoint inhibitors (ICIs) may produce durable therapeutic effects. Nevertheless, only one third of patients respond to ICIs. Breast cancer is often considered as a cold tumor, with lowly frequent mutation, decreased immune cell infiltration, and suppressive immune microenvironment (36). Inducing cell deaths other than apoptosis has been considered as novel cancer therapeutic strategies due to innate resistance to apoptosis (37). Combination of inducing pyroptosis and ICIs could display synergistically increased anti-cancer activity (28,38,39). However, most evidences of interaction between immunity and pyroptosis are derived from animal and cellular model. Here, pyroptosis-related gene signature was negatively correlated to almost all steps of cancer immunity cycle (40). Also, we found that low RS was characterized by high HLA, immune cell infiltrations and immune checkpoints. This suggested that pyroptosis exerted an impact on immuno-oncology.
In-depth exploring pyroptosis mechanisms, up-and downstream pathways may offer novel insights into breast cancer therapy. Our GSEA results suggested that high RS was in relation to carcinogenic pathways including cell cycle, ERBB signaling pathway, mTOR signaling pathway, TGF-beta signaling pathway, ubiquitin mediated proteolysis and WNT signaling pathway. Also, low RS was significantly related to immune-related pathways including autoimmune thyroid disease, cytokinecytokine receptor interaction, primary immunodeficiency, and ribosome. These data indicated the interactions of pyroptosisrelevant RS with above pathways in breast cancer progression. Epidemiologic studies have confirmed that age is a risk factor of breast cancer (41). Personalized medicine is based on individual evaluation of risk. By including two independent risk factors age and this RS, we generated a nomogram for prediction of 5-, 8-, and 10-year survival probabilities. The predictive accuracy was confirmed by comparing observed survival duration.
There are several limitations in our study. Firstly, due to the limited clinical features of patients, we cannot carry out subgroup analyses by stratifying more factors. Secondly, the pyroptosis-relevant RS was constructed and verified based on retrospective cohorts. In conclusion, our findings revealed that pyroptosis induction might be a novel strategy for breast cancer immunotherapy, characterized by high compatibility and extensive clinical applicability. In future studies, prognostic implications of pyroptosis will be observed in a larger breast cancer cohort. Also, interactions of pyroptosis with cancer immunity will be further verified in cellular and animal models.

CONCLUSION
Collectively, our data characterized expression patterns and mutations of pyroptosis-related genes. A three-gene regression model including IL18, GSDMC, and TIRAP was regarded as an independent risk factor of breast cancer prognoses. The RS was distinctly cancer immunity cycle, HLA immune cell infiltration and immune checkpoints in breast cancer. Our data suggested that pyroptosis combining with immunotherapies might be a potential therapeutic strategy.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.