- 1Department of Gastroenterology, Affiliated Yueqing Hospital, Wenzhou Medical University, Wenzhou, China
- 2Department of Laboratory Medicine, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, China
- 3Department of Hematopathology, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, China
- 4Department of Rheumatology, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, China
Triple-negative breast cancer (TNBC) is associated with poor prognosis and invalid therapeutical response to immunotherapy due to biological heterogeneity. There is an urgent need to screen for reliable indices, especially immunotherapy-associated biomarkers that can predict patient outcomes. Pyroptosis, as an inflammation-induced type of programmed cell death, is shown to create a tumor-suppressive environment and improve the chemotherapeutic response in multiple tumors. However, the specific therapeutic effect of pyroptosis in TNBC remains unclear. In this study, we present a consensus clustering by pyroptosis-related signatures of 119 patients with TNBC into two subtypes (clusterA and clusterB) with distinct immunological and prognostic characteristics. First, clusterB, associated with better outcomes, was characterized by a significantly higher pyroptosis-related signature expression, tumor microenvironment prognostic score, and upregulation of immunotherapy checkpoints. A total of 262 differentially expressed genes between the subtypes were further identified and the Ps-score was built using LASSO and COX regression analyses. The external GEO data set demonstrated that cohorts with low Ps-scores consistently had higher expression of pyroptosis-related signatures, immunocyte infiltration levels, and better prognosis. In addition, external immunotherapy and chemotherapy cohorts validated that patients with lower Ps-scores exhibited significant therapeutic response and clinical benefit. Combined with other clinical characteristics, we successfully constructed a nomogram to effectively predict the survival rate of patients with TNBC. Finally, using the scRNA-seq data sets, we validated the landscape of cellular subtypes of TNBC and successfully constructed an miRNA-Ps-score gene interaction network. These findings indicated that the systematic assessment of tumor pyroptosis and identification of Ps-scores has potential clinical implications and facilitates tailoring optimal immunotherapeutic strategies for TNBC.
Introduction
As the most common tumor with high morbidity and mortality in women, breast cancer (BC) has a poor prognosis and exacerbates a critical social burden worldwide (Global Burden of Disease Cancer Collaboration et al., 2017). Triple-negative breast cancer (TNBC), accounting for 10%–17% among all BCs, is a special subtype characterized by negative human epidermal growth factor receptor 2 (HER2), progesterone receptor (PR), and estrogen receptor (ER) (Lin et al., 2012). Due to the absence of the corresponding receptors, patients with TNBC fail to benefit from endocrine targeted therapy and HER2-targeted agents; hence, chemotherapy and surgery remain the most common treatment for patients with TNBC (Bergin and Loi, 2019). Recently, the rapid rise of immunotherapy with the combination of cisplatin or other platinum drugs, including anti-programmed cell death (PD)-1 and PD-ligand 1 (PD-L1) agents, has brought a new therapeutic landscape for patients with TNBC who did not benefit from conventional chemotherapy, radiation, or surgery (Nolan et al., 2017). However, in clinical practice, some patients with TNBC are still lacking an effective therapeutical response to immunotherapy due to genetic and biological heterogeneity (Vikas et al., 2018). Therefore, it is crucial to identify novel subtypes and screen reliable biomarkers (especially immunotherapy-related biomarkers) that can predict outcomes of patients with TNBC.
In the clinical setting, the TNM stage system is acknowledged as the most frequently used tool to predict the prognosis of patients with TNBC, which majorly depends on the inherent anatomical abnormity, including tumor size, lymph node situation, and distant metastatic status (Park et al., 2019). However, the occurrence of biological and tumor genetic heterogeneity makes it challenging for the TNM system to predict disease progression and prognosis (Park et al., 2019). Pyroptosis is a form of pro-inflammatory programmed cell death (PCD) that cleaves the gasdermin D (GSDMD) protein by classical or nonclassical pathways and triggers the production and release of cytokines (including inactive cytokines like IL-18 and IL-1β) to induce a strong inflammatory response (Yang et al., 2016). Pyroptosis is reported to create a tumor-suppressive environment by releasing inflammatory factors; therefore, inducing pyroptosis in tumors via chemotherapeutic drugs could produce antitumor effects (Shi et al., 2015). In vitro, Nathalia et al. demonstrate that omega-3 fatty acids can induce pyroptosis in TNBC cells via inducing the active CASP1 increase, further leading to the cleavage of GSDMD, formation of membrane pores, and the release of IL-1β (Pizato et al., 2018). However, the exact contribution of pyroptosis on the therapeutical response of immunotherapies and its role in the prognosis of TNBC remains unclear.
The classification of patients with TNBC based on transcriptome profiles via next-generation sequencing is considered a novel method to quickly indicate biological characteristics and help screen for the most appropriate treatment strategies (He et al., 2018). Besides conventional expression profiles, various biological signatures are also applied to identify novel molecular subtypes for the prognosis of TNBC, such as autophagy-related signatures (Kim et al., 2012), N6-methyladenosine (Wu et al., 2021), immune cell infiltration (Harano et al., 2018), etc. In this study, we aimed to build a novel scoring model (called Ps-scores) based on pyroptosis-related signatures to identify clustering subtypes of TNBC and correlate the characteristics of each subtype with prognosis, immunotherapy, and immune cell infiltration in patients with TNBC. Combining the Ps-scores and other classical clinical features, the predicted model was established to improve prognostic risk stratification and facilitate the decision making of treatments for patients with TNBC. Moreover, using single-cell RNA sequencing (scRNA-seq) technology, we successfully validated the potential cellular subtypes of TNBC and expounded the predominant expression characteristics of Ps-score-related genes in each cluster. Finally, the targeted miRNAs were predicted by combining multiple databases with differentially expressed miRNAs (DEmiRNAs), whereas the miRNA-Ps-score signature interaction network was further constructed to visualize the potential regulatory relationship. These results imply the potential links among the pyroptosis-related scores, immune microenvironment, prognosis, and response to immunotherapy for patients with TNBC. Our findings provide new insight into the prognostic signatures of TNBC and will help develop promising strategies for TNBC immunotherapy.
Materials and Methods
TNBC Data set Preparation and Preprocessing
Transcriptome profiling data (FPKM value) of 1217 BC samples with their corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA) data sets (https://portal.gdc.cancer.gov/). Through screening the “negative HER2, PR, and ER” status based on clinical data, we finally identified 119 patients with TNBC and comprehensive clinical information. Other microarray data sets of 819 patients with TNBC and prognostic information were also downloaded from the Gene Expression Omnibus (GEO) data sets (https://www.ncbi.nlm.nih.gov/geo/), including 120 ER-negative BC in GSE16446, 198 HER2-negative BC in GSE25065, 310 HER2-negative BC in GSE25055, 107 TNBC in GSE58812, and 84 TNBC in GSE157284. In addition, corresponding miRNA sequencing data sets and mutation files were also obtained from the TCGA-BRCA to investigate the miRNA regulatory mechanism, and transcriptome profiling of 179 normal breast tissues was obtained from the Genotype-Tissue Expression (GTEx) database as normal controls (Carithers and Moore, 2015). The “ComBat” algorithm of the “sva” package was further applied to remove the nonbiological technical biases due to batch effects between different data sets (Leek et al., 2012). To remove the false positives caused by batch effects, we selected several stable internal reference genes (HPRT1, PPIA, RPS13, TBP, GAPDH, and HMBS) to perform the PCA analysis, which are reported as valid reference genes for human BC cell lines by Liu et al. (2015). Moreover, the IMvigor210 data sets (Hoffman-Censits et al., 2016), including 316 metastatic urothelial carcinomas (mUCs) with immunotherapy, were applied to investigate the therapeutic reaction, and the scRNA-seq data of 1534 cells from six patients with TNBC (GSE 118389) were used to validate molecular subtypes for TNBC. Detailed information on these data sets is listed in Supplementary Table S1.
Identification of Pyroptosis-Related Signatures
According to previous studies, the caspase family, especially caspase1/4/5/8 (CASP1/4/5/8) was reported to specifically cleave GSDMD to further activate pyroptosis (Shi et al., 2015; Orning et al., 2018). In addition, Zhang et al. (2020), also found that CASP3 and granzyme B (GZMB) could convert cell apoptosis into pyroptosis through cleaving gasdermin E (GSDME). Granzyme A (GZMA) was also considered to be essential in inducing cell pyroptosis by cleaving gasdermin B (GSDMB) (Zhou et al., 2020). Moreover, inflammasome-associated families, such as absent in melanoma 2 (AIM2) and nucleotide-binding domain and leucine-rich repeat receptor (NLR), are demonstrated to induce the pyroptosis process through activating CASP1 and the release of IL1β and IL18 (Man and Kanneganti, 2015; Karki and Kanneganti, 2019). Therefore, based on the published research, a total of 33 pyroptosis-related signatures were chosen, and 24 genes were retained for subsequent analysis after filtering out the signatures with low expression (sum FPKM value of all samples less than one).
Consensus Cluster Analysis for Pyroptosis-Related Signatures in TNBC
Based on the expression of pyroptosis-related signatures, we performed hierarchical clustering analysis and applied the “ConsensusClusterPlus” R package (Wilkerson and Hayes, 2010) to conduct unsupervised clustering based on Euclidean distance and Ward’s linkage methods 1000 repeated times to ensure the classification stability. During the process, the clusters from 2 to 9 were performed, respectively, and the optimal clustering model was determined based on the consensus cumulative distribution function (CDF) plot. Moreover, we performed multiple comparisons among different pyroptosis-subtypes, including for the tumor microenvironment (TME), prognosis, and vital clinical-pathological phenotypes to explore their characteristics. The R packages “survival” (Therneau and Lumley, 2015) and “survminer” (Kassambara et al., 2017) were used to perform Kaplan–Meier survival analysis and draw survival curves between pyroptosis subtypes.
Identification of Differentially Expressed Genes (DEGs) and Functional Enrichment Analysis
To identify the DEGs between pyroptosis subtypes, the empirical Bayesian algorithm was applied through the “Limma” R package (Smyth, 2005), and the significance cutoff was set as adjusted p < .05 and absolute fold-change >1. To clarify the biological function and characteristics of pyroptosis clusters, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed by using the “ClusterProfiler” R package (Yu et al., 2012), and the results were visualized using the “ClueGO” plugin in Cytoscape v3.7.1 (Bindea et al., 2009).
TME Cell Infiltration and Gene Set Variation Analysis
To evaluate the immune cell infiltration (ICI) characteristics of TNBC tissues, we used the “CIBERSORT” R package (Chen et al., 2018) to quantitatively analyze the infiltration levels of different immune cells with the LM22 signatures by 1000 random permutations. The tumor purity scores, ICI levels, and stromal contents in different samples were evaluated via the “ESTIMATE” algorithm (Yoshihara et al., 2013). Moreover, through the “c2. cp.kegg.v6.2. symbols” data sets downloaded from the MSigDB database, we performed GSVA using the “GSVA” R package and drew a heatmap to exhibit the different immunogenic pathways (Hanzelmann et al., 2013).
Definition of Immune Characteristics Between High and Low Ps-Score Groups
To further identify a novel index representing the characteristics of the pyroptosis subtypes, we conducted univariate Cox proportional hazards regression analysis for overall survival (OS) to preliminarily screen significant genes through using the “coxph” function in the “survival” R package. Subsequently, to remove the multicollinearity among these candidate genes, LASSO regression was applied to screen independent prognosis-related genes with the optimal penalty parameter and a minimum 10-fold cross-validation (Ranstam and Cook, 2018). After further adjustment, the multivariate Cox regression (stepwise model) was conducted to identify hub genes, and the coefficients obtained from the regression algorithm were used to acquire the Ps-score based on the following formula:
Construction and Evaluation of the Pyroptosis-Related Prediction Model
The multivariate Cox regression (stepwise model) was applied to construct the prognostic model for TNBC-combined Ps-scores and other clinical features, including age, clustering subtypes, clinical stages, and TNM stages. Variables with p-values < .05 were included in the Cox regression model, and the nomogram was further constructed to predict the probability of 1-, 3-, and 5-year survival in patients with TNBC using the “survival” R package. To evaluate and validate the prediction capability of the nomogram, we calculated the concordance index and plotted the calibration curves for 3- and 5-year survival through a bootstrapping method with 1000 resamples. To further investigate the expression of the Ps-score-related genes at the protein level, the Human Protein Atlas (HPA) (Ponten et al., 2008) was used to display the results of the immunohistochemistry (IHC) technique. The detailed information of patients is included in Supplementary Table S17.
Exploration of the Significance of Ps-Scores in Response to Clinical Immunotherapy
Based on the IMvigor210 data sets with atezolizumab treatment, we performed a comprehensive comparison between different Ps-score subgroups, including response to immunotherapy, immune phenotype, and clinical remission rate. Moreover, to evaluate the potential therapeutic value of Ps-scores in chemotherapy for TNBC, we calculated the half-maximal inhibitory concentration (IC50) of common chemotherapeutic drugs based on the Genomics of Drug Sensitivity in Cancer (GDSC) databases (Yang et al., 2013). Antitumor drugs such as 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel are recommended for BC treatment by current clinical guidelines. Differences in IC50 of these chemotherapeutic drugs between Ps-score subgroups were compared by Wilcoxon test with the results exhibited in box diagrams using the “ggpubr” R package (Whitehead et al., 2019).
Validation of Molecular Subtypes Based on scRNA-Seq Analysis
To validate molecular clusters and further seek biomarkers of each cluster, the Seurat pipeline was selected for subsequent analysis. Using the Seurat package v3.0 (Butler et al., 2018), we transformed the data matrix into a “Seurat object” through the “CreateSeuratObject” R function and performed the necessary quality control. The violin diagram exhibited the number of sequencing reads per sample and the expression of mitochondrial genes (Supplementary Figure S2A). Further, to remove the influence from mitochondrial and extreme genes, we kept the number of sequenced genes at 200–10,000, directly including the majority genes, and removed the cells with average gene expression <10 and mitochondrial genes >5%. Then, we conducted data standardization through the “NormalizeData” function with the method of “LogNormalize” and used the top 1500 variable counts to perform PCA using the “RunPCA” function (Supplementary Figures S2B,C). Subsequently, t-distributed statistical neighbor embedding (tSNE) was applied to visualize the density clustering, and the “SingleR” package was applied for cell-subtype annotation based on the marker genes (Aran et al., 2019). Moreover, the pseudotime trajectory analysis was further performed using the “Monocle” package v2.0 to expound the potential inner relationship among these cell clusters (Trapnell et al., 2014).
Prediction of Potential miRNA Targets for Prognosis-Associated Signatures
The miRNAs targeting Ps-score-related genes were predicted based on the following databases: TargetScan (http://www.targetscan.org/), starBase (http://starbase.sysu.edu.cn/starbase2/index.php), miRTar (https://mirtarbase.cuhk.edu.cn/), and miRDB (http://www.mirdb.org/). The expression of TNBC miRNA was downloaded from the TCGA-BRCA data sets, and the DEmiRNAs were further identified by the “Limma” R package. Subsequently, we identified the intersection of predicted miRNA by four databases and DEmiRNAs as regulated miRNAs for each hub gene and further visualized the miRNA–mRNA interaction network using Cytoscape v3.7.1.
Results
Overview of Genetic and Biological Characteristics of Pyroptosis-Related Signatures in TNBC
After a series of rigorous screening and quality control steps, a total of 24 pyroptosis-related signatures remained for subsequent analysis in our study (Supplementary Table S2). Combined with the normal tissues in GTEx data sets, we first compared the expression of pyroptosis-related signatures between patients with TNBC and normal controls. We found that most pyroptosis-related genes were significantly upregulated in TNBC groups, including CASP1/3/5/8, GSDMA/C, NLRC4/P3/P7, IL18, IL1β, and TNF (Figure 1A). PCA indicated that the expression of these pyroptosis-related signatures could be used to divide the TNBC samples and controls into two distinct clusters (Figure 1B). In addition, the PCA of internal reference genes revealed that nonsignature genes failed to discriminate TNBC and control cohorts, indicating that the separation created a true distinction based on the pyroptosis-related DEGs in patients with TNBC (Supplementary Figure S1A). Moreover, the KEGG functional enrichment analysis revealed that these pyroptosis-related genes were predominantly focused on infectious diseases, immune response, and cellular signal conditioning mechanisms, including the NLR signaling pathway, p53 signaling pathway, TNF signaling pathway, and apoptosis (Figure 1C, Supplementary Table S3). In terms of genetics, 148 of the 203 samples (72.91%) manifested pyroptosis-related signatures in mutations and the NLR families, especially NLRP3 and NLRP7, exhibiting the highest frequency of mutations (Figure 1D, Supplementary Table S4). Moreover, the top 10 pyroptosis-related genes with the most frequent mutations were located on the 24 human chromosomes (Figure 1E).
FIGURE 1. Characteristics of pyroptosis-related signatures in patients with TNBC. (A) The expression of pyroptosis-related signatures between normal tissues and TNBC tissues from GTEx and TCGA data sets; Wilcox test, *p < .05, **p < .01, ***p < .001; ns, not statistically significant. (B) PCA showing pyroptosis-related signatures sorted TNBC and control tissues into two clusters. (C) The KEGG functional analysis of pyroptosis-related signatures. (D) The landscape of mutation profiles in patients with breast cancer from TCGA-BRCA cohort. (E) The location of the top 10 pyroptosis-related signatures with the most frequent mutations on the 24 human chromosomes.
Identification of a TNBC Cluster Pattern Based on Pyroptosis-Related Signatures
Based on the expression of pyroptosis-related signatures, we used an unsupervised clustering method to identify the subtypes of patients with TNBC and identified k = 2 as the optimum clustering model from k = 2 to k = 9 clustering with the least area under the consensus CDF curve for 69 patients in clusterA and 50 patients in clusterB (Figure 2A, Supplementary Table S4; Supplementary Figure S1B). To further clarify the intrapatient heterogeneity of patients with TNBC, we performed a comparison in the clinical differences between subtypes and found that patients in clusterB were negatively associated with severe clinical stages, including the pathological and TMN stages. Furthermore, the survival analysis showed that patients in clusterB had a longer median survival time than those in clusterA with more patients in clusterB also receiving radioactive treatments, indicating that patients with TNBC in clusterB might have a better prognosis (Figure 2B). There was no significant difference in other clinical indexes including age, sex, M stage, and the ratio of pharmaceutical and surgical therapies (Supplementary Figure S1C). Notably, the expression of pyroptosis-related gene signatures was significantly increased in patients in clusterB compared with that of the clusterA cohort (Figure 2C). In terms of the immune infiltration scores, adaptive immune response-related lymphocytes (including memory B cells, activated memory CD4+ T cells, plasma cells, CD8+ T cells, and gamma delta T cells) were significantly increased in patients in clusterB compared with the clusterA cohort. However, innate immunity and immunoregulation-related cells were significantly infiltrated in clusterA cohorts, including neutrophils, activated natural killer cells, resting memory CD4+ T cells, and regulatory T cells (Tregs) (Figure 2D, Supplementary Table S5). Higher stromal scores and immune scores with lower tumor purity were also detected in patients in clusterB compared with the clusterA groups (Figures 3A,B).
FIGURE 2. Identification of pyroptosis-related subtypes in TNBC. (A) Consensus clustering matrix for k = 2 in patients with TNBC. (B) Kaplan–Meier curves of OS for patients with TNBC divided into two subtypes. ClusterB was negatively associated with severe clinical stages, including pathological and TMN stages. (C) Heatmap showing the expression of pyroptosis-related signatures upregulated in ClusterB subtypes. (D) Boxplots show the difference in immune cell infiltration between ClusterA and ClusterB. Wilcox test, *p < .05, **p < .01, ***p < .001.
FIGURE 3. Identification of DEGs based on pyroptosis-related clusters. (A–B) Box plot showing higher stromal and immune scores with lower tumor purity detected in patients of clusterB than the clusterA group. (C) Volcano plots displaying the up- and downregulated DEGs between two subgroups in TNBC cohorts. (D) Bubble diagram showing the results of KEGG enrichment analysis of the subtypes. (E) LASSO coefficient profiles of 12 prognostic related genes and 10-times cross-validation for tuning parameter selection in the LASSO model. (F) Forest map displaying the HR and p-value of six hub genes after multivariate Cox regression analysis. (G) Heatmap showing the distinct expression of six hub genes between pyroptosis-related clusters.
Identification of DEGs Based on Pyroptosis-Related Clusters
Considering the biological characteristics of immune subtypes in TNBC, we conducted a DEG analysis between the two subtypes. Through comparing clusterA with clusterB groups, a total of 262 DEGs (including 13 clusterA- and 249 clusterB-related genes) in TNBC were identified (Figure 3C, Supplementary Table S6). To further interpret biological processes and pathways of pyroptosis-related subtypes, these DEGs were chosen to perform KEGG functional analysis. The results showed that clusterA-related genes were not enriched in any significant pathways while the clusterB-related signatures were predominantly enriched in immune activation-associated pathways, including natural killer cell–mediated cytotoxicity, the toll-like receptor signaling pathway, chemokine signaling pathway, cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, TNF signaling pathway, Th17 cell differentiation, and Th1 and Th2 cell differentiation (Figure 3D, Supplementary Table S7).
Development of Ps-Score and Characteristic of Ps-score-related Subgroups
After successively including the 262 DEGs in univariate Cox regression, LASSO regression, and multivariate Cox regression analysis as candidate prognosis-associated genes, we identified six hub genes (including CFB, IFITM1, EPSTI1, MARCO, CXCL13, and CCL5) from the Ps-score signatures based on their β coefficients (Figures 3E,F, Supplementary Table S8). In addition, the expression of these hub genes was higher in clusterB subgroups. Based on the IHC data from the HPA database, the expression of these hub genes at the protein level was further validated in BC, especially CFB, ESPIT1, and IFITM1 (Figure 3G, Figure 6C). Based on the expression of these genes and their corresponding β coefficients, the Ps-score was defined by the following formula:
FIGURE 4. Characteristics of Ps-scores and correlation to pyroptosis in TNBC. (A–B) Kaplan–Meier curves of overall survival (OS) for the patients with TNBC in the high- and low-Ps-score groups in the TCGA (A) or GEO data sets (B). (C, E) The box plots show that ClusterB possessed a lower level of Ps-scores associated with pharmaceutical and surgical therapies, and lower pathological stages. (D) Box diagram displaying the expression of pyroptosis-related signatures significantly increased in the low Ps-score groups. (F) Heatmap used to visualize pathways analyzed by GSVA showing the active biological pathways in distinct pyroptosis-related clusters.
To investigate the biological characteristics of the Ps-scores, we compared the expression of pyroptosis-related genes and ICI between different Ps-score groups and further explored the correlation between significant clinical phenotypes and the Ps-scores. The results reveal that clusterB possessed a lower level of Ps-scores associated with pharmaceutical and surgical therapy as well as lower pathological stages in patients with TNBC (Figure 4C&E). Interestingly, the expression of pyroptosis-related signatures was significantly increased in the low Ps-score groups, including the CASP, GSDM, and NOD families, as well as inflammatory factors, suggesting the potential role of pyroptosis activation in the low-score of TNBC cohorts with better prognosis (Figure 4D). Moreover, ROC analysis showed 1-, 3-, and 5-year AUC values of the Ps-scores for predicting the prognosis of patients with TNBC of 0.867, 0.867, and 0.906, respectively, in the TCGA sets (Figure 5A). Furthermore, immune infiltration analysis revealed that substantial immune cells were significantly inhibited in the high Ps-score groups, including that of CD8+ T cells, follicular helper T cells, activated CD4+ memory T cells, and plasma cells (Figure 5D). The correlation analysis also indicates that the Ps-scores were significantly positively associated with the levels of tumor purity (R = 0.57, p < .001) but negatively associated with stromal and immune scores (R = –0.48, p < .001 and R = –0.50, p < .001, respectively) (Figure 5E).
FIGURE 5. Establishment and evaluation of the Ps-score-related prognostic model for patients with TNBC. (A) Time-dependent receiver operating curves of 1-, 3-, and 5-year survival for patients with TNBC using Ps-scores. (B) Combined nomogram for predicting the probability of 1-, 3-, and 5-year survival for patients with TNBC. (C) Calibration curve of the established nomogram with 3- and 5-year survival, respectively. (D) Immune infiltration analysis revealed that substantial immune cells were significantly inhibited in high Ps-score groups. (E) Correlation analysis shows the Ps-scores significantly positively associated with the levels of tumor purity (R = 0.57) and negatively associated with stromal and immune scores (R = –0.48 and R = –0.50, respectively).
Based on the Ps-scores and some primary clinical characteristics, multivariate Cox regression analysis was conducted to construct a nomogram that could accurately predict the probability of the 1-, 3-, and 5-year survival for patients with TNBC. The Ps-scores, age, pyroptosis-related cluster, and clinical stages were considered as related predictors for the prognosis of patients with TNBC and incorporated into the nomogram with significant regression coefficients and p-values (Figure 5B, Supplementary Table S10). From the nomogram, we could observe that the Ps-score contributed the most to the total score with a 0.716 concordance index (Figure 5B, Supplementary Table S11). Calibration curves exhibited that the nomogram had a good prediction capacity in both 3- and 5-year OS for patients with TNBC (Figure 5C).
Significance of Ps-Scores in the Prediction of Response to Immunotherapy and Common Chemotherapeutics
The alluvial diagram visualized the status changes in the different characteristics of patients (Figure 6A). We found that patients with a low Ps-score in clusterB subtypes had a higher ratio of survival status. When using the TCGA and other external GEO data sets, the results of GSVA demonstrated the coincident negative correlation between Ps-scores and immunoregulation-related pathways, such as the toll-like receptor signaling pathway, antigen processing and presentation, rig I-like receptor signaling pathway, T cell receptor signaling pathway, NOD-like receptor signaling pathway, and JAK-STAT signaling pathway (Figure 6B, Supplementary Table S12).
FIGURE 6. The significance of Ps-scores in the prediction of response to immunotherapy and common chemotherapeutics for TNBC. (A) Alluvial diagram visualizing the status changes from different characteristics of patients. (B) Correlation analysis of Ps-scores and GSVA pathways in TCGA and GEO data sets. The size of nodes represents the correlation coefficient, and the red or green represents positive or negative correlation, respectively. (C) The immunohistochemical results of Ps-score-related genes from the HPA database. (D) The box plots displayed patients with better responses to PD-L1 treatment exhibited higher Ps-scores using GSE157284. (E) The relative expression of immune-checkpoint genes was significantly increased in the high Ps-scores groups of GSE157284 data sets. (F) The relative expression of immune-checkpoint genes was significantly increased in the high Ps-scores groups of IMvigor210 cohorts. (G) The box plots indicated patients with better responses to PD-L1 treatment exhibited higher Ps-scores using IMvigor210 cohorts. (H–I) Ps-score in different ACRG subtypes and the rate of CR after immunotherapy in IMvigor210 cohorts. (J) Kaplan-Meier curves of OS for the patients with TNBC in high and low Ps-score groups of IMvigor210 data sets. (K) Comparison of IC50 value of 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel in high and low Ps-score groups using GDSC databases.
To further explore the role of Ps-scores in predicting the therapeutic benefit in TNBC, we first calculated the Ps-scores of patients who accepted anti-PD-L1 immunotherapy from the GSE157284 and IMvigor210 cohorts before assigning them into high or low Ps-score groups. From the GSE157284 data set, patients with an effective response rate to anti-PD-L1 therapy showed lower Ps-scores, whereas the relative expression of immune-checkpoint genes was significantly increased in the high Ps-scores groups (Figures 6D,E, Supplementary Table S13). Congruously, we also found that the low Ps-score group had higher expression of immune-checkpoint genes in the IMvigor210 cohort and effective anti-PD-L1 responders also exhibited lower Ps-scores (Figures 6F,G). In metastatic urothelial cancers of the IMvigor210 data sets, distinct immunological subtypes might result in opposite therapeutic responses. Therefore, we also compared the Ps-score levels among these subtypes and found that the lowest Ps-score was in the inflamed subtypes (Figure 6H). Moreover, the rate of complete remission (CR) after immunotherapy was also increased in the low Ps-score cohort compared with the high Ps-score cohort, with the low Ps-score group validated to have a better prognosis for metastatic urothelial cancers (Figures 6I,J). All these results suggested that the Ps-score might serve as a significant indicator in immunotherapy decision making for cancers.
Besides checkpoint blocker therapy, we also investigated the potential associations between the Ps-scores and the curative efficacy of common chemotherapeutics in treating BCs. The IC50 value was calculated for five common anti-BC chemotherapeutic drugs obtained from the GDSC databases, including 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel (Supplementary Table S14). Notably, all the drugs exhibited lower IC50 values in the low Ps-score groups, indicating patients with low Ps-scores might obtain a better curative efficacy from common chemotherapy (Figure 6K). Collectively, these outcomes indicate that Ps-scores could be associated with the response to immunotherapy and common chemotherapy.
Validation of Cellular Subtypes in TNBC Through scRNA-Seq Analysis
To validate the potential subtypes of patients with TNBC, the GSE118389 data set along with 1534 cells from six TNBC tissues were used to identify concrete cellular subtypes and corresponding marker genes. A total of 12 distinct cellular clusters were identified through tSNE analysis (Figure 7A) with the marker-genes of each cluster listed in Supplementary Table S15. Moreover, the results of cell-type annotation using the “SingleR” package indicated these cell clusters fell into six cellular subtypes, including epithelial cells, erythrocytes, CD8+ T cells, fibroblasts, endothelial cells, and monocytes, of which the epithelial cells were the most common cell type with six subtypes (Figure 7B). Notably, of the six Ps-score-related genes, five genes were identified as marker genes, and their expression in each cellular subtype is shown in a bubble diagram (Figure 7C). The expression of CCL5 and IFITM1 were generally increased in nearly all cellular subtypes and significantly higher than the expression of other signatures. Although the expression of remanent genes was relatively lower in these cells, significant cellular specificity was found in these Ps-score-related genes. For example, MARCO was particularly expressed in the epithelial cell subtype 4 and monocytes while CFB was particularly expressed in the fibroblast subtype 2 as well as epithelial cell subtypes 3 and 4. Interestingly, EPSTI1 was significantly expressed in immune-related cellular subtypes including monocytes and CD4+ T cells, consistent with the results of IHC. In addition, the pseudotime trajectory analysis also revealed the arrangement of different cellular subtypes that formed a certain rule based on its spatial relationships (Supplementary Figure S2E). Concretely, the trajectory analysis of epithelial cells revealed that a small quantity of epithelial cell subtype 3 was distributed at the start of the trajectory while mixed cells from epithelial cell subtypes 5 and 6 were distributed at the end of the trajectory. Moreover, epithelial cell subtype 2 was uniformly located behind epithelial cell subtypes 1, whereas epithelial cell subtype 4 nearly existed throughout the trajectory (Figure 7D).
FIGURE 7. Validation of cellular subtypes in TNBC through scRNA-seq analysis. (A) Visualization of tSNE colored according to cell types for TNBC single-cell transcriptomes. (B) Heatmap revealing the scaled expression of top 10 marker-genes for each cluster defined in (A). (C) Dot plots showing the expression of indicated Ps-score-related genes for each cell cluster. (D) Monocle pseudospace trajectory analysis revealing the cellular lineage progression of Epithelial cells subtypes in patients with TNBC colored according to different Epithelial cell clusters. (E) The miRNA-hub genes interaction network based on the prediction of four databases. The round node represented the miRNA and the rhomb represented the Ps-score-related genes. (F) Volcano plots displayed the upregulated and downregulated DEmiRNAs between two subgroups in patients with TNBC. (G) Venn diagram of the predicted miRNA and DEmiRNAs.
Prediction and Validation of miRNAs Interacted With Hub Ps-Score Genes
To further explore the potential regulatory role of miRNAs targeted to these Ps-score-related genes in TNBC, we found 223 probable miRNAs and successfully constructed the miRNA-hub gene interaction network based on the prediction of the TargetScan, starBase, miRTar, and miRDB databases (Figure 7E, Supplementary Table S16). In addition, to validate the regulatory role of miRNAs in TNBC, we also identified 37 DEmiRNAs between pyroptosis-related clusters, before screening the top 20 miRNAs serving as vital regulatory factors (Figures 7F,G). Finally, we successfully simplified the interaction network with 20 miRNAs and four targeted genes to verify that ESPIT1 and CXCL13 are the most active targets regulated by massive DEmiRNAs (Figure 7H).
Discussion
As a malignant tumor with high mortality, TNBC is known for its poor prognosis, which stems from ineffective therapeutic response to immunotherapy due to tumor biological heterogeneity. Recently, the IMpassion130 trial demonstrated that the combination of atezolizumab, a PD-L1 inhibitor, and nab-paclitaxel could prolong the OS in PD-L1- patients with positive advanced TNBC, heralding the emergence of immunotherapy as an effective treatment strategy for TNBC (Schmid et al., 2018). In addition, The U.S. Food and Drug Administration (FDA) also approved SP142, a PD-L1 IHC assay, as an auxiliary test to identify eligibility for atezolizumab therapy in patients with advanced TNBC. However, the results of IHC staining are still short of high interlaboratory reproducibility with subjective judgment. For example, the IHC levels of PD-L1 were investigated in a total of 443 patients with BC by (Wang et al., 2017) but only ∼16% of these tumors exhibited positive PD-L1 levels. Although the PD-L1 test plays a potentially significant role in the management of multiple advanced carcinomas, objective standardization for this test has not been achieved; hence, its current use in the clinical practice poses a twofold risk to patients: false positive could result in potentially toxic therapies resulting in unforeseen complications, such as miscarriages, or PD-L1 false negatives would benefit from therapy but are excluded from receiving treatment (Reisenbichler et al., 2020). Therefore, identification of novel immunotherapy-related subtyping and reliable objective prognostic indicators for immunotherapy in TNBC is urgently needed.
In contrast to apoptosis, pyroptosis usually occurs in abnormal cells infected by microbes as a positive programmed cell death process, thus inducing the release of pro-inflammatory cytokines and activating an inflammatory response (Bedoui et al., 2020). Prompted by microbes, pyroptosis can also be converted from apoptosis and play various roles in multiple tumors. Pyroptosis has shown antitumor effects through inhibiting the tumor growth in liver and gastric cancers while showing both suppression or promotion effects in BC (Zaki et al., 2010; Chen et al., 2012; Shao et al., 2021). Shi et al. (2015) demonstrated that the activation of the NLRP3 (NOD-, LRR-, and pyrin domain-containing 3) inflammasome was integral for the activation of pyroptosis by recruiting CASP1, further leading to cleavage of GSDMD. In the present study, we explored all the signatures and pathways directly related to pyroptosis in TNBC and detected that the NLRP3/CASP1/GSDMD pathway-related pyroptosis was activated in patients with TNBC, implying that pyroptosis might participate in the mechanism of TNBC, which was associated with the prognosis of TNBC.
The classification of patients based on pathognomonic gene expression profiles is considered a proven method and applied to various studies of TNBC, including autophagy-related signatures (Kim et al., 2012), N6-methyladenosine (Wu et al., 2021), and immune cell infiltration (Harano et al., 2018). In this study, we first proposed a pyroptosis-related molecular subtype based on clustering pyroptosis-related signatures with distinct clinical and immunological characteristics. Interestingly, the characteristics of the two molecular subtypes manifested in significant homogeneity. We detected that patients in ClusterB presented a longer median survival time than those in ClusterA, whereas patients in ClusterB also negatively associated with serious clinical stages, including pathological and TMN stages, suggesting these pyroptosis-related signatures were also significantly associated with different survival risks in patients with TNBC. Our results also reached some consensus: (1) nearly all the pyroptosis-related signatures exhibited higher expression in patients in ClusterB; (2) ClusterB was a specific subtype with a better prognosis and slighter clinical pathological phenotypes; (3) ClusterB was identified as an immune-activated phenotype with higher TME immune scores and infiltration levels of adaptive immune response-related immune cells.
To further explore the potential biological functional features of the pyroptosis-related subtypes in TNBC, we investigated the DEGs between the two subtypes and performed KEGG function enrichment analysis. Consistent with the immunological signatures of subtypes, functional enrichment analysis revealed that immune-activation associated pathways, including natural killer cell-mediated cytotoxicity, the toll-like receptor signaling pathway, chemokine signaling pathway, cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, TNF signaling pathway, and Th17 cell differentiation as well as Th1 and Th2 cell differentiation were significantly enriched in the ClusterB cohorts. Of these pathways, the activation of multiple immune pathways is reported to suppress metastatic spread in TNBC (Zanker et al., 2020) and could be the potential mechanism for a better prognosis of patients in Cluster B.
Furthermore, to increase the clinical application value and create better clinical practicability, we successfully constructed a novel pyroptosis-related scoring tool (Ps-score) to determine the prognostic risk of TNBC based on six hub genes from two clusters. High-expression of these risk signatures at the protein levels was confirmed by IHC from the HPA database, and the Ps-scores effectively stratified patients with TNBC from the TCGA and GEO data set, respectively, into high- and low-risk groups. Survival analysis revealed that the low-score groups had longer OS than patients with high Ps-scores, and ROC curves exhibited a great predictive capacity of Ps-scores for the 1-, 3-, and 5-year survival of TNBC. In addition, the Ps-scores were significantly decreased in the ClusterB cohorts, indicating that the Ps-score could reflect the heterogeneity of patients with TNBC. Moreover, the Ps-score also represented patients with different clinical outcomes and was associated with the response to immunotherapy. The patients with high Ps-scores exhibited worse clinical prognosis and lower expression of immune checkpoints. TME ICI analysis also demonstrated that the Ps-score was significantly negatively correlated with the infiltration levels and tumor purity, suggesting its value in immunotherapy. Finally, a correlation analysis of the Ps-score and pathways based on GSVA results demonstrate a coincident negative relation in both TCGA and other GEO data sets, implying the activation of the various immune-related processes might occur more frequently in low Ps-score cohorts. Notably, besides immune-related pathways activated in low-score patients, massive activation of signal regulatory pathways was observed in high Ps-score groups, including TGF-β1, Wnt, Notch, and the MAPK signaling pathway, which is reportedly involved in the mechanism of TNBC and as target pathways for drug treatments (Giltnane and Balko, 2014; Kim et al., 2016; Pohl et al., 2017; Giuli et al., 2019).
Furthermore, using the IMvigor210 data sets, we also speculated that the Ps-score might be applicable to estimate the clinical response to immunotherapy in other tumors as well. Besides immunotherapy, common chemotherapeutic drugs also demonstrated lower IC50 values in the low Ps-score cohorts, including 5-fluorouracil, cisplatin, docetaxel, doxorubicin, and paclitaxel from the GDSC database, implying that these chemotherapeutic drugs would be more effective in patients with TNBC with low Ps-scores. Overall, these findings from external data sets validated the potential benefits of using the Ps-score system and indicated its role in predicting curative responses to common chemotherapies and immune checkpoint therapies.
Finally, the scRNA-seq analysis demonstrated the authentic existence of cellular subtypes with their marker genes in patients with TNBC and clearly showed the distribution of Ps-score-related genes in each subtype. Admittingly, TNBC originated from epithelial cells, and the results of the scRNA-seq also demonstrated multiple subtypes of epithelial cells, reflecting different clusters of tumor cells. Interferon inducible transmembrane 1 (IFITM1) is reported to promote the progression of TNBC through regulating integrin, NFκB, and IL6 gene expression and might serve as a novel therapeutic target for patients with IFITM1+ TNBC (Provance et al., 2021). Of the Ps-score-related genes, our analysis also detected that IFITM1 exhibited relatively high expression in epithelial cell subtypes 3–5, consistent with the above patients with IFITM1+ TNBC. In previous studies, the epithelial-stromal interaction 1 (EPSTI1) is also shown to modulate the extrinsic apoptotic pathway in TNBC cell lines, which highlighted its potential as a therapeutic target for patients with TNBC (Capdevila-Busquets et al., 2015). Interestingly, EPSTI1 is overexpressed in monocytes and CD8+ T cells, suggesting EPSTI1 might participate in the process of extrinsic apoptosis with the activation of the immune response. Moreover, the pseudotime trajectory analysis displayed the distribution of tumor epithelial subtypes and demonstrated the existence of inner heterogeneity and potential cellular differentiation in patients with TNBC. For the common subtypes of BCs, microRNA profiles from different breast cells were applied to distinguish and reflect different subtypes, including luminal A, luminal B, and basal and malignant myoepithelioma, indicating that the expression of genes in cells could directly reflect the different subtypes in BCs (Bockmeyer et al., 2011). Despite the differences in cellular and individual subtypes, pyroptosis-related signature genes distinguished both subtypes of TNBC based on their differential expression. Combined with the differential expression of miRNAs between pyroptosis-related clusters, we ultimately constructed a miRNA–mRNA interaction network, including 20 miRNAs and four hub genes and found that EPSTI1 and CXCL13 were the central nodes with the most miRNA regulation.
Our study has the limitation that the high-throughput sequencing data sets for initial analysis were relatively insufficient as it was simply obtained from public databases. The corresponding results and conclusions remain to be further investigated through more external congeneric research and should be validated via functional experiments in vivo and in vitro. Furthermore, several conclusions of this study require further research to confirm its reproducibility, improve the clinical application of pyroptosis-related clusters, and elaborate on the role of Ps-scores in predicting the response to immunotherapy for TNBC.
Conclusion
Our study is the first to propose molecular subtypes based on clustering pyroptosis-related signature expression with distinct clinical and immunological signatures in patients with TNBC. Moreover, we identified and validated a Ps-score system as an effective tool to predict the OS and immunotherapy efficacy in patients with TNBC. Finally, we preliminarily explored the cellular subtypes using scRNA-seq data sets to demonstrate the heterogeneity of TNBC and successfully construct an interaction network to expound the regulatory miRNA targeted Ps-score-related signatures. The various transcriptomic analyses facilitated the screening of significant genetic signatures of TNBC to provide a new clinical application of Ps-scores in predicting prognosis and chemo-immunotherapeutic response for patients with TNBC.
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 authors.
Author Contributions
CL and JP contributed to data acquisition, analysis and drafting of the manuscript. YJ and YW participated in the process of data acquisition and analysis. ZJ contributed to figures presentation and revision of the manuscript. XC contributed to the design of the study. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by the General Scientific Projects of the Zhejiang Education Department (Y202147905).
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.
Acknowledgments
We would like to acknowledge the TCGA, GEO, GTEx, and HPA database for providing data. We would also like to acknowledge Dr. Jing Luo (Rheumatology Department, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou; Beijing Tsinghua Changgung Hospital, School of Clinical Medicine, Tsinghua University, Beijing, China) for the excellent suggestions in statistics and technical assistance of this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.788670/full#supplementary-material
References
Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-Based Analysis of Lung Single-Cell Sequencing Reveals a Transitional Profibrotic Macrophage. Nat. Immunol. 20, 163–172. doi:10.1038/s41590-018-0276-y
Bedoui, S., Herold, M. J., and Strasser, A. (2020). Emerging Connectivity of Programmed Cell Death Pathways and its Physiological Implications. Nat. Rev. Mol. Cel. Biol. 21, 678–695. doi:10.1038/s41580-020-0270-8
Bergin, A. R. T., and Loi, S. (2019). Triple-negative Breast Cancer: Recent Treatment Advances. F1000Res 8, F1000 Faculty Rev–1342. doi:10.12688/f1000research.18888.1
Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: a Cytoscape Plug-In to Decipher Functionally Grouped Gene Ontology and Pathway Annotation Networks. Bioinformatics 25, 1091–1093. doi:10.1093/bioinformatics/btp101
Bockmeyer, C. L., Christgen, M., Müller, M., Fischer, S., Ahrens, P., Länger, F., et al. (2011). MicroRNA Profiles of Healthy Basal and Luminal Mammary Epithelial Cells Are Distinct and Reflected in Different Breast Cancer Subtypes. Breast Cancer Res. Treat. 130, 735–745. doi:10.1007/s10549-010-1303-3
Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating Single-Cell Transcriptomic Data across Different Conditions, Technologies, and Species. Nat. Biotechnol. 36, 411–420. doi:10.1038/nbt.4096
Capdevila-Busquets, E., Badiola, N., Arroyo, R., Alcalde, V., Soler-Lopez, M., and Aloy, P. (2015). Breast Cancer Genes PSMC3IP and EPSTI1 Play a Role in Apoptosis Regulation. PLoS One 10, e0115352. doi:10.1371/journal.pone.0115352
Carithers, L. J., and Moore, H. M. (2015). The Genotype-Tissue Expression (GTEx) Project. Biopreserv. Biobank. 13, 307–308. doi:10.1089/bio.2015.29031.hmm
Chen, L. C., Wang, L. J., Tsang, N. M., Ojcius, D. M., Chen, C. C., Ouyang, C. N., et al. (2012). Tumour Inflammasome‐Derived IL‐1β Recruits Neutrophils and Improves Local Recurrence‐Free Survival in EBV‐induced Nasopharyngeal Carcinoma. EMBO Mol. Med. 4, 1276–1293. doi:10.1002/emmm.201201569
Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12
Giltnane, J. M., and Balko, J. M. (2014). Rationale for Targeting the Ras/MAPK Pathway in Triple-Negative Breast Cancer. Discov. Med. 17, 275–283.
Giuli, M. V., Giuliani, E., Screpanti, I., Bellavia, D., and Checquolo, S. (2019). Notch Signaling Activation as a Hallmark for Triple-Negative Breast Cancer Subtype. J. Oncol. 2019, 8707053. doi:10.1155/2019/8707053
Global Burden of Disease CancerCollaboration Fitzmaurice, C., Allen, C., Barber, R. M., Barregard, L., Bhutta, Z. A., et al. (2017). Global, Regional, and National Cancer Incidence, Mortality, Years of Life Lost, Years Lived with Disability, and Disability-Adjusted Life-Years for 32 Cancer Groups, 1990 to 2015: A Systematic Analysis for the Global Burden of Disease Study. JAMA Oncol. 3, 524–548. doi:10.1001/jamaoncol.2016.5688
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Harano, K., Wang, Y., Lim, B., Seitz, R. S., Morris, S. W., Bailey, D. B., et al. (2018). Rates of Immune Cell Infiltration in Patients with Triple-Negative Breast Cancer by Molecular Subtype. PLoS One 13, e0204513. doi:10.1371/journal.pone.0204513
He, Y., Jiang, Z., Chen, C., and Wang, X. (2018). Classification of Triple-Negative Breast Cancers Based on Immunogenomic Profiling. J. Exp. Clin. Cancer Res. 37, 327. doi:10.1186/s13046-018-1002-1
Hoffman-Censits, J. H., Grivas, P., Van Der Heijden, M. S., Dreicer, R., Loriot, Y., Retz, M., et al. (2016). IMvigor 210, a Phase II Trial of Atezolizumab (MPDL3280A) in Platinum-Treated Locally Advanced or Metastatic Urothelial Carcinoma (mUC). Am. Soc. Clin. Oncol. 34, 355. doi:10.1200/jco.2016.34.2_suppl.355
Karki, R., and Kanneganti, T.-D. (2019). Diverging Inflammasome Signals in Tumorigenesis and Potential Targeting. Nat. Rev. Cancer 19, 197–214. doi:10.1038/s41568-019-0123-y
Kassambara, A., Kosinski, M., and Biecek, P. (2017). Package ‘Survminer’ [J]. Drawing Survival Curves Using ‘ggplot2’. (R Package Version 03 1).
Kim, S., Jung, W. H., and Koo, J. S. (2012). Differences in Autophagy-Related Activity by Molecular Subtype in Triple-Negative Breast Cancer. Tumor Biol. 33, 1681–1694. doi:10.1007/s13277-012-0424-1
Kim, S., Lee, J., Jeon, M., Lee, J. E., and Nam, S. J. (2016). Zerumbone Suppresses the Motility and Tumorigenecity of Triple Negative Breast Cancer Cells via the Inhibition of TGF-β1 Signaling Pathway. Oncotarget 7, 1544–1558. doi:10.18632/oncotarget.6441
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 28, 882–883. doi:10.1093/bioinformatics/bts034
Lin, N. U., Vanderplas, A., Hughes, M. E., Theriault, R. L., Edge, S. B., Wong, Y.-N., et al. (2012). Clinicopathologic Features, Patterns of Recurrence, and Survival Among Women with Triple-Negative Breast Cancer in the National Comprehensive Cancer Network. Cancer 118, 5463–5472. doi:10.1002/cncr.27581
Liu, L. L., Zhao, H., Ma, T. F., Ge, F., Chen, C. S., and Zhang, Y. P. (2015). Identification of Valid Reference Genes for the Normalization of RT-qPCR Expression Studies in Human Breast Cancer Cell Lines Treated with and without Transient Transfection. PLoS One 10, e0117058. doi:10.1371/journal.pone.0117058
Man, S. M., and Kanneganti, T.-D. (2015). Regulation of Inflammasome Activation. Immunol. Rev. 265, 6–21. doi:10.1111/imr.12296
Nolan, E., Savas, P., Policheni, A. N., Darcy, P. K., Vaillant, F., Mintoff, C. P., et al. (2017). Combined Immune Checkpoint Blockade as a Therapeutic Strategy for BRCA1-Mutated Breast Cancer. Sci. Transl Med. 9, eaal4922. doi:10.1126/scitranslmed.aal4922
Orning, P., Weng, D., Starheim, K., Ratner, D., Best, Z., Lee, B., et al. (2018). Pathogen Blockade of TAK1 Triggers Caspase-8-Dependent Cleavage of Gasdermin D and Cell Death. Science 362, 1064–1069. doi:10.1126/science.aau2818
Park, Y. H., Lee, S. J., Cho, E. Y., La. Choi, Y., Lee, J. E., Nam, S. J., et al. (2019). Clinical Relevance of TNM Staging System According to Breast Cancer Subtypes. Ann. Oncol. 30, 2011. doi:10.1093/annonc/mdz223
Pizato, N., Luzete, B. C., Kiffer, L. F. M. V., Corrêa, L. H., de Oliveira Santos, I., Assumpção, J. A. F., et al. (2018). Omega-3 Docosahexaenoic Acid Induces Pyroptosis Cell Death in Triple-Negative Breast Cancer Cells. Sci. Rep. 8, 1952. doi:10.1038/s41598-018-20422-0
Pohl, S.-G., Brook, N., Agostino, M., Arfuso, F., Kumar, A. P., and Dharmarajan, A. (2017). Wnt Signaling in Triple-Negative Breast Cancer. Oncogenesis 6, e310. doi:10.1038/oncsis.2017.14
Pontén, F., Jirström, K., and Uhlen, M. (2008). The Human Protein Atlas-A Tool for Pathology. J. Pathol. 216, 387–393. doi:10.1002/path.2440
Provance, O. K., Geanes, E. S., Lui, A. J., Roy, A., Holloran, S. M., Gunewardena, S., et al. (2021). Disrupting Interferon-Alpha and NF-kappaB Crosstalk Suppresses IFITM1 Expression Attenuating Triple-Negative Breast Cancer Progression. Cancer Lett. 514, 12–29. doi:10.1016/j.canlet.2021.05.006
Ranstam, J., and Cook, J. A. (2018). LASSO Regression [J]. Br. J. Surg. 105 (10), 1348. doi:10.1002/bjs.10895
Reisenbichler, E. S., Han, G., Bellizzi, A., Bossuyt, V., Brock, J., Cole, K., et al. (2020). Prospective Multi-Institutional Evaluation of Pathologist Assessment of PD-L1 Assays for Patient Selection in Triple Negative Breast Cancer. Mod. Pathol. 33, 1746–1752. doi:10.1038/s41379-020-0544-x
Schmid, P., Adams, S., Rugo, H. S., Schneeweiss, A., Barrios, C. H., Iwata, H., et al. (2018). Atezolizumab and Nab-Paclitaxel in Advanced Triple-Negative Breast Cancer. N. Engl. J. Med. 379, 2108–2121. doi:10.1056/NEJMoa1809615
Shao, W., Yang, Z., Fu, Y., Zheng, L., Liu, F., Chai, L., et al. (2021). The Pyroptosis-Related Signature Predicts Prognosis and Indicates Immune Microenvironment Infiltration in Gastric Cancer. Front. Cel Dev. Biol. 9, 676485. doi:10.3389/fcell.2021.676485
Shi, J., Zhao, Y., Wang, K., Shi, X., Wang, Y., Huang, H., et al. (2015). Cleavage of GSDMD by Inflammatory Caspases Determines Pyroptotic Cell Death. Nature 526, 660–665. doi:10.1038/nature15514
Smyth, G. K. (2005). Limma: Linear Models for Microarray Data, Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Springer, 397–420.
Sullivan, L. M., Massaro, J. M., and D'Agostino, R. B. (2004). Presentation of Multivariate Data for Clinical Use: The Framingham Study Risk Score Functions. Statist. Med. 23, 1631–1660. doi:10.1002/sim.1742
Trapnell, C., Cacchiarelli, D., Grimsby, J., Pokharel, P., Li, S., Morse, M., et al. (2014). The Dynamics and Regulators of Cell Fate Decisions Are Revealed by Pseudotemporal Ordering of Single Cells. Nat. Biotechnol. 32, 381–386. doi:10.1038/nbt.2859
Vikas, P., Borcherding, N., and Zhang, W. (2018). The Clinical Promise of Immunotherapy in Triple-Negative Breast Cancer. Cmar 10, 6823–6833. doi:10.2147/CMAR.S185176
Wang, Z.-Q., Milne, K., Derocher, H., Webb, J. R., Nelson, B. H., and Watson, P. H. (2017). PD-L1 and Intratumoral Immune Response in Breast Cancer. Oncotarget 8, 51641–51651. doi:10.18632/oncotarget.18305
Whitehead, M. J., McCanney, G. A., Willison, H. J., and Barnett, S. C. (2019). MyelinJ: an ImageJ Macro for High Throughput Analysis of Myelinating Cultures. Bioinformatics 35, 4528–4530. doi:10.1093/bioinformatics/btz403
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26, 1572–1573. doi:10.1093/bioinformatics/btq170
Wu, J., Cai, Y., Zhao, G., and Li, M. (2021). A Ten N6-Methyladenosine-Related Long Non-Coding RNAs Signature Predicts Prognosis of Triple-Negative Breast Cancer. J. Clin. Lab. Anal. 35, e23779. doi:10.1002/jcla.23779
Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111
Yang, F., Cao, L., Sun, Z., Jin, J., Fang, H., Zhang, W., et al. (2016). Evaluation of Breast Cancer Stem Cells and Intratumor Stemness Heterogeneity in Triple-Negative Breast Cancer as Prognostic Factors. Int. J. Biol. Sci. 12, 1568–1577. doi:10.7150/ijbs.16874
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zaki, M. H., Vogel, P., Body-Malapel, M., Lamkanfi, M., and Kanneganti, T.-D. (2010). IL-18 Production Downstream of the Nlrp3 Inflammasome Confers protection against Colorectal Tumor Formation. J.I. 185, 4912–4920. doi:10.4049/jimmunol.1002046
Zanker, D. J., Spurling, A. J., Brockwell, N. K., Owen, K. L., Zakhour, J. M., Robinson, T., et al. (2020). Intratumoral Administration of the Toll-like Receptor 7/8 Agonist 3M-052 Enhances Interferon-Driven Tumor Immunogenicity and Suppresses Metastatic Spread in Preclinical Triple-Negative Breast Cancer. Clin. Transl Immunol. 9, e1177. doi:10.1002/cti2.1177
Zhang, Z., Zhang, Y., Xia, S., Kong, Q., Li, S., Liu, X., et al. (2020). Gasdermin E Suppresses Tumour Growth by Activating Anti-tumour Immunity. Nature 579, 415–420. doi:10.1038/s41586-020-2071-9
Zhou, Z., He, H., Wang, K., Shi, X., Wang, Y., Su, Y., et al. (2020). Granzyme A from Cytotoxic Lymphocytes Cleaves GSDMB to Trigger Pyroptosis in Target Cells. Science 368, eaaz7548. doi:10.1126/science.aaz7548
Glossary
AIM2 absent in melanoma 2
BC breast cancer
CASP1/3/4/5/8 caspase1/3/4/5/8
CDF cumulative distribution function
CR complete remission
DEGs differentially expressed genes
DEmiRNAs differentially expressed miRNAs
EPSTI1 epithelial-stromal interaction 1
ER estrogen receptor
FDA Food and Drug Administration
GDSC Genomics of Drug Sensitivity in Cancer
GEO Gene Expression Omnibus
GSDMB gasdermin B
GSDMD gasdermin D
GSDME gasdermin E
GSVA gene set variation analysis
GTEx Genotype-Tissue Expression
GZMA granzyme A
GZMB granzyme B
HER2 human epidermal growth factor receptor 2
HPA Human Protein Atlas
IC50 half-maximal inhibitory concentration
ICI immune cell infiltration
IFITM1 interferon inducible transmembrane 1
IHC immunohistochemistry
KEGG Kyoto Encyclopedia of Genes and Genomes
mUC metastatic urothelial carcinoma
NLR nucleotide-binding domain and leucine-rich repeat receptor
OS overall survival
PCA principal component analysis
PCD programmed cell death
PD programmed cell death
PD-L1 PD-ligand 1
PR progesterone receptor
scRNA-seq single-cell RNA sequencing
TCGA The Cancer Genome Atlas
TME tumor microenvironment
TNBC triple-negative breast cancer
Tregs regulatory T cells
tSNE t-distributed statistical neighbor embedding.
Keywords: pyroptosis, triple-negative breast cancer, molecular subtype, tumor microenvironment, prognosis, single-cell RNA sequencing
Citation: Li C, Pan J, Jiang Y, Wu Y, Jin Z and Chen X (2022) Characterization of Pyroptosis-Related Subtypes via RNA-Seq and ScRNA-Seq to Predict Chemo-Immunotherapy Response in Triple-Negative Breast Cancer. Front. Genet. 13:788670. doi: 10.3389/fgene.2022.788670
Received: 03 October 2021; Accepted: 15 February 2022;
Published: 21 March 2022.
Edited by:
Luis Zapata, Institute of Cancer Research (ICR), United KingdomReviewed by:
Feng Gao, The Sixth Affiliated Hospital of Sun Yat-sen University, ChinaYuchen Yang, Sun Yat-sen University, China
Copyright © 2022 Li, Pan, Jiang, Wu, Jin and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhenlin Jin, amluemhlbmxpbjE2M0AxNjMuY29t; Xupeng Chen, MTM1MDY1NTEyMDVAcXEuY29t
†These authors have contributed equally to this work