The expressed gene RNA sequencing (RNA-seq) datasets and corresponding clinical data for BC were extracted from The Cancer Gene Atlas (TCGA) database (https://portal.gdc.cancer.gov/) for breast tumor and normal breast tissue samples. The models were validated by downloading GSE42568 and GSE88770 datasets from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/). 268 ARGs were obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) using the keyword Apoptosis-homo sapiens: hsa04210 pathway (https://www.kegg.jp/entry/hsa04210) and previous related literature [12, 13]. CellMiner was utilized to retrieve RNA-seq expression data and compound activity information for the developmental therapeutics program (DTP) NCI-60 dataset [14] (Fig. 1).

Fig. 1.

Research process. GEO, gene expression omnibus.

ARGs differentially expressed in BC were identified using the “limma” package (version 3.14) in R software (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria). “clusterProfiler” in R Gene Ontology (GO) and KEGG pathway analyses were performed using [13, 15, 16, 17]. The heatmap, ggpubr, and ggplot packages were used to visualize the genes.

Using univariate and multivariate Cox regression, the model was constructed by identifying apoptotic genes associated with BC prognosis. Eight apoptotic genes that were significantly associated with prognosis were finally retained. Risk scores for all patients were calculated as follows: $Riskscore={\sum }_{i=1}^n(Coe{f}_i\times Exp{r}_i)$(n, Expr${}_{\text{i}}$, and Coef${}_{\text{i}}$ represented gene number, level of gene expression, and Coefficient values calculated from multiple Cox regression analysis, respectively). A median risk score was used to categorize patients into high- risk group (HRG) and low-risk group (LRG).

Based on the CIBERSORT algorithm, 22 immune cells were identified in each BC specimen and the sum of the composition of each immune cell was equal to 1 [18]. The levels of infiltration of different immune cells and multiple immune pathway scores were analyzed using single sample gene enrichment analysis (ssGSEA) in the “GSVA” R package (version 3.14) [19]. Subsequently, the immune pathways of the two groups were then examined.

Immunotherapy targets commonly used in clinical immunotherapy were used to assess the sensitivity of the HRG to immunotherapy. A cohort of advanced renal hyaline cell carcinoma immunotherapy treated with PD-1 blockade was included (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7499153/), excluding patients with incomplete responses; finally, 311 patients were included in the analysis.

Age, stage, tumor size (T), lymph node status (N), metastasis (M), and risk scores were evaluated in terms of their prognostic value through univariate and multivariate Cox regression analyses in the ‘survival’ R package. (https://cran.r-project.org/web/packages/survival). Based on Kaplan-Meier (K-M) analysis and log-rank test, overall survival (OS) was compared between the two groups. A time-dependent Receiver Operator Characteristic (ROC) curve analysis was performed using the “timeROC” package (https://cran.r-project.org/web/packages/timeROC).

A clinical nomogram was generated using the “rms” R package (https://cran.R-project.org/package=rms) to predict individual survival. In order to assess the agreement between actual and predicted survival, calibration curves were used. Using external GEO data, the predictive power of the risk model was validated.

Gene expression data and the Z-scores for drug sensitivity were downloaded from Cell Miner (National Cancer Institute) for 60 different tumor cell lines and included 23,808 gene expression profiles and 352 drugs approved by the FDA or in clinical trials [14]. The relationship between the model genes and drug treatment response was investigated. Pearson’s correlation coefficients between genes and drugs were calculated and screened for statistical significance at p 0.01.

Sections of BC tissue and adjacent tissue were embedded in paraffin and fixed in formalin. Sections were dewaxed after baking at 65 °C for 1 h, passed through a gradient alcohol series (Absolute ethanol, 90% ethanol, 80% ethanol, and 70% ethanol), rehydrated and washed. The paraffin-embedded sections were boiled in an antigen repair solution for 10 min for antigen repair at a temperature of 96–98 degrees and then incubated with the following primary antibodies TUBA1C (1:50; D123443; Sangon Biotech, Shanghai, China), BCL2A1 (1:50; D220158; Sangon Biotech), TUBA3D (1:100; D110022; Sangon Biotech), EMP1 (1:100; D162994; Sangon Biotech), F2 (1:100; A00044; Boster Bio, Pleasanton, CA, USA), GSN (1:50; D160423; Sangon Biotech), PMAIP1 (1:100; A9801; ABclonal, Woburn, MA, USA), TP53AIP1 (1:100; D262005; Sangon Biotech) overnight at 4 ℃. The secondary antibodies at 37 ℃ for 30 min. Furthermore, samples were back stained with hematoxylin, dehydrated with an ethanol-xylene series, and covered for 30 minutes with neutral balm (Biosharp; BL704A). Gene staining intensity was rated as follows: no staining (negative), 0; light yellow (weakly positive), 1; brownish-yellow (positive), 2; and tan brown (strongly negative), 3. The final score of each gene was determined by a qualified pathologist blinded to the study to determine its expression as low or high.

Our statistical analyses were performed using R software (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria), the Perl language for processing the data matrix, and GSEA4.1 for analyzing the KEGG enrichment pathway in the two groups in GEO. p 0.05 was considered statistically significant.

The expression of 268 ARGs obtained from the KEGG pathway and previous literature was in 1096 breast tumor samples and 112 non-tumor breast tissue samples in TCGA, using the limma package. 74 differential ARGs were screened based on false discovery rate 0.05, and Log fold change $\ge$1, and 29 down-regulated and 45 up-regulated genes were obtained (Fig. 2A,B).

Fig. 2.

Analysis of apoptosis-related genes (ARGs) and functions differentially expressed in breast cancer (BC) and normal tissues. (A,B) Volcano and box line plots of differentially expressed apoptosis-associated genes (p 0.05). (C,D) Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathway analysis.

In the KEGG pathway enrichment analysis of ARGs, five categories of related pathways were identified (Fig. 2C). Apoptosis pathways, TNF signaling pathways, mitogen-activated protein kinase (MAPK) signaling pathways, NOD-like receptor signaling pathways, and NOD-regulatory pathways were enriched. In addition, ARGs were enriched in breast cancer-related pathways. In the GO analysis (Fig. 2D), endogenous and exogenous apoptosis-related pathways and apoptosis regulatory processes were enriched in biological processes (BP), immune T cell activation and interleukocyte adhesion. Molecular functions (MF) were enriched for protein heterodimeric activity and cytokine receptor binding.

From the 74 different ARGs, 17 were identified as being related to prognosis through univariate Cox analysis. Furthermore, a multivariate analysis was performed to eight ARGs significantly associated with prognosis, and the expression coefficients of 8 independent risk genes were calculated using multivariate Cox regression ratios. The following formula was developed to calculate the prognosis of each patient with BC using the eight genes:

$Riskscore=\left(-0.18513\times expressionlevelofPMAIP1\right)$ $+\left(-1.08828\times expressionlevelofTP53AIP1\right)$ $+\left(-0.12899\times expressionlevelofTUBA3D\right)$ $+\left(0.381992\times expressionlevelofTUBA1C\right)$ $+\left(-0.3809\times expressionlevelofBCL2A1\right)$ $+\left(0.446263\times expressionlevelofEMP1\right)$ $+\left(-0.31589\times expressionlevelofGSN\right)$ $+\left(2.312882\times expressionlevelofF2\right)$

With the increasing number of fatalities, the risk score increased (Fig. 3D,E). K-M survival analysis found significant OS in the LRG compared to the HRG (Fig. 3B). For the area under the curve (AUC) at 3, 5, and 10 years, the survival time-dependent ROC curve predictive model showed excellent results with values of 0.73, 0.736, and 0.783 (Fig. 3A).

Fig. 3.

Construction of risk scores and their relationship to clinical characteristics. (A) The Receiver Operator Characteristic (ROC) curves to assess the prognostic ability of the risk model. (B) Two groups Kaplan-Meier (K-M) survival curves. (C) Forest plot for univariate Cox regression analysis. (D) Survival distribution of BC patients in The Cancer Gene Atlas (TCGA). (E) Risk score distribution of BC patients. (F) Forest plot of multivariate Cox regression analysis. (G) Calibration curve checking the accuracy of the model. (H) Nomogram predicts survival in BC Patients.

Cox analyses were conducted on age, staging, and T, tumor size; N, lymph node status; M, metastasis (TNM) as univariate and multivariate variables. Both clinical variables could be used as predictors in the univariate Cox regression analysis (Fig. 3C); however, the multifactorial Cox regression analysis showed that age (hazard ratio [HR] 1.029 [1.014–1.046] p 0.01) and risk score (HR 1.160 [1.111–1.212] p 0.01) can be used as independent predictors of clinical outcomes (Fig. 3F). Furthermore, a nomogram was created through multivariate Cox analysis to predict the 3, 5, and 10-year overall survival rates based on the eight model genes. Our calibration curves verify the accuracy of the nomograms (Fig. 3G,H).

A new expression matrix was obtained by combining the GSE42568 (104 BC specimens) and GSE88770 (117 BC specimens) datasets to validate model accuracy. Heatmaps of ARGs (Fig. 4A) showed LRG and HRG for patients with BC based on risk scores; as the risk score increased, the mortality rate also increased (Fig. 4B). According to the K-M curves (Fig. 4C), survival rates were significantly higher in the LRG, whereas the modeled AUC values at 3, 5, and 10 years were 0.795, 0.748, and 0.753, respectively (Fig. 4D). In the LRG, immune pathways were enriched (Fig. 4E).

Fig. 4.

GEO dataset validate. (A) Heatmap of 8 risk genes. (B) Survival distribution of BC patients in GEO dataset. (C) K-M survival curves. (D) The time-dependent ROC curves. (E) Analysis of KEGG-related pathways of risk model in GEO dataset. GEO, gene expression omnibus; BC, breast cancer; ROC, receiver operator characteristic; KEGG, kyoto encyclopedia of genes and genomes.

To better understand their immunological differences, we performed CIBERSORT, ssGSEA, and immune-related pathway analysis. According to CIBERSORT analysis, the LRG had significantly higher levels of plasma cells, CD8 T cells, activated CD4 memory T cells, activated natural killer (NK) cells, T-cells regulatory (T regs), and $\gamma$$\delta$T (Fig. 5A). The expression of immune cells in 28 was analyzed via ssGSEA (Fig. 5B). Overall, the levels of different immune cells were typically elevated in the LRG. For example, the LRG enriched its immune function, for instance, by enhancing its NK cell, NK T-cell and activated CD4 T, CD8 T and B-cell. Regarding the 13 immune-related pathway scores, the LRG had higher scores (Fig. 5C).

Fig. 5.

Immunological analysis of two groups. (A) Cibersort analyses the proportion of immune cells. (B) single sample gene set enrichment analysis (ssGESA) analysis of immune cell expression. (C) Scores of 13 Immune-related pathways in both groups. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns, no significant.

To assess the sensitivity of the risk models to immune checkpoints, we used commonly used immune checkpoint targets in the two groups to verify CD274 (PD-L1), CTLA4, PDCD1, and PDCD1LG2 (PD-L2) expression (Fig. 6A–D,G). A PD-1 blocking treatment cohort, employed to verify the immunotherapy response of the model, revealed that longer life in the LRG (Fig. 6E) and a higher immunotherapy response than the HRG (Fig. 6F). These data indicate the excellent response of this model to immunotherapy.

Fig. 6.

Comparison of immune checkpoints between the two groups and validate PD-L1 immunotherapy response. (A–D,G) Comparison of immunological checkpoints. (E) Survival curves of risk models in the PD-L1 treatment cohort. (F) Results of PD-L1 treatment in both groups of patients. *p 0.05; ***p 0.001.

Interestingly (Fig. 7), with an increase gene expression, both GSN and BCL2A1 showed a strong correlation with vemurafenib, dabrafenib, and PD-98059. PMAIP1 showed a strong correlation with palbociclib. These results indicate that these genes are sensitive to these drugs; however, EMP1 and GSN were negatively correlated with palbociclib but also to oxaliplatin, which indicated the emergence of drug resistance with high expression of this gene. EMP1 was resistant to dexrazoxane.

Fig. 7.

Drug sensitivity analysis of model genes.

Immunohistochemistry (IHC) techniques were used to test expression levels of risk model genes in BC tissues and adjacent tissues. Our findings validated the predicted results (Fig. 8) that PMAIP1, TUBA3D, TUBA1C, F2, and BCL2A1 were highly expressed in BC tissues, however, GSN, EMP1, and TP53AIP1 expression showed contrasting results.

Fig. 8.

Image of gene expression in BC tissue and adjacent tissue Immunohistochemical. (A) Representative images of eight ARGs expression in BC tissue and their adjacent normal tissues (right magnification 400$\times$, left partial zoom); BC tissue (right magnification 400$\times$, left partial zoom). Scale bar: 50 µm. (B) Heatmap of the expression of eight ARGs in 10 paired BC tissue and their adjacent tumor tissues. A, adjacent; T, tumor.