A novel pyroptosis-related signature predicts prognosis and response to treatment in breast carcinoma

Schematic diagram of the study design

Fifty-two candidate pyroptosis-related genes with relevance scores greater than 1.0 were identified from the GeneCards database, and CASP8, GZMB and ZBP1 were identified from the literature, for a total of 55 candidate genes to be finalized (Supplementary Table 1). The theory of this relevance score is explained at https://www.elastic.co/guide/en/elasticsearch/guide/current/scoring-theory.html. To identify potential candidates and create a robust signature encompassing 15 pyroptosis-related genes to predict survival, the least absolute shrinkage and selection operator (LASSO) method was utilized (Figure 1A). Subsequently, we validated the model from The Cancer Genome Atlas (TCGA) dataset as the training group in two GSE datasets as validation groups. To further confirm the model’s predictive capability, a meta-analysis was performed on the three datasets (Figure 1B). In addition to assessing the prognosis of patients, we hope that our predictive model will be useful as a clinical guideline. We therefore analyzed the relationship between clinicopathological features and the pyroptosis risk score (PRS) and explored the relationship between PRS and immunotherapy. For patients with different risks, according to our model, we also predicted several different chemotherapeutic agents that could serve as clinical treatment options. Finally, we established a nomogram to evaluate and quantify the prognosis of BC patients based on PRS and other clinical characteristic variables (Figure 1C).

Figure 1. Schematic diagram of the study design. (A) A new prognostic model related to pyroptosis has been established from the literature and databases. (B) The stability and generalisability of the model was verified in an external independent database. (C) Prediction and application of PRS for clinical response to treatment. LASSO, least absolute shrinkage and selection operator; PRS, Pyroptosis-related Risk Scores; AUC, area under the curve.

Development of a pyroptosis-related gene signature for prognosis prediction

To identify the genes among these 55 candidate genes that are most relevant to the prognosis of BC patients, the LASSO Cox regression model was used to identify genes with the greatest prognostic value. To overcome overfitting, tenfold cross-validation was used, with 0.01 selected as the best λ value (Figure 2A). Finally, 15 predicted genes (NLRC4, IRF3, ANO6, GSDMC, TP53, FGF21, IL36B, DHX9, FOXO3, IL36G, IL18, GJA1, MST1, GZMB and GBP1) were identified to have nonzero LASSO coefficients and were incorporated into the gene signature model (Figure 2B). Figure 2C depicts the LASSO coefficient distribution for the genes in the signature. To evaluate the prediction effectiveness of the 15 pyroptosis-related gene-based signatures in the TCGA cohort, a risk score was developed (Supplementary Table 2). Patients with a PRS greater than 0.168 were classified as high-risk, whereas those with a PRS less than 0.168 were classified as low-risk. Figure 3A, 3C show the distributions of risk scores, survival time, and survival status. As shown in Figure 3B, patients with high and low PRSs have very different prognoses. Patients with high-risk scores had a significantly shorter survival time (P < 0.001; hazard ratio (HR) = 6.1; confidence interval (CI) = 4-9.4). The expression of the 15 genes in the gene signature was examined further and is depicted in a heatmap (Figure 3D). The heatmap showed that there were differences in the expression of some of the predicted genes between the two groups. Then, the area under the curve (AUC) was calculated using time-dependent receiver operating characteristic (ROC) curve analysis. The AUCs of the signature for predicting overall survival (OS) at 3, 5 and 7 years reached 0.7211, 0.6748 and 0.7208, respectively. We selected four of the 15 characterized genes with the highest AUCs for ROC analysis and calculated the AUCs of the four genes and PRS (Figure 3E). This also suggested that, when compared to a single gene, the PRS was the most reliable predictor of OS. We examined the relationships of the PRS with clinicopathological characteristics and clinical data to investigate the link between the PRS and clinical data. PRS was discovered to be substantially related to age, tumor stage, histologic subtype, and histologic grade (Table 1).

Figure 2. LASSO regression identifies prognosis-related genes. (A) The adjustment parameter (lambda) in the LASSO model was selected for 10-fold cross-validation by the minimum criterion. Partial likelihood deviation curves were plotted against lambda. Dotted vertical lines were drawn at the optimal values by using the minimum criterion and 1 standard error of the minimum criterion (1-SE criterion). (B) LASSO coefficient profiles of the 55 pyroptosis-related genes. A coefficient profile plot was produced against the log (lambda) sequence. A vertical line was drawn at the value selected using 10-fold cross-validation, where the optimal lambda resulted in 15 nonzero coefficients. (C) Distribution of the LASSO coefficients of the 15 immune-related gene signatures. The horizontal coordinate indicates LASSO coefficients, genes with negative coefficients in this regression indicate prognostic protective genes (green marker), and positive numbers indicate poor prognostic genes (red marker).

Figure 3. Prognostic analysis of the 15 pyroptosis-related gene signature model in TCGA cohort. (A) The distribution and cutoff value of the risk score in the TCGA cohort. (B) Kaplan–Meier curves for the OS of patients in the two groups in the TCGA cohort. (C) The distributions of OS status, OS and the risk score in the TCGA cohort. (D) Heatmap showing differences in the expression of 15 pyroptosis-related genes between high-PRS group and low-PRS group. (E) tROC analysis showed that the PRS was an accurate variable for survival prediction. The four genes shown in the figure have the four highest AUCs among the 15 signature genes. TCGA, The Cancer Genome Atlas; PRS, pyroptosis-related risk score; High-PRS, high pyroptosis-related risk score; Low-PRS, low pyroptosis-related risk score; OS, overall survival; tROC, time-dependent receiver operating characteristic; AUC, area under the curve.

In summary, low-risk patients are more prone to pyroptosis, or pyroptosis activity is greater in the low-risk group; that is, low-risk patients survive longer because the cancer tissue in their bodies displays higher levels of pyroptosis, resulting in more cancer cells dying. However, it is worth noting that such results do not mean that the differential expression of several genes in BC induces or suppresses pyroptosis because some of the 15 genes in this signature are not differentially expressed genes.

Validation of the 15 pyroptosis-related gene signatures in the Gene Expression Omnibus (GEO) cohort

The pyroptosis-related gene signature was subsequently validated in two separate external validation sets to establish its stability and generalizability in various populations. We selected two GSE datasets from the GEO database as validation sets. Patients in the GSE21653 and GSE20685 cohorts were also assigned to one of two groups: high-PRS or low-PRS. The assign formula is the same as in the TCGA cohort. The cutoff value used in both validation cohorts was 0.168, which was consistent with that used in the training cohort. Similar to the training group, in both validation cohorts, we also analyzed the distributions of the risk score (Figure 4A, 4B) and survival status (Figure 4C, 4D). The results were highly consistent with the TCGA group. The gene expression heatmaps are shown in Figure 4E, 4F. Both external validation sets’ analytical results were very consistent with those of the training set. In the two validation cohorts, Kaplan–Meier analysis revealed that low PRS predicted better OS than high PRS. (P =0.0011, Figure 4G; P =0.028, Figure 4H). Next, we conducted a meta-analysis of genetic variables in the training cohort as well as two validation cohorts separated into two groups. The meta-analysis found that all patients with greater PRS had a poorer prognosis than those with lower PRS, as illustrated in Figure 5A. (fixed effect model overall HR = 4.14, 95% CI 3.07–5.58; random effect model HR = 3.84, 95% CI 2.24–6.57). The meta-analysis showed medium heterogeneity, and we attribute this to the fact that the TCGA cohort contains a much larger number of cases than the two validation cohorts. The results of univariate and multivariate Cox regression analyses in the TCGA database are depicted using forest plots (Figure 5B, 5C), and the full data are given in Supplementary Table 3. In the TCGA cohort, the risk score was discovered to be an independent prognostic predictor. (HR=2.9, CI 2.1–4.1, P < 0.01; HR=2.33, CI 1.52–3.56, P < 0.01). Furthermore, in both cohorts, tumor stage and age were independent prognostic indicators (P < 0.01). These results demonstrate the high stability and generalizability of the 15 pyroptosis-associated gene prediction model in BC patients.

Figure 4. Validation of the gene signature in the external validation sets. The distribution and cutoff value of the risk score in the groups in the GSE20685 (A) and GSE21653 (B). The distributions of OS status, OS and the risk score in the GSE20685 (C) and GSE21653 (D). The expression heatmap of the 15 pyroptosis-related prognostic genes between the high and low-PRS groups in the GSE20685 (E) and GSE21653 (F). Kaplan-Meier curves for the OS of patients in the high-risk group and low-risk group in the GSE20685 (G) and GSE21653 (H) cohorts. TCGA, The Cancer Genome Atlas; GSE, gene expression omnibus series; OS, overall survival; PRS, pyroptosis-related risk score; High-PRS, high pyroptosis-related risk score; Low-PRS, low pyroptosis-related risk score.

Figure 5. Verification of the generalizability and stability of the prediction model in TCGA cohort. (A) Meta-analysis of the TCGA training set and 2 external validation sets. (B) Univariate Cox regression analyses of OS in the TCGA training cohort. (C) Multivariate Cox regression analyses of OS in the TCGA training cohort.

Functional and heterogeneity analyses of the 15 pyroptosis-related genes in the high- and low-PRS groups

To determine exactly what role these genes play in terms of their expression and function, further analyses are needed. We analyzed the expression and functions of these 15 pyroptosis-related genes. Some pyroptosis-related genes, including ANO6, TP53, GSDMC, FGF21 and IL36B, had high expression levels in BC patients (Figure 6A). We also found some genes that were not significantly different in expression, including DHX9, NLRC4, IRF3 and GJA1 (Figure 6A). We believe that these genes, which are not significantly differentially expressed but have an impact on prognosis, play a part in the regulation of certain functions. GO and KEGG pathway analyses were also utilized to investigate the probable functions of these genes in the two groups. Interestingly, the TCGA cohort’s pyroptosis-related genes were shown to be enriched in various cancer-related molecular pathways, including immune-related pathways such as T cell activation, T cell regulation, and primary immunodeficiency (Figure 6B, 6C). Reversed-phase protein arrays (RPPAs) were utilized to assess the key pathways in the two groups to further explore tumor heterogeneity between the two groups of patients and to examine the variations in tumorigenesis processes between the two groups. The analysis results showed that tumor purity, proliferation scores, and apoptosis scores differed between the high- and low-risk groups (Figure 6D–6F). All three scores may indicate that high-risk patients have more malignant tumor cells that are prone to metastasis and have more difficulty undergoing apoptosis. The pathway scores, which are protein expression signatures of pathway activity, associated with tumor lineage were from an RPPA in prior work in the literature [18]. Given the highly heterogeneous character of breast cancer, the prognosis of different subtypes of BC patients differs considerably. As a result, we investigated the prognosis of patients in various risk categories among different subtypes of BC patients. The results showed that in patients with HER2 and LUMINAL subtypes, the low-risk group had a better prognosis than the high-risk group (Figure 6H, 6I), while the prognosis of patients with triple-negative breast cancer (TNBC) was not (Figure 6G).

Figure 6. Functional and heterogeneity analysis of 15 pyroptosis-related genes between high- and low-PRS groups. (A) Volcano plot of differentially expressed genes between the low- and high-risk groups. Orange indicates the 15 pyroptosis-related gene signature. (B) The most significant or shared GO enrichment terms in the TCGA cohort. (C) The most significant or shared KEGG pathways in the TCGA cohort. (D) The Boxplots show differences in tumor purity between the high and low-PRS group. (E) The Boxplots show differences in proliferation score between the high and low-PRS group. (F) The Boxplots show differences in apoptosis score between the high and low-PRS group. (G) Results of survival analysis of high- and low- risk groups in TNBC breast cancer patients. (H) Results of survival analysis of high- and low- risk groups in LUMINAL breast cancer patients. (I) Results of survival analysis of high- and low- risk groups in HER2 breast cancer patients. The Kruskal–Wallis test was performed to calculate the P-value. GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PRS, pyroptosis-related risk score; High-PRS, high pyroptosis-related risk score; Low-PRS, low pyroptosis-related risk score.

Response to chemotherapy and immunotherapy of high- and low-PRS patients

Systemic chemotherapy remains the foundation of BC treatment. However, some refractory, advanced tumors lack more chemotherapy options. Therefore, we used two chemotherapy drug response databases (PRISM and CTRP) to identify new highly sensitive chemotherapy candidates for high PRS patients. First, the AUC values and PRSs were analyzed by Spearman correlation analysis to identify chemicals with negative correlation coefficients (Spearman r < -0.01 for CTRP and -0.01 for PRISM). Next, an examination of drug response differences between the high- and low-PRS groups was performed to find drugs with lower estimated AUC values in the high-PRS group [log2-fold change (FC) >0.07]. It is essential to note that a lower AUC indicates higher drug sensitivity. Three CTRP-derived compounds (cucurbitacin-I, PI-103 and SB743921) and three PRISM-derived compounds (including temsirolimus, mitoxantrone and ispinesib) were found to be potentially sensitive in the high-PRS group. In the high-PRS group, all of these drugs displayed a negative connection with the PRS and lower calculated AUC values. (Figure 7A for CTRP and Figure 7B for PRISM). For the chemotherapy drug analysis, we also analyzed the responses to some of the known and widely used chemotherapy drugs, cisplatin and fluorouracil. The results showed that these two common chemotherapy drugs also differed significantly in high- and low-risk populations. The low-risk populations had lower AUC values, indicating greater sensitivity to the two common chemotherapeutic agents. (Supplementary Figure 1). Multiple-perspective studies were undertaken to assess the therapeutic potential of these compounds in BC to further evaluate whether the compounds we found had potential clinical utility. First, a detailed assessment of the literature in PubMed was conducted to discover experimental and clinical evidence of the six prospective medications in the treatment of BC. Second, the log2 FC values of the differential mRNA and protein expression levels of genes relevant to drug targets between tumor and normal tissues were determined, and a larger log2 FC value suggested better potential for BC therapy. Third, CMap analysis was utilized to confirm compounds whose gene expression patterns were oppositional to the BC-specific expression patterns (the expression of some genes was raised in tumor tissue but lowered after treatment with particular drugs). Lower CMap scores imply that these chemicals might have therapeutic potential for BC (Figure 7C). The approaches described above are based on published literature [19].

Figure 7. Response to chemotherapy and immunotherapy of high- and low-IRS patients. (A, B) The bubble plot shows the degree of negative correlation between AUC values and PRS scores of three CTRP-derived compounds (A) and three PRISM-derived compounds (B). Longer lines in the graph indicate a stronger negative correlation, predicting greater drug sensitivity. Boxplots indicate the results of differential drug response analysis between the PRS high and PRS low groups for the three CTRP-derived compounds (A) and the three PRISM-derived compounds (B). Note that lower values on the y-axis of boxplots imply greater drug sensitivity. (C) Evidence from multiple sources to identify the most promising therapeutic agents for the high-IRS group. Six compounds come from CTRP and PRISM are shown on the diagram, respectively. (D) Proportion of immune-responsive and nonimmune-responsive populations in the PRS high- and low-expression groups. (E, F) Survival analysis graph showing differences in survival between PRS-high and IRS-low groups in nonimmune-responsive (E) and immune-responsive (F) patients. (G) Differences in immune infiltration scores of 24 immune cell types in high and low PRS groups in TCGA database. * means p <0.05, ** means p <0.01, *** means p<0.005, **** means p<0.001.

In addition to chemotherapy, immunotherapy is a new direction for BC treatment. Additionally, in our previous GO and KEGG analyses, we found that these genes appear to be inextricably linked to the immune system. Therefore, we sought to determine whether there was a difference in the response to immunotherapy between the high- and low-risk groups. Figure 7D depicts the predicted response to immunotherapy of patients with high and low PRS. In the high-PRS group, only 63% of patients were expected to react to immunotherapy, but in the low-PRS group, this proportion was 80%. Based on these findings, we hypothesize that patients with low PRS may be more susceptible to immunotherapy. Patients with a high PRS had a poorer prognosis in both groups when survival analysis was performed individually in the immunotherapy-responsive and immunotherapy-nonresponsive groups (Figure 7E, 7F). These findings imply that patients with high PRS are not especially sensitive in terms of immunological response. This may be because of the different expression of genes associated with immunity. To further explore the relationship between the PRS and immune status, the immune infiltration scores in TCGA were calculated by using the Immune Cell Abundance Identifier (ImmuCellAI), which is used to precisely assess the abundance of 24 immune cell types, including 18 T cell subsets. (Figure 7G). It is encouraging to note that most of the immune cells, especially the various types of T cells, demonstrated significant differences in the different PRS groups (adjusted P < 0.05). These findings imply that the 15 pyroptosis-related genes we selected are potentially related to immune infiltration.

Developing a prediction nomogram for BC patients

Independent risk indicators were used to develop a risk assessment nomogram to provide a clinically relevant strategy for predicting the OS probability of BC patients (Figure 8A). These indicators included age, tumor stage, PR status, HER2 status, ER status and the risk score related to pyroptosis. Figure 8B–8D show calibration plots for 3-, 5-, and 7-year survival probabilities in the TCGA cohort, indicating that the nomogram had a high degree of accuracy. This indicates that our nomogram has a high predictive value.

Figure 8. Predictive nomogram for PRS and clinical features. (A) The 15 pyroptosis-related prognostic model for predicting 3-, 5-, and 7-year OS in BC patients. The independent risk factors were used to build a risk estimation nomogram to predict the probability of OS in BC patients. (B) The calibration plots for 3-year survival probabilities in the TCGA cohort. (C) The calibration plots for 5-year survival probabilities in the TCGA cohort. (D) The calibration plots for 7-year survival probabilities in the TCGA cohort.

Data collection

The training set consisted of the RNA data collected and clinical information from 1007 BC patients acquired from The Cancer Genome Atlas (TCGA). In addition, the GSE21653 (n=266) and GSE20685 (n=327) datasets, all of which contain gene expression files from the Gene Expression Omnibus (GEO) database, were utilized as external validation sets. Only patients with overall survival (OS) times of more than 30 days were included in all three datasets. Candidate pyroptosis-related genes were identified mostly from the literature and the GeneCards database (https://www.genecards.org/) using the term pyroptosis, and genes with relevance scores greater than 1.0 were selected. In the list of candidate genes, CASP8 [57], GZMB [58] and ZBP1 [59] were screened by reading extensive literature and were not in the Genecards database screening results. All candidate pyroptosis-related genes are provided in Supplementary Table 1. The infiltration scores of 24 immune cells were obtained from ImmuCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) [60] and normalized in R software (version 4.0.3). ImmuCellAI is a website that estimates the abundance of 24 immune cells using gene expression datasets such as RNA-Seq and microarray data. It may be used to quantify the difference in immune cell infiltration across groups and predict patients’ responses to immune checkpoint blockade treatment. It works by obtaining a reference expression profile for each cell type from the GEO database and curating a specific gene set from a publication as a gene signature. Then, the total expression deviation of the gene signatures in the input expression profiles from the reference expression profiles of the 24 immune cell types was calculated. Deviations were then assigned to the corresponding immune cell type based on the enrichment fraction of its genetic profile, which was calculated using the single sample gene set enrichment analysis (ssGSEA) algorithm [60]. Drug sensitivity data of cancer cell lines (CCLs) were gathered from the Cancer Therapeutics Response Portal (CTRP v.2.0, released October 2015) and the PRISM Repurposing dataset (19Q4, released December 2019).

Study design

The training set was the TCGA cohort, while the validation sets were the GSE20685 and GSE21653 cohorts. Fifty-five potential pyroptosis-related genes were found mostly through research articles and search results of the GeneCards database, with pyroptosis as the keyword, and genes with correlation scores greater than 1.0 were selected. The 55 genes were then evaluated using the R package glmnet (version: 4.0-2) for least absolute shrinkage and selection operator (LASSO) Cox regression to filter down the possible pyroptosis-related genes. The 55 pyroptosis-related genes were found to have nonzero coefficients in the model, and the samples were divided into high- and low-risk groups based on the cutoff value of 0.1677 obtained using the survminer R package’s surv cutpoint function (Version: 0.4.3). The formula of the pyroptosis risk score (PRS) was as follows:

PRS = sum of coefficients × normalized expression level of pyroptosis-related genes.

Finally, independent risk variables found by multivariate Cox regression analysis were used to create a nomogram for predicting patient OS. Finally, to evaluate the performance of the nomogram, calibration plots were created. The concordance index (C-index) was used to examine the consistency between the model’s projected probability and the actual outcomes. The R package rms was used to plot the nomogram and calibration plots (Version: 4.0.2).

Gene set enrichment analysis (GSEA) and functional enrichment analysis of pyroptosis-related genes

GSEA was used to discover differences in the mechanisms in BC patients to analyze the potential pathways of pyroptosis-related genes. Based on the risk score, genes with P < 0.05 and |log2-fold change (FC)| ≥1 were considered to be substantially differentially expressed between the high- and low-risk groups. The clusterProfiler R tool was used to perform GO and KEGG analysis on these pyroptosis-related genes [61].

Estimation of chemotherapy and immunotherapy responses

Data on immunotherapy responses (anti-PD1 or anti-CTLA4) were obtained from ImmuCellAI, a gene set signature-based technique developed for predicting immunotherapy responses with high accuracy by evaluating gene expression data [60]. Through this tool, we divided BC patients into immunotherapy-responsive and nonresponsive groups. The chi-square test was then used to analyze the immunotherapy response group and the PRS group. The infiltration scores of the 24 immune cells from ImmuCellAI were also used to compare the differences between the high- and low-PRS groups.

Proteomic analysis is one of the commonly used tools in tumor pathophysiology. The TCGA project team has made extensive use of reverse-phase protein array (RPPA) technology to perform tumor proteomic studies. This microarray allows analysis of the expression and trends of multiple marker proteins in the sample. By integrating RPPA microarray data from TCGA and several independent tumor research projects, The Cancer Proteome Atlas (TCPA) is provided to facilitate researchers to view and analyze the visualized microarray data. The article was published in nature methods [62]. The data we used in this study were the breast cancer RPPA microarray data compiled by Ciriello G et al. and can be downloaded directly from the supplemental data in the article [18].

The connectivity map (CMap) is a resource that uses cellular responses to perturbation to find relationships between diseases, genes, and therapeutics. The list of up- and downregulated differentially expressed genes obtained from the experimental analysis was compared with the database reference dataset using CMap; a correlation score (-100 to 100) was calculated based on the enrichment of differentially expressed genes in the reference gene expression profile; a positive number indicated that the up- and downregulated differentially expressed genes were similar to the reference gene expression profile; a negative number indicated that the up- and downregulated differentially expressed genes may be opposite to the reference gene expression profile; finally, the reference gene expression profile was ranked according to the correlation score [63].

CTRP, which contains sensitivity data for 481 chemicals in 835 CCLs, and the PRISM Repurposing dataset (19Q4, published December 2019), which contains sensitivity data for 1448 chemicals in 482 CCLs, were used to gather drug sensitivity data for CCLs. The AUC value is presented in these two datasets as a measure of drug sensitivity, with lower AUC values suggesting more sensitivity to treatment. Each TCGA sample’s drug sensitivity was assessed using the R package pRRophetic (Version: 4.15-1), which has a built-in ridge regression model that was used to predict the chemotherapy response of clinical samples based on their expression profiles [19, 64]. The AUC values and PRSs were then examined using Spearman correlation to identify compounds with negative correlation coefficients (Spearman r < −0.25 for CTRP and −0.30 for PRISM). Finally, an evaluation of drug response differences between the high-PRS (highest decile) and low-PRS (lowest decile) groups was conducted.

Statistical analysis

R version 4.0.3 was used for all statistical studies (2020-10-10). The Mann–Whitney U test and the Pearson chi-square test were used to compare continuous and categorical variables between the training cohort and validation cohort, respectively. GSEA was performed to test the TCGA cohort findings with the “c5.go.v7.2.entrez.gmt” gene set using the clusterProfiler R package (Version: 3.18.0; https://bioconductor.org/packages/clusterProfiler/) [65]. Using the meta R package (Version: 4.15- 1; https://cran.r-project.org/web/packages/meta/index.html), a meta-analysis was performed to assess the prognostic significance of the gene signature across all datasets. The primary predictive variables of OS were identified using multivariate and univariate Cox regression models (P < 0.05). The nomogram and calibration curve were plotted using the rms R package. The log-rank test was used to compare Kaplan–Meier survival curves. The volcano plot and heatmap were created using the ggplot2 R package. Pearson’s r correlation and Spearman’s rank-order correlation were used to calculate correlations between two continuous variables. The missing AUC values were imputed using K-nearest neighbor (KNN) imputation. P < 0.05 was regarded statistically significant.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.