Pyroptosis is a newly discovered form of cell programmed necrosis, but its role and mechanism in cancer cells remain unclear. The aim of this study is to systematically analyze the transcriptional sequencing data of breast cancer (BC) to find a pyroptosis-related prognostic marker to predict the survival of BC patients.
The original RNA sequencing (RNA-seq) expression data and corresponding clinical data of BC were downloaded from The Cancer Genome Atlas (TGCA) database, followed by differential analysis. The pyroptosis-related differentially expressed genes (DE-PRGs) were employed to perform a computational difference algorithm and Cox regression analysis. The least absolute shrinkage and selection operator (LASSO) was utilized to avoid overfitting. A total of 4 pyroptosis-related genes (PRGs) with potential prognostic value were identified, and a risk scoring formula was constructed based on these genes. According to the risk scores, the patients could be classified into high- and low-risk score groups. The potential molecular mechanisms and properties of PRGs were explored by computational biology and verified in Gene Expression Omnibus (GEO) datasets. In addition, the quantitative real time PCR (RT-qPCR) and Human Protein Atlas (HPA) were performed to validate the expression of the key genes.
A PRGs signature, which was an independent prognostic factor, was constructed, and could divide patients into high- and low-risk groups. The results from the prognostic analysis indicated that the survival was significantly poorer in the high-risk group than in the low-risk group both in TCGA and in GEO, indicating that the signature is valuable for survival prediction and personalized immunotherapy of BC patients.
The pyroptosis-related biomarkers were identified for BC prognosis. The findings of this study provide new insights into the development of the efficacy of personalized immunotherapy and accurate cancer treatment options.
Breast cancer is the most common malignancy in females and one of the three most common tumors worldwide . Breast cancer is a heterogeneous subtype of tumor with poor prognosis [2,3,4]. Due to the development of more effective and superior medical diagnostic and imaging techniques, the mortality rate of BC has been greatly reduced, but the prognosis of patients with BC is still poor [1, 5]. The lack of effective features and diagnostic tools to predict prognosis or long- term survival in patients remains a major obstacle to improve detection and treatment strategies for BC .
Pyroptosis distincting from apoptosis, is accompanied by inflammation and immune response, and it is a new form of Gasdermin (GSDM) family-mediated programmed cell death [7, 8].. When bacteria, fungi or parasites invade, immune cells, such as lymphocytes and neutrophils actively kill pathogens through a series of signal transduction. Many pathogens invade and hide in host cells in order to avoid detection by antiseptic substance and phagocytes in body fluids. To eliminate these pathogens, the solution is to clean out them together with the infected cells. Killing infected cells can be done by cell intrinsic mechanisms like necroptosis, apoptosis, and pyroptosis . Gasdermin family member contains gasdermin A, B, C, D, E, and DFNB59 [6, 10]. GasderminD (GSDMD) represents a big gasdermin family, with a new membrane pore forming activity. GSDMD, the substrate of both caspase-11/4/5 and caspase-1, is by far the best researched . Pyroptosis was originally not correctly appreciated for decades because it was similar to apoptosis, and it was condemned as a special type of apoptosis through caspase-1 . Caspase-1 and caspase-11/4/5 induced pyroptosis by cleavage GSDMD, release its gasdermin-N structure domain, the domain structure has the activity of punching holes in the membrane, eventually leading to cell swelling and osmotic lysis. The cells undergo the morphological changes described above . After the demonstration of the GSDMD-mediated pathway, other pyroptosis mechanisms, such as caspase-3/8-mediated pathway and granzyme-mediated pathway, have been clarified by several studies. Chemotherapy can induce caspase-3-mediated cleavage of GSDME, and form N-GSDME terminal, which can cause pyroptosis of tumor cells . Caspase-8 specifically cleaves GSDMC to produce N-GSDMC, and forms pores in cell membrane to induce pyroptosis . Recent studies have proved that GzmB can further activate anti-tumor immune response and inhibit tumor growth by activating caspase-3/GSDME or directly cracking GSDME and inducing pyroptosis [16, 17]. Pyroptosis shows different morphology compared with apoptosis, it is lytic, featuring cell swelling under microscope [18, 19]. Recently, pyroptosis has become a research hotspot in the occurrence and development of tumors. Besides, it is reported to be closely related to gastric cancer, colorectal cancer, hepatocellular carcinoma, breast cancer, and lung cancers [20,21,22,23,24,25], but its role and mechanism in cancer cells remain unclear.
More and more studies have shown that tumor microenvironment (TME) and tumor stemness are closely related to BC occurrence and development [26, 27], infiltration of numerous inflammatory cells in BC, and the density of CD8+ T cells is highly related to the immune escape of BC [28, 29]. PD-1 and PD-L1 constitute an essential inhibitory mechanism which causes T cell exhaustion in tumor microenvironment. That’s the main reason why PD-L1 has drawn increasing attention of researchers [30,31,32]. But the underlying mechanisms of pyroptosis in breast cancer microenvironment progression and immune response remain unclear. This study mainly aimed to explore PRGs in BC, and systematically investigate the association between the pyroptosis-related gene signature and immune microenvironment, immune cell infiltration, cancer chemoresistance, cancer stem cells (CSCs). These results supported the feasibility of constructing tumor prognostic risk models using PRGs. At present, few studies have reported the prognostic value of PRGs in BC in recent years.
In this study, the mRNA expression profiles and clinical data of BC patients were first downloaded from TCGA to identify differentially expressed genes (DEGs), especially pyroptosis-related genes. Then, a prognostic signature with these genes was constructed and its reliability was validated in the GEO database. Finally, this gene signature was proved to be able to predict the prognosis of BC and assess the patient’s tumor microenvironment and other states, thereby contributing to clinical treatment.
Material and methods
The RNA sequencing (RNA-seq) expression data and clinicopathological information of female breast cancer patients from 1053 breast cancer tissue samples and 111 nontumor tissue samples were downloaded from the TCGA BC dataset (https://portal.gdc.cancer.gov/), and were used as training cohort. Probes were transformed to corresponding Entrez gene names referring to the annotation files. 33 genes associated with pyroptosis were identified from previous literature . In order to get more breast cancer datasets, the GSE42568 and GSE86166 datasets, which were obtained from the gene expression omnibus (GEO: https://www.ncbi.nlm.nih.gov/geo/) database, were used as testing cohort. Batch normalization was applied by using ‘sva’ and ‘limma’ R package . A total of 470 breast cancer samples were obtained. The detailed flow-process diagram of this study is shown in Fig. 1.
Construction of a pyroptosis-related gene signature
First, the expression level of pyroptosis-related genes was extracted from the total gene expression list. If a gene appeared more than once in the same sample, the ‘limma’ of Bioconductor R package was utilized for averaging operations . Second the limma was utilized to identify DE-PRGs between breast cancer tissue samples and normal breast tissue samples. The false discovery rate (FDR) threshold was set at FDR < 0.05 for DE-PRGs calling. A protein-protein interaction (PPI) network of proteins encoded by DE-PRGs of high-risk epidermal growth factor receptor 2-positive (HER2+) and triple-negative/basal-like molecular subtypes was visualized using String (http://string-db.org). To establish the pyroptosis-related gene risk model, univariate Cox regression analysis was performed on the pyroptosis-related genes. A total of 4 prognostic related differential genes were obtained by the intersection of DE-PRGs and prognostic genes.
To avoid overfitting, the least absolute shrinkage and selection operator (LASSO) was utilized to select variables with high prognostic value . Next, 1000 LASSO iterations were performed for prognostic model construction using the ‘glmnet’ package in R, and their regression coefficients were obtained. Finally, the formula of the risk score was composed as follows, and risk scores were computed: Risk score = ∑ni = ∑Coefi × xi, where xi represents the normalized expression level of target gene i and Coefi represents the regression coefficient. According to the median risk score in TCGA dataset, 1014 patients in the data set were divided into high-risk and low-risk groups after samples with a survival time of zero were removed. The Kaplan-Meier (K-M) plot was used to evaluate survival differences between the high- and low-risk groups. To analyze the distribution differences between different groups, PCA was performed using the ‘prcomp’ function in the STATS package in R. A t-SNE analysis was implemented using the R package Rtsne (https://github.com/jkrijthe/Rtsne).
Univariate and multivariate cox regression analysis
Univariate cox regression analysis was presented for assessment of the prognostic values of the risk score and clinical features (Age, Stage, T classification, N classification, M classification). Then, multivariate cox regression analysis was used to determine which prognostic factors could independently predict the survival of patients. Adjusted p < 0.05 is considered to be statistically significant using the ‘survival’ package.
Functional enrichment and pathway analysis
To further investigate the biological processes associated with the pyroptosis-related genes, BC patients were divided into the high- and low-risk groups based on the median risk score in TCGA and testing cohort, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for all selected DEGs between the two risk . Cohorts were performed with the ‘clusterProfiler’ package in BioConductor using |log2FC| ≥ 1 and FDR < 0.05 as thresholds in TCGA. Considering the relatively small sample size in testing cohort, the threshold was set as FDR < 0.05.
Estimation of TME cell infiltration, immuneScore, stromalScore, PD-L1, tumor stemness, drug sensitivity
The scores of 16 tumor-infiltrating immune cells and 13 immune-related functions for samples were determined by single-sample gene-set enrichment analysis (ssGSEA). The ImmuneScore, StromalScore, and ESTIMATEScore were calculated using the ‘ESTIMATE’ package [38, 39]. Correlations between the risk signature and the key immune regulators, PD-L1 and PD-L2 were evaluated. The DNA index is a score based on methylation data and the RNA index is a score based on transcriptome data in TCGA, they can reflect the amount of stem cells [40, 41]. The NCI-60 database and information on 216 FDA-approved chemotherapy drugs were obtained from the CellMiner interface (https://discover.nci.nih.gov/cellminer). Spearman correlation analyses were used to measure the relationship among the risk score, ImmuneScore, StromalScore, PD-L1 and PD-L2 expression, tumor stemness, and drug sensitivity.
Analysis based on human protein atlas database
The HPA database covers all pathological and gene expression data collected from a large number of studies using different cell lines and tissue types . Immunohistochemistry images in this database were implemented in the present work to examine 4 PRGs levels within diverse tissues along with their localization in cells.
Cell culture and reagents
Human breast cancer cell lines MCF7, MCF-10A, HCC1937, MDA-MB-231 were provided by Affiliated Hospital of Qingdao University. HCC1937 and MDA-MB-231 ware cultured in DMEM (Invitrogen, USA), MCF7 cells were maintained in RPMI 1640 (Gibco, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco, USA) and 1% penicillin-streptomycin (PS, 100 μg/ml) (Enpromise, Hangzhou, China), MCF-10A was maintained in Medium-F12 (DMEM/F12) (Gibco, USA) supplemented with 5% horse serum (Gibco, USA), 1% penicillin/streptomycin, 0.5 μg/ml hydrocortisone, 100 ng/ml cholera toxin (Sigma, USA), 10 μg/ml insulin (Gibco, USA), and 20 ng/ml recombinant human EGF (Invitrogen, USA). All cells were cultured in a humid environment of 37 °C and 5% CO2.
RNA isolation and quantitative real-time polymerase chain reaction PCR (RT-qPCR)
Total RNA was extracted from cells using TRIzol reagent (Invitrogen, USA). Complementary DNA (cDNA) was synthesized using the total RNA and a PrimeScript RT reagent kit (Takara). TB- Green assays (Takara) were used to perform the RT-qPCR on a Roche LightCycler® 480 instrument. The data was calculated through the 2-ΔΔCt strategy, normalizing with GAPDH. The primer sequences used for qRT-PCR in this study are listed in Table 1.
We used R software (version 4.0.3) to perform all statistical analyses. The Student’s t-test was used to compare gene expression levels between BC samples and non-cancer samples. Heatmaps of the LASSO analysis genes were plotted using the ‘heatmap’ R package. R packages ‘survival’ and ‘survminer’ were used for survival analysis [43, 44]. The OS for the two risk groups was evaluated by Kaplan-Meier (K-M) survival curves and log-rank test. Bioconductor R package ‘GSVA’ was used to compare ssGSEA enrichment scores for immune cells and immune-related pathways between the two groups (i.e. high- or low-risk groups) . Unless otherwise stated, p < 0.05 was considered statistically significant. P values were showed as: ns not significant; *P < 0.05;**P < 0.01; ***P < 0.001.
Identification of prognostic PRGs in the breast TCGA cohort
Threshold was set at FDR < 0.05 to compare PRGs expression level between breast cancer and normal tissue. A total of 22 differently expressed genes was obtained. In the univariate Cox regression model, we found 5 genes were associated with a significant OS. Subsequently, a total of 4 intersection genes (GPX4, GSDMD, GSDMC, IL18) was selected as hub genes for further analyses (Fig. 2a). The PPI network of HER2+ and basal-like molecular subtypes that indicates tight interplay of pyroptosis-related genes are shown in Fig. 2b and c. Additionally, the prognosis of 4 genes was shown in Fig. 2d and the expression profiles of the 4 genes were showed in a heatmap (Fig. 2e). In order to compare four genes expression level in clinical cases, and to explore the clinical significance of the signature. The mRNA expression level was showed in Fig. 2f, and GPX4 showed a low expression trend, while IL18 showed a high expression.
Construction of a prognostic PRGs signature
According to the result of the LASSO, we ended up with 4 key PRGs which related to prognosis for building the prognostic signature of breast cancer. Figure 3a shows the risk score distribution of patients. Figure 3c shows the survival status of high and low risk group of patients in TCGA database. With the increasing of risk score, the patient’s survival time reduced, on the contrary, the death risk increased. High-risk group of patients than low-risk groups have a greater probability of death incidents. As we can see from t-SNE mappings and PCA (Fig. 3e and g), patients have formed two different clusters. The result of Kaplan-Meier plot shows that the survival rate of high-risk group was obviously lower than low-risk group (Fig. 3i).
Validation of the prognostic model in the validation cohort
To further evaluate the accuracy of the PRGs signature in predicting BC prognosis, it was validated by us using the same methods in the testing cohort. The validation set is segregated into high (N = 235) and low (N = 235) risk groups, Each patient’s survival outcome, risk status were demonstrated in Fig. 3b, d, f, h. K-M analysis shows that patients in the high-risk group also had a worse prognosis than those in the low-risk group (P < 0.05, Fig. 3j). Similarly, the survival rate of high-risk group was significantly lower than that of low-risk group in basal-like molecular subtypes and luminal subtype (P < 0.05, Fig. 4a and b).
Independent prognostic value of the 4-gene signature
We observed clinical factors and gene signature prognostic significance through univariate and multivariate regression. Samples have been Chosen with complete clinical information. In 867 cases of patients, according to the age, clinical stage, histological grade and clinical pathologic factors, risk parameters were explored for patients. We have defined these variables indicated significant differences in univariate analysis and stage, age N-classification, M-classification showed significant differences in multivariable analysis. Risk parameters were important prognostic values of p < 0.05 (Fig. 5a, b and Table 2).
Functional enrichment analyses
GO and KEGG functional enrichment analyses were performed on risk-related DEGs to investigate the potential functions. The result of GO analyse indicated that the DEGs equally concentrated in membrane raft, membrane microdomain, membrane region and external side of plasma membrane (Fig. 6a, b). And KEGG functional enrichment analysis suggested that the DEGs were mainly related to Cytokine-cytokine receptor interaction, Hematopoietic cell lineage, etc. both in TCGA and GEO database (Fig. 6c, d). To further explore the relationship between BC prognosis and immune status, we quantified immune cell infiltration score and immune-related function using ssGSEA. The correlations between ssGSEA scores and different risk groups showed that the scores of iDCs, aDCs, CD8 + Tcells, T helper cells, NK cells, Macrophages, Th2-cells, Treg were higher in the low-risk group (Fig. 7a, b). Meanwhile, APC costimulation, cytolytic_activity, inflammation−promoting, parainflammation, etc. were significantly different between the low- and high-risk groups in both TCGA and GEO database (Fig. 7c, d).
Associations with immunity, tumor stemness, and drug sensitivity
The constructed risk signature was significantly negatively correlated with the immune and stromal scores (Fig. 8a, b). Furthermore, there are no significant correlation between PyroptosisScore and DNAss (Fig. 8c), but positively correlated with RNAss (Fig. 8d). When studying the relationship with immune checkpoints, considering the role of PD-L1 (also known as CD274) and PD-L2 (also known as PDCD1LG2) in immune microenvironment and immune escape, we analyzed the difference in the expression of these two proteins between high and low risk groups. The results showed that the protein expression levels of PD-L1 and PD-L2 in the high risk group were significantly lower than those in the low risk group (Fig. 8e, f), and the protein expression level was negatively correlated with the risk score (Fig. 8g, h). As shown in Table 3, the PRGs are resistant to most drugs (p < 0.05). For example, Paclitaxel, Vinorelbine, Gemcitabine and Epirubicin are commonly used drugs to treat breast cancer, the expression levels of IL18 were negatively associated with tumor cell sensitivity to Paclitaxel, Vinorelbine and Epirubicin. In contrast, sensitivity to the chemotherapy drug Gemcitabine was positively associated with the expression levels of GSDMD in tumor cells (Fig. 9).
Expression levels of key genes in the clinical samples and RT-qPCR
To confirm the bioinformatics prediction, we detected the expression of the four PRGs by RT-qPCR. It showed a low mRNA expression trend of GPX4, GSDMD and GSDMC in all three types of breast cancer cells, while the mRNA expression of IL18 varied (Fig. 10a). The immunohistochemistry results were studied by using the HPA database in normal breast tissue and tumor tissue to explore the clinical significance of the signature. The results showed the expression and distribution of GPX4, GSDMD, GSDMC, IL18 in breast cancer and normal tissues (Fig. 10b).
Breast cancer is one of the most common tumors in women. At present, the treatment of breast cancer includes surgery, chemotherapy and radiotherapy. Although the cure rate has been greatly increased, so far, there is still no effective method to accurately predict the prognosis of patients with BC. Previous studies have found many tumor molecular markers of BC. The detection and targeted treatment of estrogen receptor, progesterone receptor and human epidermal growth factor receptor 2 have been widely used in clinical practice . But the lack of accurate prognostic biomarkers are still the main problem with improve clinical outcomes in patients with breast cancer. In recent years, with the correct understanding of pyroptosis, the latest research has found that pyroptosis plays an important role in the occurrence and development of tumors [47, 48]. Therefore, the study of biomarkers related to pyroptosis is expected to treat breast cancer more accurately.
In this study, 4 pyroptosis-related genes were found differently expressed in breast cancer and related to prognosis according to TCGA. So, we firstly establish a pyroptosis-related genes signature in the context of BC before evaluating its prognostic value and clinical significance. Then, patients with BC were classified according to the expression of pyroptosis-related genes. In order to obtain more samples and verify the feasibility of signature, two sets of GEO data were downloaded, normalized and integrated. Our gene signature was found to be able to predict prognosis in BC patients with high accuracy in training and testing cohorts. In addition, the risk score was identified as an excellent independent prognostic factor characterized by good sensitivity and specificity.
Then we further explore the relationship between pyroptosis and BC, The high-risk group was also rich in biological processes related to malignant progression, and there were significant differences between the two risk model subgroups of BC patients in both the training and testing cohorts. In the high-risk group, almost all immune cell infiltration and immune function were suppressed. Given the critical roles of these immune cells in stimulating anti-tumor immunity , it is reasonable to conclude that anti-tumor immunity was significantly reduced in the high-risk group of BC patients. In addition, the ESTIMATE algorithm showed that the stromal and immune cell scores were both inversely associated with risk scores, confirming poor immune cell infiltration in the high-risk subgroup. Triple negative breast cancer (TNBC) has a poor prognosis and high mortality compared to other breast cancers , so intensive efforts have been made to develop treatments targeting TNBC. Cancer immunotherapy targeting PD-L1 has improved outcomes for TNBC , and shown efficacy in other several cancers . PD-L1 and PD-L2 are key regulators of immune responses [52, 53]. This study also verified that PD-L1 and PD-L2 were significantly different in the two risk subgroups, and both were negatively correlated with risk score. The levels of nearly all immune checkpoints were significantly lower in the high-risk subgroup than in the low-risk subgroup, indicating that the immune response was dramatically altered in this group. Comprehensive analysis of immune cells, immune function, immune-related markers and PRGs confirmed the important role of pyroptosis in immune regulation in TME landscape. CSCs are considered as the major cause to tumor initiation, recurrence, metastasis, and drug resistance, driving poor clinical outcomes in patients . In this study, the risk signature was positively correlated with the stem cell score, confirming that our newly constructed gene signature was a risk factor for BC. Some researchers have been studying the role of pyroptosis-related genes in the development of cancer. And the downregulation of GSDMD was found to attenuate tumor proliferation via the intrinsic mitochondrial apoptotic pathway and inhibition of EGFR/Akt signaling and predicted a good prognosis in non-small cell lung cancer . However, it has also been proved that down-regulation of GSDMD promotes gastric cancer proliferation by regulating cell cycle-related proteins and over-expression of GSDMC is a prognostic factor for predicting a poor outcome in lung adenocarcinoma [20, 55]. So, more experiments should be done to confirm our findings.
Despite the prognostic value of the risk signature, there are several limitations in this study. First of all, this was a retrospective analysis, thus, prospective studies are needed to confirm the results. Secondly, there is a lack of experimental analysis to validate the results of bioinformatics analyses. In the future, more functional studies are needed to understand pyroptosis-related genes and their role in BC development.
In conclusion, 4 pyroptosis-related genes were found associated with BC prognosis in this study. The signature was proved to be independently associated with OS in TCGC cohort and GEO validation cohort. More significantly, it was found extremely valuable in functional analysis, tumor microenvironment, and drug sensitivity, providing insight for predicting the prognosis of BC. But the specific potential mechanism between pyroptosis-related genes and tumor immunity is still unclear and deserves further study. Our work will help shed light on the role of pyroptosis in tumgenesis, particularly in the areas of immune response, tumor microenvironment and drug resistance, which are crucial for the development of personalized cancer therapies.
Gong W, et al. STING-mediated syk signaling attenuates tumorigenesis of colitis‑associated colorectal cancer through enhancing intestinal epithelium pyroptosis. 2022;28:572-85.
Shen Z, et al. Metformin inhibits hepatocellular carcinoma development by inducing apoptosis and pyroptosis through regulating FOXO3. 2021;13.
Zhang Z, et al. Caspase-3-mediated GSDME induced Pyroptosis in breast cancer cells through the ROS/JNK signalling pathway. J Cell Mol Med. 2021;25:8159-68.
Gao J, et al. Downregulation of GSDMD attenuates tumor proliferation via the intrinsic mitochondrial apoptotic pathway and inhibition of EGFR/Akt signaling and predicts a good prognosis in non-small cell lung cancer. 2018;40(4):1971–84.
The Medical Ethics Committee of The Affiliated Hospital of Qingdao University approved this study and is in compliance with the Helsinki Declaration. Written consent was obtained from all participants.
Consent for publication
The authors declare that they have no competing interests
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Gong, M., Liu, X., Zhao, X. et al. A pyroptosis-related gene signature predicting survival and tumor immune microenvironment in breast cancer and validation.
BMC Cancer22, 1005 (2022). https://doi.org/10.1186/s12885-022-09856-y