We applied the common datasets GSE5364 (23), GSE29044 (24), and GSE42568 (25) from the Gene Expression Omnibus (GEO https://www.ncbi.nlm.nih.gov/geo/) database, and obtained gene expression data and corresponding clinical data from the TCGA-BRCA project in The Cancer Genome Atlas (TCGA https://portal.gdc.cancer.gov/) database.

We searched and extracted mitophagy-related genes from the GeneCards database (https://www.genecards.org/) using the keyword “mitophagy.” The GeneCards relevance score integrates evidence from curated databases, functional annotations, and published literature to quantify gene–process associations. Genes with a relevance score > 1 were selected to retain those with moderate to strong evidence supporting their involvement in mitophagy, as commonly adopted in previous bioinformatics studies, while reducing potential noise from weakly associated genes. This threshold allowed us to maintain sufficient gene coverage for network construction and survival analyses while minimizing the inclusion of weakly supported candidates.

For three data sets from the GEO database (GSE5364, GSE29044, GSE42568), we used GEO2R to identify differentially expressed genes and calculated the significance FDR for each gene using the P.adjust function. Finally, based on the significant differential genes (adj.P<0.05, Log Fold Change>1.5) and the availability of data, we achieved visualization of volcano and box plots.

We selected significant differentially expressed genes from each of the three datasets in the GEO database. Use the Metascape database (https://metascape.org/) to perform correlation and enrichment analysis on three groups of significantly differentially expressed genes. We used the following criteria for the analysis: (1) Three gene lists were used to input genes in sequence; (2) For GO annotations, three main types of gene functions were extracted: biological processes (BP), molecular functions (MF), and cellular components (CC); (3) In functional enrichment analysis, the default parameters for filtering are: minimum overlap of 3, enrichment factor of 1.5, and P-value of 0.01. Meanwhile, the Metascape database also constructed and visualized protein-protein interaction networks at input sites based on STRING software.

Firstly, we extracted significantly differentially expressed genes (DEGs) that were consistently identified across three GEO datasets (GSE5364, GSE29044, and GSE42568), resulting in a total of 177 common DEGs. These genes represent robust BC–associated transcriptional alterations shared among independent cohorts. Subsequently, the 177 DEGs were intersected with mitophagy-related genes retrieved from the GeneCards database to identify genes that are both functionally associated with mitophagy and transcriptionally dysregulated in BC. This integrative filtering strategy yielded 70 overlapping genes, which were defined as mitophagy-related genes in the context of BC and were subjected to subsequent analyses in this study. After that, import 70 intersection sites into the STRING database (https://string-db.org) and generate a PPI network. When constructing a PPI network, PPI network pairs with a moderate confidence score > 0.4 were exported, and then visualized using Cytoscape (version 3.8.2). Among them, we applied Molecular Complex Detection (MCODE http://apps.cytoscape.org/apps/mcode) to find densely connected areas and subnets within the main network; In addition, five topological analysis methods of cytoHubba (https://apps.cytoscape.org/apps/cytohubba) were applied to rank and screen hub sites using gene attributes.

We use the Sangerbox database (http://www.sangerbox.com/tool) to obtain somatic mutation data related to BC samples in “maf” format from the TCGA database. We used the “Maftools” software package in R software to draw the waterfall diagram of somatic mutation and visualize the mutation sites.

In order to screen potential prognostic related sites for BC among 70 Mitophagy related sites, we extracted 12 sites from the “MCODE1” network with the highest score. Use the Kaplan Meier Plotter database (http://kmplot.com/analysis/) to draw the K-M survival curve, and the ROC Plotter database (www.rocplot.org) to draw the ROC curve; Analyze each site and further screen 12 Mitophagy related sites based on K-M survival data and ROC-AUC (26).

In the Sangerbox database, we identified Mitophagy-related tumor subtypes (High-Risk, Low-Risk) by consistency clustering using the “ConcentrusClusterPlus” software package from R software. Then, we evaluated the clustering results with k between 2 and 10 and repeated the process of obtaining the ideal number of clusters 1000 times to ensure stable results and achieved data visualization.

We applied Gene Set Enrichment Analysis (GSEA) to detect whether the enrichment of GO and KEGG pathways is statistically significant in high-risk and low-risk subtypes. The transcriptome data from TCGA database were divided into two groups in high-risk and low-risk subtypes and imported into GSAE software. Both p<0.05 and FDR<0.25 are considered criteria for significant gene sets (27).

In order to identify the difference of immune characteristics between high risk and low risk subtypes, we loaded the transcriptome expression data into CIBERSORT database (https://cibersort.stanford.edu/), and repeated the results 1000 times to determine the relative percentage of 22 immune cell types. To further confirm the accuracy of the data, we used the MCPCounter Package from Rstudio to detect infiltration differences in 8 immune cell families and two stromal cell families in two tumor subtypes. Finally, we analyzed and calculated the differences in the scores of Stromal score, Immune score, and ESTIMATE score between the two subtypes using the ESTIMATE database (https://bioinformatics.mdanderson.org/estimate/).

The immunohistochemical (IHC) images of normal human tissues and BC tissues were obtained from human protein profiles (HPA, https://www.proteinatlas.org/). Used to identify the protein expression differences of the ultimately screened 8 Mitophagy-related sites between human normal tissue samples and BC samples.

The NCI-60 cell line panel is a drug efficacy screening developed by the Developmental Therapeutics Program (DTP) of the National Cancer Institute (NCI) in the United States. Thousands of compounds have been applied to the NCI-60. CellMiner (https://www.nci.nih.gov/cellminer/) is a web-based genomic and pharmacological toolset for exploring transcriptional and drug profiles in the NCI-60 cell line panel. Based on the corresponding data from CellMiner, the Corrplot R software package and Spearman method (p < 0.05) were used to analyze the relationship between the expression of PBK, NEK2 and drug sensitivity.

MDA-MB-231 and MCF-7 was purchased from the American Type Culture Collection (Manassas, VA, USA), and has been authenticated for mycoplasma contamination and short tandem repeat profile analysis. In an incubator set at 37 °C with 5% CO₂, cells were grown in DMEM combined with 10% fetal bovine serum and 100 U/ml penicillin-streptomycin combination (Gibco, USA).

Antibodies against PBK(4942T), a-Tublin(3873T), b-Actin(4967T) were acquired from Cell Signaling Technology (Danvers, MA, USA); Antibodies against NEK2(ab279717) were acquired from Abcam (Cambridge, MA, USA). Ki67(9129S) was acquired from Cell Signaling Technology (Danvers, MA, USA).

The cell lysis assay kit (KGP 2100) was purchased from Keygen Biotech (Nanjing, China). Cells were added to lysis buffer and lysed on ice for 5 minutes. The cells were scraped off using a cell suspension apparatus and transferred to a 1.5 ml centrifuge at 4°C for 15 minutes at 13000 rpm. The supernatant was tested for protein concentration using a detergent compatible Bradford protein assay kit (P0006 C). For protein extraction of liver cancer tissue and xenograft tumor tissue, 300 μ l of whole protein lysate was added every 10 mg of tissue, and thoroughly ground in a protein homogenizer. The supernatant was centrifuged under the same conditions as described above to detect protein concentration. The proteins were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis, and then moved to the polyvinylidene fluoride (PVDF) membrane. Prepare PVDF membrane by blocking it with 5% BSA, and then add primary antibody overnight at 4°C; Add secondary antibody at room temperature for 1 hour, scan with Li COR Odyssey infrared imaging system, and finally analyze with Image J software.

5x10⁴ cells were cultured in 96-well cluster plates. CCK-8 reagents were used to dilute at a ratio of 1:9; add CCK-8 reagent 100 μl per well in the dark for 2 h in a 37°C, finally read the absorbance at 450 nm wavelength.

Cells were seeded in confocal dishes at a density of 1 × 10⁵ cells per well and treated with the indicated drugs or reagents. After treatment, cells were fixed with 100% methanol at −20°C for 15 min and permeabilized with 0.1% Triton X-100 for 5 min. Following PBS washes, cells were blocked with 5% bovine serum albumin (BSA) in PBS at room temperature for 1 h. Subsequently, cells were incubated with primary antibodies at 4°C overnight. After thorough PBS washing, cells were incubated with fluorescently labeled secondary antibodies at room temperature for 2 h. Nuclei were counterstained with DAPI. Fluorescence images were captured using an Olympus FV1000 confocal laser scanning microscope (Tokyo, Japan). The degree of protein colocalization was quantified using Image-Pro Plus software, which calculated the overlap coefficient between the two fluorescence signals.

All statistical analyses were performed using SPSS software (version 21.0; IBM Corp., Armonk, NY, USA). Continuous variables were presented as mean ± standard deviation or median (interquartile range), as appropriate, and compared using the Student’s t-test or Mann–Whitney U test. Categorical variables were analyzed using the chi-square test or Fisher’s exact test. Survival differences were evaluated using Kaplan–Meier analysis and compared by the log-rank test. Univariate and multivariate Cox proportional hazards regression analyses were conducted to identify independent prognostic factors, with hazard ratios (HRs) and 95% confidence intervals (CIs) reported. Logistic regression analysis was used to assess the independent association of PBK and NEK2 expression with high-risk classification, and odds ratios (ORs) with 95% CIs were calculated. All statistical tests were two-sided, and a P value < 0.05 was considered statistically significant.

Firstly, three gene expression programs (GSE5364, GSE29044, and GSE42568) from the GEO database were included in our study. We used the website analysis tool GEO2R from GEO database to analyze the transcriptome expression profile data of three datasets to screen DEGs (Log Fold-Change ≥ 1.5 or Log Fold-Change ≤ -1.5 and adj.P<0.05). Among them, in GSE5364, we obtained 687 DEDs, including 401 upregulated sites and 286 downregulated sites; In GSE29044, we obtained 723 DEDs, including 232 upregulated sites and 491 downregulated sites; In GSE42568, we obtained 2049 DEDs, including 966 upregulated sites and 1083 downregulated sites. The DEGs of the three datasets are shown in Figures 1A–C by volcanic maps, respectively; The data availability and standardization results are shown in the box plot in Figures 1D–F. Subsequently, we intersected three sets of DEGs and obtained 177 common differentially expressed genes (Figure 1G). Finally, we downloaded mitophagy-related sites from the GeneCards database, with the inclusion criteria of reliability score>1. And by intersecting 177 common differentially expressed genes with all obtained mitophagy-related sites, 70 mitophagy related sites were ultimately obtained (Figure 1H).

In order to understand the potential biological effects and site interactions of DEGs. We used the Metascape database to conduct pathway enrichment and correlation analysis on DEGs. The pathway enrichment analysis in the Metascape database shows that these DEGs are mainly enriched in the following pathways: circulatory system process, negative regulation of cell population proliferation, negative regulation of cell differentiation, cell population proliferation, positive regulation of transferase activity, regulation of epithelial cell proliferation, response to growth factor, positive regulation of locomotion, positive regulation of locomotion, enzyme-linked receptor protein signaling pathway, response to hormone, tube morphogenesis, actin filament-based process from GO terms; and Signaling by Rho GTPases, Signaling by Nuclear Receptors, Signaling by Receptor Tyrosine Kinases, Extracellular matrix organization from Reactome dataset (Figures 2A, B). Afterwards, we defined the interaction and similarity relationship of DEGs using the Overlap diagram (Figures 2C, D), and the results showed a significant correlation between the significant sites in the three BC datasets. Finally, we constructed a DEGs interaction network through three dimensions (enrichment pathway level, dataset level, and P-value level) to better understand the potential interaction relationships of DEGs (Figures 2E–G).

To further elucidate the regulation and interaction relationships of 70 mitophagy-related sites, we imported the obtained 70 mitophagy-related sites into the STRING database for protein-protein interaction analysis, and visualized the results in Cytoscape (version 3.8.2). The STRING database identified 70 protein products and 376 edges PPI networks (Figure 3A). After that, we used the MCODE plugin in Cytoscape to identify two important interaction modules in the interaction network: MCODE1 and MCODE2 (Figures 3B, C). Among them, MCODE1 has obvious interactions, and through the prediction of interaction networks and subnetworks, MCODE1 is believed to play an important role in the development of BC. In addition, we used the cytohubba plugin. The top ten hub genes obtained through five algorithms of MCC, DMNC, MNC, Degree and EPC were shown in Supplementary Table 1, and the overlapping hub genes in the five algorithms are verified through the Venn diagram (Figure 3D).

We included 70 mitophagy related sites in the study of somatic mutation. The results showed that in the mutation spectrum data of BC, the mutation degree of 70 mitophagy related sites varied from 9.5% to 0.6% (Figure 4), which may be a potential factor leading to the dual role of mitophagy in the prognosis of BC.

The PPI network shows that the interaction between MCODE1 is significant; And through the prediction of MCODE interaction networks and subnetworks, MCODE1 is believed to play an important role in the development of BC. To verify the impact of MCODE1 sites on the prognosis of BC patients and extract potential prognostic biomarkers, we extracted 12 sites from MCODE1 (CDC20, TPX2, TTK, PBK, TOP2A, TK1, NEK2, FOXM1, CEP55, CDK1, CCNB2, BIRC5), and analyzed the prognostic role of each locus in BC based on the Kaplan Meier Plotter database; The inclusion criteria for effective effects are: confidence interval (CI) excluding 1 and logrank-P<0.05. The results showed that 12 sites met the inclusion criteria and were considered valuable for further ROC validation (Figures 5A–L).

For the reason to further confirm the diagnostic and therapeutic efficacy of the 12 sites (CDC20, TPX2, TTK, PBK, TOP2A, TK1, NEK2, FOXM1, CEP55, CDK1, CCNB2, BIRC5) selected in “3.5” and build a disease diagnosis model based on mitophagy related biomarkers in BC, we use the transcriptome data of BC in the ROC Plotter database, at the same time link gene expression levels with response to therapy, and analyze the prognostic effects of diagnosis and treatment at each site to ultimately determine prognostic biomarkers. The inclusion criteria we adopted were: AUC>0.5, P value<0.05. The final analysis results indicate that 9 sites (CDC20, TPX2, PBK, TOP2A, NEK2, FOXM1, CDK1, CCNB2, CEP55) have potential roles in the diagnosis and treatment of BC, and have the value of becoming biomarkers for BC prognosis (Figures 6A–L). Finally, we analyzed the expression of 9 prognostic related sites in control tissue, BC tissue, and BC metastatic tissue, and found that their expression levels in BC tissue were significantly higher than those in control tissue (Figures 6M, N).

Firstly, we analyzed the expression differences between BC and control groups for 70 mitophagy related sites based on the transcriptome gene expression profile of TCGA-BRCA project, and the difference results were shown in the form of heatmap (Figure 7A). After analyzing the complete body differences, we used Consistency clustering analysis to determine the mitophagy-related clusters of BC samples. After conducting k-means clustering analysis, two subgroups (C1, C2) in the TCGA queue were identified as having different mitophagy gene expression patterns (Figures 7B–E). Subsequently, in order to verify the survival status of patients in the C1 and C2 subgroups, we downloaded survival data from the TCGA database for patients in two subgroups and conducted K-M survival analysis on both subgroups. The survival curve showed that the survival level of the C2 subgroup was significantly higher than that of the C1 subgroup (Figure 7G). Therefore, we define the C1 subgroup as a high-risk group and the C2 subgroup as a low-risk group. Finally, we analyzed the expression of 70 mitophagy-related sites in the high-risk and low-risk groups, and found significant differences in expression between the two groups (Figure 7F). There was a significant block clustering in the heatmap.

To further determine the potential functional differences caused by differential expression sites in high-risk groups and low-risk groups, we performed a GSEA analysis based on different pathway sets (GO and KEGG) for the transcriptome expression data of the two subgroups. The analysis results show that in GO terms, the condensed chromosome, the condensed chromosome centromeric region, the meiosis i cell cycle process, the chromosome centromeric region, nuclear chromosome segregation are the five major enrichment pathways for high-risk groups (Figure 8A); endocrine process, the regulation of osteoblast differentiation, artery development, the positive regulation of osteoblast differentiation, the renal system process are the five major enrichment pathways in the low-risk group (Figure 8B). In KEGG terms, spliceosome, cell cycle, proteasome, RNA degradation, oocyte meiosis are the five major enrichment pathways of high-risk groups (Figure 8C); focal adhesion, ecm receptor interaction, complement and coagulation cascades, hypertrophic cardiomyopathy hcm, dilated cardiomyopathy are the five major enrichment pathways in the low-risk group (Figure 8D).

We further analyzed the tumor microenvironment and immune invasion level of different subgroups. Firstly, we used the ESTIMATE tool to score the immune levels of the high-risk and low-risk groups. The results showed that the Stromal Score, Immune Score, and ESTIMATE-Score of the low-risk group were higher than those of the high-risk group, and the results were statistically significant (Figure 9A). Subsequently, we used the CIBERSORT tool to analyze and validate the infiltration levels of 22 immune cells between the high-risk and low-risk groups. The results showed that, the infiltration level of: cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), NK cells resting, Macrophages M0, Macrophages M1, Dendritic cells increased in the high-risk group; The infiltration level of: B cells native, T cells CD4 memory resting, Monocytes, Dendritic cells resting, Mast cells resting increased in the low-risk group; And the results are all statistically significant. The analysis results are presented in stacked bar charts and box plots (Figures 9B, C). Finally, we also used MCPCounter tool for immune infiltration analysis between high-risk and low-risk groups, and the analysis results showed that, the infiltration level of Monocytic lineage increased in the high-risk group; The infiltration level of: B lineage, Myeloid dendritic cells, Neutrophils, Endothelial cells and Fibroblasts increased in the low-risk group; and the results are statistically significant (Figure 9D).

We obtained IHC images of the 9 prognostic related sites (CDC20, TPX2, PBK, TOP2A, NEK2, FOXM1, CDK1, CCNB2, CEP55) screened in the Human Protein Atlas database, and further investigated their protein expression differences in normal and BC tissues (Figure 10). The results showed that the protein expression of 9 prognostic related sites in BC tissues was higher than that in normal tissues, which was consistent with the above analysis results of transcriptome data.

To further validate the biological functions of the mitophagy-related prognostic biomarkers identified above, we selected PBK and NEK2, two genes that exhibited the most significant prognostic value and expression differences, for functional experiments in breast cancer cells (Triple-negative breast cancer cells MDA-MB-231). Western blot analysis confirmed that siRNA-mediated knockdown of PBK and NEK2 markedly reduced their protein expression levels compared with the control group (Figures 11A–D). Subsequently, CCK-8 assays demonstrated that silencing either PBK or NEK2 significantly inhibited the proliferation of breast cancer cells (Figures 11E, F). Moreover, confocal microscopy analysis revealed that the proliferation activity and nuclear density of breast cancer cells were notably decreased after PBK or NEK2 interference (Figures 11G, H). PBK is negatively correlated with the compound Lificguat and positively correlated with the compound GW-5074 and BMS-777607 (P<0.05), while NEK2 is negatively correlated with the compound Lexibulin and positively correlated with the drug Dexrazoxane (P<0.05) (Figure 11I). Furthermore, in order to further verify the stability of the results, we further verified the functions of PBK and NEK2 in the Non triple negative breast cancer cell line MCF-7, and got the same results (Figures 12A–H).

Clinicopathological characteristics of patients in the high- and low-risk groups are summarized in Supplementary Table 2, with categorical variables compared using the χ² test. Survival status did not differ significantly between groups (P = 0.512). Significant differences were observed in T stage (P < 0.001), M stage (P = 0.002), and clinical stage (P < 0.001). The high-risk group showed a higher proportion of T2 tumors, more M0 cases, and increased proportions of stage IIA and IIB, while fewer early-stage tumors were observed, indicating more advanced tumor characteristics. N stage showed a marginal difference (P = 0.051). Continuous variables were analyzed using independent-samples t tests (Supplementary Table 3). Patients in the high-risk group were slightly younger than those in the low-risk group (56.61 vs. 58.74 years; mean difference = −2.13; P = 0.009). In addition, PBK and NEK2 expression levels were significantly higher in the high-risk group (both P < 0.001). Logistic regression analyses were performed to assess risk-associated factors (Supplementary Table 4). In univariate analysis, PBK and NEK2 were strongly associated with high-risk status (both P < 0.001). After mutual adjustment, both PBK (OR = 3.23, 95% CI: 2.43–4.28) and NEK2 (OR = 3.37, 95% CI: 2.52–4.51) remained independently associated with high-risk status, supporting their roles as robust contributors to risk stratification.