Construction of a novel five programmed cell death-related gene signature as a promising prognostic model for triple negative breast cancer

Introduction

Breast cancer is the most common malignancy among women worldwide and a leading cause of cancer-related death (Giaquinto et al., 2022). Triple negative breast cancer (TNBC) accounts for approximately 10–20% of breast cancer cases and lacks expression of the estrogen receptor, progesterone receptor, and human epidermal growth factor receptor-2. Substantial strides in regular screening, surgical resection, adjuvant chemotherapy, precise radiotherapy, and immune regulators have resulted in decreasing the 5-year mortality of TNBC over the last three decades (Korde et al., 2021; Lee, 2023). However, TNBC is a more aggressive subtype of breast cancer that usually progresses rapidly, develops drug resistance, metastasis, and relapses, and remains a challenge for clinicians to treat. These, along with the fact that TNBC represents a heterogeneous group of diverse subtypes with own biological and molecular characteristics (Jiang et al., 2019), highlight the need to identify more sensitive and effective prognostic biomarkers as surrogates of clinical and pathological features (Duffy et al., 2015).

Cell death includes accidental cell death and programmed cell death (PCD). Unlike accidental cell death caused by extreme physical, chemical, or mechanical damage, PCD is a conserved mechanism of cell suicide that is controlled by various pathways, including apoptosis, ferroptosis, cuproptosis, and so on (Wang et al., 2024; Chen et al., 2024). Apoptosis is the most extensively studied form of PCD and is morphologically characterized by cell shrinkage, organelles destruction, cytoplasmic condensation, and nuclear fragmentation (Saraste & Pulkki, 2000). Ferroptosis is characterized by disordered intracellular iron ion flow and significant accumulation of reactive oxygen species and lipid peroxide levels without the need for caspases (Stockwell et al., 2017). Cuproptosis is a new form of PCD triggered by accumulation of intracellular copper level and aggregation of mitochondrial lipoylated proteins (Tsvetkov et al., 2022). With the deepening of PCD research, various PCD modes have been discovered and studied, such as autophagy, necroptosis, pyroptosis, disulfidptosis, oxeiptosis, alkaliptosis, parthanatos, netotic cell death, entotic cell death, and lysosome-dependent cell death (Chen et al., 2024).

Mounting evidence showed that PCD contributed to carcinogenesis and cancer progression, and thus have attracted increasing attention in oncology research. Studies with regard to PCD-related gene signature as a prediction model are emerging. For example, Chen et al. (2024) screened 149 PCD-related differentially expressed genes, of which INHBA, LRRK2, HSP90AA1, HSPB8, and EIF2AK2 were identified as the hub genes of esophageal squamous cell carcinoma. Cao et al. (2023) constructed a 16 PCD-related gene model with potential in predicting prognosis and response to immune checkpoint inhibitors in cancer. Dong et al. (2024) identified seven PCD-related genes to establish the PCD-related risk score for the advanced non-small cell lung cancer model, effectively stratifying overall survival in patients with advanced non-small cell lung cancer. However, the implication of PCD-related genes in TNBC is still unknown. It is therefore necessary to elucidate the molecular significance of PCD-related genes in TNBC and their correlations with survival outcomes and treatment efficacy. In this study, five PCD-related genes associated with prognosis were identified, and a risk model was constructed to assess the prognosis and accurately stratify TNBC patients to enhance survival outcomes. The findings of the present work might offer useful insights in digging sensitive and effective biomarkers for TNBC prognosis prediction and establishing accurate prognostic model in clinical management.

Materials and Methods

Result

Clustering and identification of PCD-related genes

Information of 109 TNBC patients was downloaded from the TCGA database and compared with 1160 PCD-related genes from previous publications (Fig. 1). Univariate Cox regression analysis was used to screen PCD-related genes associated with survival, and a total of 29 genes (VDAC1, ATP6V0D2, DAPK2, ZNF385A, MILR1, SPTLC2, PLA2G15, LAMP3, ATP6V0D1, MT1G, LIPT1, MUL1, LGALS8, HPS6, PINK1, SVIP, CREB3L1, ACKR3, SERPINE1, LRSAM1, GPR137, BCL2A1, BCL2L10, NCK2, CTSD, CDKN1A, FZD9, RRP8, and TPCN1) were enrolled. CDF was applied to categorize the optimal number of clusters. When K was identified as 2, clustering results were relatively stable and CDF delta displayed the slowest decreasing trend, demonstrating that the differences were most significant when TNBC samples were divided into two groups (Figs. 2A–2B). As shown in the present work, samples were clustered into two groups: C1 and C2, as exhibited in the heat map (Fig. 2C). In addition, clustering consistency plots with other Ks (from 3 to 10) samples showed lower average intra-group consistency, when compared with the optimal number of clusters (Fig. 2D). Kaplan-Meier curve was used to analyze the survival difference between two groups (Fig. 2E). Subgroup C1 had remarkable better prognosis with respect to subgroup C2, indicating that survival difference could be related to the differentially expressed PCD-related genes between C1 and C2. Therefore, the volcano map and heat map of the differentially expressed PCD-related genes between C1 and C2 were drawn by using the limma package in R language (Figs. 2F–2G).

Construction of PCD-related prognostic model for TNBC patients

According to the preliminarily obtained gene dataset, a prognostic model was constructed via LASSO analysis by using the glmnet package in R language. After picking the optimal penalty parameter, λ associated with the minimum five fold cross-validation, a total of five PCD-related genes, namely, SEPTIN3, SCARB1, CHML, SYNM, and COL5A3, were then chosen as candidate genes (Figs. 3A, 3B). Survival information of these genes were evaluated by using the survival package in R language. In particular, SEPTIN3, SCARB1, CHML, and SYNM were associated with better prognosis; whereas COL5A3 was responsible for poor prognosis (Fig. 3C). Based on the LASSO analysis and multivariate Cox regression analysis, these five PCD-related genes were finally used for prognostic model construction. Risk score for each patient was calculated as follows: risk score = −0.0799972791070521 * SCARB1 + 0.196866693004799 * COL5A3 − 0.153621818676914 * CHML − 0.141523874769071 * SEPTIN3 − 0.0623684500004447 * SYNM. Then, 109 TNBC patients in the TCGA cohort were stratified into high-risk group and low-risk group, according to the optimal cut-off value (−0.918997797) of risk score calculated by using the maxstat package in R language. When compared with the patients in the low-risk group, patients in the high-risk group showed a significantly poor survival probability (Fig. 3D). By using the pROC package in R language, time-dependent ROC analysis was carried out to further evaluate the prediction efficiency of the constructed gene signature, with the areas under curve (AUC) of 365 days, 1,095 days, and 1,825 days being 0.90, 0.89, and 0.89 (Fig. 3E), respectively. Relationship among risk score, gene expression, and survival status was also checked. Increased expressions of SCARB1, SEPTIN3, CHML, and SYNM were observed in patients from low-risk group, and elevated level of COL5A3 was found in patients from high-risk group (Figs. 3F–3G), indicating that SCARB1, SEPTIN3, CHML, and SYNM were protective factors and COL5A3 was a risk factor.

Validation of PCD-related prognostic model for TNBC patients

To solidify the accuracy of PCD-related prognostic model, the GSE58812, GSE21653, and GSE135565 were used as validation cohorts. Patients from these datasets were therefore divided into high-risk group and low-risk group, based on the optimal cut-off value of calculated risk score as described above. Correlations among risk score, gene expression, and survival status were also evaluated. As shown in the Kaplan-Meier curves, patients in the high-risk group had a significantly worse overall survival, with respect to patients in low-risk group (Fig. 4). Moreover, patients in the low-risk group expressed increased levels of SCARB1, SEPTIN3, CHML, and SYNM; whereas patients in the high-risk group displayed elevated levels of COL5A3 (Fig. 4). The results collectively showed that the constructed PCD-related prognostic model was accurate for survival outcome prediction in TNBC patients.

Elucidation of pathways and molecular functions in the TCGA cohort

To better understand the differences in signaling pathways and molecular functions between the high-risk group and the low-risk group, gene set enrichment analysis was carried out. Volcano map of the differentially expressed genes between the high-risk group and the low-risk group were drawn by using the limma package in R language (Fig. 5A). KEGG analysis indicated that many metabolism-related molecular pathways were enriched, including protein and fat digestion and absorption, glycerolipid, glycerophospholipid, and arachidonic acid metabolism (Fig. 5B). Pathways related to tumor cell invasion and metastasis, such as the extracellular matrix (ECM)-receptor interaction, and relaxin signaling pathway, were also significantly enriched. GO analysis of biological processes, cell components, and molecular functions showed that the differentially expressed genes between the high-risk group and the low-risk group were mainly involved in cell proliferation, extracellular matrix, protein binding, and collage binding (Figs. 5C–5E). Interestingly, several cell components responsible for intercellular communication were also detected, such as extracellular region, exosomes, and extracellular vesicles. GeneMANIA database was applied to predict the genes with similar functions to these five genes. A total of 20 genes were selected, and function analysis indicated that they took part in cholesterol transport, collagen, and metabolism process (Fig. 5F). Gene set enrichment analysis demonstrated that the high-risk group was significantly enriched in coagulation, angiogenesis, P53 pathway, apoptosis, and epithelial-mesenchymal transition (Fig. S1). Collectively, the above results showed that the differences between the high-risk group and the low-risk group classified by the PCD-related prognostic model were mainly associated with tumor metabolism, cancer cell proliferation, invasion and metastasis, and tumor microenvironment.

Immune infiltration, and drug susceptibility analysis

Immune cells play a significant role in oncogenesis, progression, and prognosis of TNBC. CIBERSORT algorithms were further performed to explore the immune cell landscape of TNBC patients in the TCGA cohort. Analysis of immune cell infiltration patterns revealed substantial variations in the distribution of 22 distinct immune cell types in the high-risk group and the low-risk group (Fig. 6A). Such variations in the proportions of tumor-infiltrating immune cells might represent an intrinsic feature that could characterize individual differences. Moreover, the proportions of different subpopulations of tumor-infiltrating immune cells were weakly to moderately correlated (Fig. 6B). Seven immune infiltration algorithms were used to evaluate the association of five PCD-related genes and risk score with immune infiltration. The results indicated that risk score was positively correlated with immune infiltration; whereas COL5A3 and CHML were negatively correlated with immune infiltration (Figs. 6C, 6D). Box plots were applied to show the differences in drug sensitivity between the high-risk group and the low-risk group by using the oncoPredict package in R language. Low-risk group was more sensitive to most chemotherapy drugs, with respect to high-risk group, suggesting a potential factor for better prognosis in the low-risk group (Fig. 7A). Remarkably, the high-risk group was relatively sensitive to ABT737, Bl2536, and Daporinad, when compared to the low-risk group (Fig. 7B). This would provide a new way to study more effective chemotherapy regimen for high-risk TNBC patients.

Expression levels of PCD-related genes involved in prognostic model

Expression levels of the five PCD-related genes in the TNBC cell line and the mammary epithelial cell line were analyzed by PCR. Specifically, SCARB1, SEPTIN3, CHML, and SYNM were down-regulated, while COL5A3 was significantly upregulated in MDA-MB-231 compared with that in MCF-10A (Fig. 8). These results, along with the bioinformatics analysis data, demonstrated that the prognostic model is meaningful.

Construction of a nomogram

Univariate Cox regression analysis indicated that overall survival was significantly correlated with risk score, T, N, and stage in the TCGA cohort (Fig. 9A), and multivariate Cox regression analysis revealed that risk score and N stage were the independent prognostic factors (Fig. 9B). To build a useful predictive method, a prognostic nomogram that integrated independent prognostic factors, namely, PCD-related risk score and N stage, was generated (Fig. 9C). Based on the time-dependent ROC curves, the AUCs of nomogram achieved 0.99, 0.85, and 0.87 at 365 days, 1,095 days, and 1,825 days, respectively (Fig. 9D). Calibration curves showed the prediction value of the nomogram and demonstrated high accuracy of the predicted survival (Fig. 9E). These results demonstrated that the nomogram could be a promising method for predicting overall survival in TNBC patients.

Discussion

Breast cancer is the most common malignancy among women worldwide, and the subtype of the particular note is TNBC, which lacks expression of estrogen receptor, progesterone receptor, and human epidermal growth factor receptor-2. Both clinical and pathological features have been widely used to predict therapeutic response and long-term results. Since TNBC represents a heterogeneous group of diverse subtypes with its own biological and molecular characteristics (Jiang et al., 2019), identification of novel prognostic biomarkers and establishment of more accurate prognostic models remain imperative and warranted in TNBC research (Duffy et al., 2015).

PCD is a common gene regulated cell death mode in multicellular organisms. With the deepening of PCD research, various PCD modes are being discovered and studied (Wang et al., 2024; Chen et al., 2024). Recently, signatures based on the differentially expressed PCD-related genes have gained much attention and shown great potential in prognosis prediction of cancer (Chen et al., 2024; Cao et al., 2023; Dong et al., 2024; Cancer Genome Atlas Research Network, 2017; Cheng et al., 2023; Sha et al., 2022; Liu et al., 2023). Chen et al. (2024) screened 149 PCD-related differentially expressed genes, of which INHBA, LRRK2, HSP90AA1, HSPB8, and EIF2AK2 were identified as the hub genes of esophageal squamous cell carcinoma. Cao et al. (2023) constructed a 16 PCD-related gene model with potential in predicting prognosis and response to immune checkpoint inhibitors in cancer. Dong et al. (2024) identified seven PCD-related genes to establish the PCD-related risk score for the advanced non-small cell lung cancer model, effectively stratifying overall survival in patients with advanced non-small cell lung cancer. In the present work, a novel five-gene signature including SEPTIN3, SCARB1, CHML, SYNM, and COL5A3 was identified from a number of PCD-related genes and applied for prognosis prediction in TNBC patients. Specifically, SEPTIN3, SCARB1, CHML, and SYNM were associated with better prognosis; whereas COL5A3 was responsible for poor outcome. Risk score based on the five-gene signature was an independent prognostic factor of TNBC and patients in high-risk group showed a significantly worse overall survival, with respect to patients in the low-risk group. Results from survival curves and AUCs values indicated that the constructed PCD-related prognostic model was accurate for survival prediction in TNBC patients not only from the TCGA cohort but also from the GSE58812, GSE21653, and GSE135565 cohorts. A prognostic nomogram that integrated PCD-related risk score and clinical N stage was therefore generated. All these results demonstrated a good performance of PCD-related prognostic model and the prediction model could be a promising indicator for TNBC survival.

Gene SCARB1 encodes the scavenger receptor class B type I (SR-BI) glycoprotein and regulates cholesterol exchange between cells and high-density lipoproteins (Gutierrez-Pajares et al., 2016). Mounting studies have reported the significant role of SCARB1 gene and SR-BI protein in cancer proliferation and progression. In the context of breast cancer, over-expression of SR-BI could increase high-density lipoproteins-mediated proliferation of breast cancer cells via PI3K/AP-1 pathway (Cao et al., 2004). Down-regulation of SR-BI in breast cancer cells was associated with decreased cellular cholesterol content and reduced tumor aggressiveness (Danilo et al., 2013). COL5A3, a member of the collagen triple helical repeat family, takes part in cell growth and migration. Amrutkar et al. (2019) confirmed that COL5A3 promoted chemoresistance of pancreatic cancer cells to gemcitabine treatment. COL5A3 has also been shown to participate in breast cancer brain metastasis, assess the infiltration of cancer-associated fibroblasts, and predict immune and chemotherapy responses (Zhang et al., 2021; Song et al., 2024). CHML, also known as Rab escort protein 2, is one of the key factors for Rab proteins prenylation. Researchers have found that CHML was related to the development of urothelial carcinoma, multiple myeloma, hepatocellular carcinoma, and lung cancer (Li et al., 2008; Zhang et al., 2019; Chen et al., 2019; Dong et al., 2021). CHML could promote proliferation, inhibit apoptosis, and induce metastasis of tumor cells, and high expression of CHML is associated with poor survival (Li et al., 2008; Zhang et al., 2019; Chen et al., 2019; Dong et al., 2021). SEPTIN3, a member of the septin family primarily expressed in brain and testis, is a membrane-bound presynaptic protein connected to autophagy (Rosa et al., 2020). Wang, Yang & Gao (2023) have found that SEPTIN3 was overexpressed in TNBC and was related to poor prognosis (Yang, Wang & Gao, 2024). Moreover, SEPTIN3 was observed to favor cell growth and oncogenesis. SEPTIN3 promoted TNBC cell aggressiveness and proliferation via activation of Wnt signaling pathway (Wang, Yang & Gao, 2023). SYNM is a type IV intermediate filament that has been reported to modulate cell adhesion and motility. Aberrant promoter methylation of gene SYNM was associated with lymph node metastases and advanced tumor grade (Noetzel et al., 2010). Besides, SYNM was positively correlated with CD8 T cells and monocytes, but was negatively correlated with γδ T cells and M1 macrophages, suggesting that SYNM could affect breast cancer cells by modulating immune-infiltrating cells (Bao & He, 2022). Given that only a few publications report the roles of these genes on TNBC are currently available, more thorough researches are now needed to clarify the significance and mechanism of these five genes.

Besides resisting cell death and metastasis of malignant cells, tumor metabolism and immune microenvironment are two hallmarks of cancer (Hanahan, 2022). In the present study, KEGG pathway and GO enrichment analysis were carried out by evaluating the differentially expressed genes between high-risk and low-risk groups classified by the five-gene signature. Several metabolism-related pathways including protein and fat digestion and absorption, glycerolipid, glycerophospholipid, and arachidonic acid metabolism were significantly enriched. Moreover, immune infiltration analysis revealed substantial variations in the distribution of 22 distinct immune cell types in the high-risk group and the low-risk group. These gene enrichment analyses collectively suggested that PCD might be closely related to tumor metabolism and immune microenvironment of TNBC.

To the best of our knowledge, the five PCD-related gene signature and nomogram have not been reported previously and could be a promising prognostic model for TNBC patients. In particular, all these five genes, namely, SEPTIN3, SCARB1, CHML, SYNM, and COL5A3, were first applied in the prognosis prediction model for TNBC. The constructed nomogram combining the present PCD-related gene signature with patient’s N stage showed improved performance in predicting overall survival, and were relatively higher than some other published nomograms for TNBC prediction (Cheng et al., 2023; Sha et al., 2022). Our study has several limitations. First, the five PCD-related gene signature and prognostic model was established by analyzing retrospective data from public databases, leading to inevitable bias. Second, few publications reporting the functions and mechanisms of these genes on TNBC are currently available. Third, TNBC usually progresses rapidly, develops metastasis and relapses, so disease-free survival, progression-free survival, and short-term survival should be included as well. Fourth, genomic and transcriptomic landscape of TNBC were quite different (Jiang et al., 2019). Given that TNBC subtypes were not classified in the selected TCGA and GEO database, we are unable to classify the transcriptomic, clinicopathological, and survival information of TNBC patients in the present study. Fifth, the external validity of the prognostic model and its applicability to diverse clinical settings are not evaluated, and taking the aforementioned limitations together suggest that our findings need to be validated with prospective large cohort studies, and more thorough cellular research should be performed in this field.

Conclusion

In summary, this study established a new five-gene prognostic model and nomogram to predict overall survival of TNBC patients, which can not only be applied for patients management, but also provide a new direction for exploring therapeutic targets of TNBC.

Supplemental Information