Gene expression datasets were retrieved from the Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/geo/). Two independent microarray datasets were selected for this in silico analysis: GSE231629 (containing 109 samples) and GSE3744 (containing 47 samples). The selection of these specific datasets was based on three criteria: (1) both datasets were generated on the same Affymetrix GPL570 platform ([HG-U133_Plus_2] Array), eliminating cross-platform batch effects without requiring additional normalization; (2) each dataset provided paired or adjacent normal breast tissue controls, which are not universally available in larger repositories; and (3) the smaller sample size of GSE3744 served as a conservative validation set to ensure that identified DEGs were not artifacts of a single large cohort. We acknowledge that larger cohorts such as TCGA (n > 1000) and METABRIC (n ~ 2000) exist. However, in the present study, TCGA data were intentionally reserved for independent validation of expression patterns (via GEPIA/UALCAN) and survival analysis (via Kaplan-Meier Plotter, which incorporates TCGA and GEO datasets). This two-stage design (discovery in smaller, well-matched GEO cohorts followed by validation in large-scale repositories) minimizes the risk of overfitting and ensures generalizability. Since the present analyses were performed on microarray data that had been pre-processed and normalized by the GEO platform, and given the use of an identical platform for both datasets, no additional normalization or batch effect correction was applied. All databases were accessed between January and February 2025. To ensure the reproducibility and currency of our findings, we re-analyzed the key results using the most current versions of these platforms as of May 2026. No significant changes in the direction or magnitude of hazard ratios or DEG lists were observed, confirming the stability of our conclusions.

Differentially expressed genes between BC and adjacent normal tissues were identified separately for each dataset using the interactive web tool GEO2R (accessed January 2025). GEO2R employs the limma R package to perform linear model analysis on the normalized GEO series matrix data. Genes with an adjusted p-value (Benjamini & Hochberg false discovery rate) < 0.01 and |log2-fold change| (|log2FC|) > 114,15 were selected as DEGs.

The list of common DEGs present in both datasets was submitted to the STRING database (version 11.5; https://string-db.org) to construct a protein-protein interaction (PPI) network16. For the identification of significant interactions within the PPI network, a combined score > 0.4 from the STRING database was considered. This threshold, consistent with common practices in similar studies, indicates a medium to high level of confidence in the reliability of the predicted interactions16. The resulting network file was downloaded and imported into Cytoscape software (version 3.9.1) for visualization and further analysis. Hub genes within the PPI network were identified using the cytoHubba plugin (version 0.1). The top-ranking nodes from each of eight distinct topological algorithms (Maximum Clique Centrality [MCC], Maximum Neighborhood Component [MNC], Density of Maximum Neighborhood Component [DMNC], Degree, Closeness, Bottleneck, Eccentricity, and Edge Percolated Component [EPC]) were extracted. Genes that consistently appeared among the top 10 nodes across the majority of these algorithms were selected as candidate hub genes for downstream analysis.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for the common DEGs using the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8; https://david.ncifcrf.gov/). The analysis included Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) categories from GO. Terms with a Benjamini-Hochberg-adjusted p-value < 0.05 were considered significantly enriched.

The prognostic value of the identified hub genes was assessed using the Kaplan-Meier Plotter online tool (https://kmplot.com/analysis/, accessed February 2025, re-verified May 2026). The BC dataset was selected. Two complementary survival analyses were performed:

First, an unstratified analysis using the median cutoff option to dichotomize patients into high- and low-expression groups based on the probe with the most significant prognostic value (Jetset best probe). Hazard ratios (HR) with 95% confidence intervals and log-rank p-values were calculated and extracted directly from the tool.

Second, and more importantly, a subtype-stratified analysis was performed by splitting patients according to PAM50 molecular subtypes (basal-like, HER2-enriched, luminal A, luminal B) using the built-in subtype annotation in Kaplan-Meier Plotter. This stratification was conducted to investigate whether the prognostic effects of the hub genes were context-dependent and to resolve potential paradoxical findings from the unstratified analysis. Hazard ratios and p-values were calculated separately for each subtype.

The correlations of hub genes with clinicopathological parameters, such as tumor subtypes, were analyzed using UALCAN (https://ualcan.path.uab.edu/, accessed February 2025, re-verified May 2026), a portal for The Cancer Genome Atlas (TCGA) data.

Experimentally validated miRNA-mRNA interactions for the hub genes were retrieved from the ENCORI (starbase.sysu.edu.cn) and miRTargetLink 2.0 (https://ccb-compute.cs.uni-saarland.de/mirtargetlink2/) databases. The search criteria in ENCORI were set to: strict stringency (CLIP-Data >= 5), with only experimentally validated miRNA-mRNA interactions with strong evidence (supported by reporter assays, Western blots, or qPCR) being included.

Transcription factors (TFs) regulating the hub genes were predicted using the TRRUST database (version 2; https://www.grnpedia.org/trrust/). This database contains manually curated transcriptional regulatory networks for humans. All predicted TFs with documented regulatory relationships (activation or repression) were recorded17.

Potential pharmacological agents targeting the hub genes were identified using the Drug-Gene Interaction Database (DGIdb, version 4.2.0; http://www.dgidb.org). The search was conducted using default parameters, and all predicted interactions, along with their interaction scores and sources, were compiled for analysis.

Differential expression analysis of the two GEO datasets (GSE3744 and GSE231629) identified 2360 and 1140 DEGs, respectively, as shown in Figure 1A-B. A comparative analysis revealed 93 genes that were consistently differentially expressed in both datasets, as shown in Figure 1C. To elucidate the functional relationships among these common DEGs, a PPI network was constructed using the STRING database. The resulting network, comprising 93 nodes and 226 edges, was imported into Cytoscape for topological analysis (Supplementary material S1). Hub genes within this network were identified by integrating results from eight topological algorithms (MCC, MNC, DMNC, Degree, Closeness, Bottleneck, ECC, and EPC) using the cytoHubba plugin. An UpSet plot analysis was employed to evaluate consensus, revealing seven key hub genes: CDK1, CDC20, CDCA8, and CENPF, which showed consensus across all eight algorithms, along with TTK, FOXM1, and MYBL2, whose results are detailed in Figure 2. These hub genes exhibited the highest degree of connectivity within the PPI network, suggesting their central regulatory roles.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on the 93 common DEGs to determine their biological functions. The GO analysis revealed significant enrichment in biological processes related to cell division and chromosome segregation. Specifically, the most significantly enriched BP terms included mitotic nuclear division, chromosome segregation, and cell cycle G2/M phase transition. For CC, the DEGs were predominantly localized to the condensed chromosome, centromeric region, and spindle. Key MF terms included microtubule binding and ATPase activity, with the results presented in Table 1. KEGG pathway analysis further corroborated these findings, indicating that the common DEGs were significantly enriched in four primary pathways, including the cell cycle, oocyte meiosis, cellular senescence, and viral carcinogenesis, as shown in Table 2. The cell cycle pathway demonstrated the highest level of enrichment, encompassing multiple hub genes (CDK1, CDC20, and CDCA8).

The mRNA expression levels of the seven hub genes were analyzed using the GEPIA database, which incorporates data from the TCGA and GTEx projects. Consistent with the initial DEG analysis, all seven hub genes were significantly overexpressed in breast invasive carcinoma (BRCA) tumor tissues compared to normal breast tissues, and the results are shown in Figure 3A-G (p < 0.05). The most pronounced overexpression was observed for CDCA8 and CENPF, while TTK showed the most modest increase.

To evaluate their prognostic value, overall survival (OS) analysis was conducted using the Kaplan-Meier Plotter database. Initial unstratified analysis showed that high mRNA expression of CDCA8, MYBL2, and CENPF was significantly associated with poorer OS (Hazard Ratio [HR] > 1, log-rank p < 0.001). Conversely, high expression of CDC20, FOXM1, CDK1, and TTK genes showed a trend toward better OS (HR < 1, log-rank p < 0.05), which appeared paradoxical given their well-established oncogenic roles in promoting cell cycle progression.

To resolve this paradox and rigorously explore potential confounding by molecular heterogeneity, we performed stratified survival analysis based on PAM50 molecular subtypes (basal-like, HER2-enriched, luminal A, luminal B). This analysis revealed that the prognostic impact of these hub genes is highly context-dependent (Table 3, Figure 4A-G).

Results from the subtype-stratified analysis:

These results demonstrate that failing to account for molecular subtype heterogeneity can lead to misleading or even opposite prognostic associations. The paradoxical protective effects observed in the unstratified analysis (e.g., FOXM1, CDC20) were resolved by subtype stratification, revealing context-dependent risk associations that align with known biology. Further details are presented in Supplementary material S2.

The relationship between hub gene expression and BC molecular subtypes was investigated using the UALCAN portal with TCGA data. All seven hub genes displayed significantly elevated expression in triple-negative breast cancer (TNBC) samples compared to other subtypes (Luminal A, Luminal B, HER2-enriched), as illustrated in Figure 5A-G (p < 0.05). This suggests a potential subtype-specific role for these genes, particularly in the more aggressive TNBC subgroup.

To investigate the post-transcriptional regulation of the hub genes, miRNA-mRNA interactions were analyzed using two authoritative databases, ENCORI and miRTargetLink 2.0, considering only interactions with strong experimental evidence. Analysis of these data revealed a complex regulatory landscape, which was visualized as an interaction network (Supplementary material S3). In total, 20 unique miRNAs were identified that simultaneously targeted three or more of the hub genes, as shown in Table 4. Among these, hsa-miR-29a-3p and hsa-miR-29b-3p emerged as key regulators, predicted to target five of the seven hub genes (MYBL2, CDK1, CDC20, CDCA8, and CENPF).

Two hub genes (FOXM1 and MYBL2) are transcription factors that may act as important regulators in BC patients. FOXM1 targets 13 genes (APOE, AR, BRIP1, BTG2, CCNB1, CDC25A, CDC25B, CDC6, CDH1, CDKN2A, CYP3A4, FGB, and KRT15), and MYBL2 targets 8 genes (CCNA1, CCND1, MYBL2, MYC, COL1A1, PSG1, NCL, and SP1). The transcription factors that may regulate the expression of these target genes were identified based on the TRRUST database (Supplementary Table S4). Furthermore, our analysis revealed that three other hub genes (CDK1, CDC20, and TTK), while not themselves transcription factors, are under strong regulation by upstream transcription factors, highlighting their pivotal position within the regulatory network.

Potential drug compounds were identified and ranked according to the highest aggregate interaction score provided by the DGIdb database. A screening for potential pharmacological interventions was performed using the DGIdb database. Among the seven hub genes, only CDK1 and TTK had documented interactions with known drugs or investigational compounds in the database, as shown in Table 5. For CDK1, 22 potential inhibitors were listed, including palbociclib and ribociclib (CDK4/6 inhibitors with known CDK1 affinity). For TTK (also known as MPS1), three investigational kinase inhibitors (BAY-1217389, CFI-402257, and MPS1-IN-1) were identified. No direct, high-confidence drug interactions were reported for the remaining five hub genes in the queried database.

BC: Breast cancer; BP: Biological Process; CC: Cellular Component; DEGs: Differentially expressed genes; ER: Estrogen receptor; GEO: Gene Expression Omnibus; GO: Gene Ontology; HER2: Human epidermal growth factor receptor 2; HR: Hazard Ratio; KEGG: Kyoto Encyclopedia of Genes and Genomes; MF: Molecular Function; OS: Overall survival; PAM50: Prediction of Analysis of Microarray 50 molecular subtype classifier; PPI: Protein-protein interaction; PR: Progesterone receptor; TFs: Transcription factors; TNBC: Triple-negative breast cancer; WGCNA: Weighted Gene Co-expression Network Analysis.

The paper has been extracted from Sedigheh Akhtartavan's PhD thesis supported by the Research Council of Semnan University of Medical Sciences (2083).

Not applicable.

This study was supported by the Research Council of Semnan University of Medical Sciences (Grant No. 2083).

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

Not applicable.

Not applicable.

None.

The authors declare that they have no competing interests.