Datasets: Publicly available transcriptomic and clinical data were obtained from The Cancer Genome Atlas (TCGA-BRCA, n=1,098 tumor/113 normal) and Gene Expression Omnibus (GSE42568: 121 tumor/17 normal; GSE9893: 155 tumor; GSE7390: 198 tumor). The METABRIC cohort (n=1,904) served as an independent validation set. Raw RNA-seq and microarray data were normalized using DESeq2 (TCGA) and limma (GEO), with batch effects corrected via ComBat. Clinical endpoints included overall survival (OS) and relapse-free survival (RFS).

Batch Effect Correction: To address technical variations across multiple datasets, we performed systematic batch effect correction. All datasets were individually normalized using z-score scaling, followed by identification of common genes across platforms. Batch effects were removed using the ComBat algorithm from the sva R package, with each dataset treated as a separate batch. Correction effectiveness was validated through principal component analysis (Supplementary Figure S1A). All subsequent analyses were conducted using the batch-corrected data.

Inclusion Criteria: (1) Histologically confirmed breast cancer; (2) Complete survival data; (3) RNA sequencing/microarray available.

Differential Expression Analysis: DEGs between tumor and normal tissues were first identified, utilizing DESeq2 (|log2FC| > 1, FDR < 0.05) for the TCGA cohort and limma (adj. p < 0.05) for the GEO dataset. Subsequently, these DEGs were subjected to Weighted Gene Co-expression Network Analysis (WGCNA) to identify clinically relevant modules. Using a soft-thresholding power of β = 6 to achieve a scale-free topology (R² > 0.85), we constructed co-expression networks and identified gene modules via dynamic tree cutting. The turquoise and blue modules, which showed the strongest correlations with clinical outcomes, were prioritized. Hub genes were rigorously defined within these modules by applying dual thresholds of module membership (MM) > 0.8 and gene significance (GS) > 0.2, culminating in the selection of 102 high-confidence hub genes for the construction of the prognostic model.

Feature Selection: The intersection of DEGs and WGCNA hub genes yielded 102 candidate genes, which were further refined to 51 genes shared across all cohorts. Algorithm Optimization: Six machine learning methods (StepCox, LASSO, Ridge, Elastic Net, RSF, SVM) were evaluated using 10-fold cross-validation. The StepCox-RSF hybrid model was ultimately selected based on its superior performance in integrated Brier scores and C-index. Random Survival Forest (RSF) Modeling: The RSF model was implemented using the randomForestSRC R package. Key hyperparameters, including the number of trees (ntree = 1000), the minimum node size (nodesize = 5), and the number of variables tried at each split (mtry = 2), were optimized via 10-fold cross-validation with the goal of minimizing the C-index error. Variable importance (VIMP) and minimal depth analysis were then applied to the optimized RSF model to identify the final 3-gene signature. Model performance was rigorously assessed using time-dependent ROC analysis (timeROC package) and Kaplan-Meier survival analysis.

Validation Cohorts: Prognostic scores were calculated as: Score = Σ(Genei Expression × RSF Coefficienti) Risk stratification used median cutoffs. Immune Microenvironment Analysis: CIBERSORT: Estimated immune cell fractions from bulk RNA-seq (LM22 signature matrix). ESTIMATE: Computed stromal/immune scores using estimate R package. Checkpoint Molecules: PD-L1 (CD274), LAG3, TIGIT expression quantified via normalized TPM/RSEM. Functional Enrichment: DEGs between subtypes underwent GO and KEGG analysis using clusterProfiler (FDR <0.05).

Single-Cell Data Processing: Public scRNA-seq datasets (GSE176078, GSE161529) were analyzed using Seurat (v4.0). Cells were clustered (resolution=0.8) and annotated via marker genes. Gene expression patterns were visualized using UMAP.

Murine breast cancer cell line 4T1 (American Type Culture Collection, ATCC), were maintained in Dulbecco’s Modified Eagle Medium (DMEM; Gibco) supplemented with 10% fetal bovine serum (FBS; HyClone) and 1% penicillin-streptomycin (Sigma-Aldrich). Cells were incubated at 37 °C in a humidified atmosphere containing 5% CO₂ and passaged at 80–90% confluence using 0.25% trypsin-EDTA (Gibco). Cell line authenticity was verified by short tandem repeat (STR) profiling, and mycoplasma contamination was routinely tested using a MycoAlert Kit (Lonza).

For stable EGR3 overexpression, lentiviral particles were generated by co-transfecting HEK293T cells with the pLVX-EGR3 plasmid (Clontech) and packaging plasmids (psPAX2 and pMD2.G) using Lipofectamine 3000 (Invitrogen). Viral supernatants were harvested at 48 and 72 hours post-transfection, concentrated via centrifugation (50,000 × g, 4 °C, 2 hours), and titrated using a Lenti-X qRT-PCR Titration Kit (Takara). Target cells (4T1) were transduced at a multiplicity of infection (MOI) of 10 in the presence of 8 μg/mL polybrene (Sigma-Aldrich). Stable clones were selected with 2 μg/mL puromycin (InvivoGen) for 7 days, and overexpression efficiency was confirmed by quantitative RT-PCR and Western blotting.

Cells (2 × 10³/well) were seeded in 96-well plates. Proliferation was assessed at 48 and 72 hours using a Cell Counting Kit-8 (CCK-8; Dojindo) following the manufacturer’s protocol. Briefly, 10 μL CCK-8 reagent was added to each well, incubated for 2 hours at 37 °C, and absorbance was measured at 450 nm using a microplate reader (BioTek).

Cell migration was evaluated using Transwell chambers (Corning) with 8-μm pores. Serum-starved cells (5 × 10⁴ in 200 μL serum-free DMEM) were added to the upper chamber, while the lower chamber contained DMEM with 20% FBS as a chemoattractant. After 24 hours, non-migrated cells were removed with a cotton swab. Migrated cells on the lower membrane surface were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and imaged under a light microscope (Nikon Eclipse). Cell numbers were quantified using ImageJ software (NIH).

All animal procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of Chongqing Medical university and conducted in compliance with ARRIVE guidelines. Female BALB/c mice (6–8 weeks old; Charles River Laboratories) were randomized into two groups (n = 6/group). For tumor induction, 4T1 cells (1 × 10⁶ cells in 100 μL PBS) transfected with EGR3 or empty vector were injected into the fourth mammary fat pad. Tumor dimensions (length and width) were measured every 3 days using calipers, and tumor volume was calculated as: Volume=0.5×length×width². Mice were euthanized at day 21 post-injection, and tumors were excised, weighed, and processed for downstream analyses.

Peripheral blood samples were collected from healthy adult volunteers under an institution-approved protocol. PBMCs were isolated using Ficoll-Paque PLUS (Cytiva) density gradient centrifugation according to the manufacturer’s instructions. Briefly, fresh blood was diluted 1:1 with phosphate-buffered saline (PBS). The diluted blood was carefully layered over an equal volume of Ficoll-Paque solution and centrifuged at 400 × g for 30 minutes at room temperature with the brake off. The PBMC layer at the interface was carefully aspirated, washed twice with PBS, and resuspended in complete RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin. Cell viability, assessed by Trypan Blue exclusion, was consistently >95%.

To investigate the immunomodulatory effect of EGR3-expressing tumor cells on T cell function within a more physiologically relevant immune context, we established a co-culture system using human PBMCs. Control (Vector) and EGR3-overexpressing (EGR3-OE) MDA-MB-231 human triple-negative breast cancer cells were seeded in 12-well plates and allowed to adhere overnight. The following day, freshly isolated PBMCs from healthy donors were added to the tumor cells at a tumor cell to PBMC ratio of 1:10. The co-culture was maintained in complete RPMI-1640 medium for 48 hours. PBMCs cultured alone in complete medium served as a baseline control. To evaluate T cell cytokine production, brefeldin A (1:1000, BioLegend) was added to the culture for the final 4–6 hours to inhibit protein transport, allowing for intracellular accumulation of cytokines. After co-culture, cells were harvested, stained for surface markers (CD8), fixed, permeabilized, and then stained intracellularly for IFN-γ and TNF-α prior to flow cytometry analysis of CD8⁺ T cell function.

Tumors were minced and enzymatically dissociated into single-cell suspensions using a Tumor Dissociation Kit (Miltenyi Biotec) and the gentleMACS Octo Dissociator (Miltenyi Biotec). Red blood cells were lysed with ACK buffer (Gibco). Cells were stained with fluorochrome-conjugated antibodies against CD8a (APC, clone 53-6.7), IFN-γ (PE, clone XMG1.2), and TNF-α (FITC, clone MP6-XT22; all from BioLegend) for 30 minutes at 4 °C in the dark. After washing, data were acquired on a BD FACSymphony™ flow cytometer (BD Biosciences) and analyzed using FlowJo v10.8 software (TreeStar). Gating strategies excluded doublets and dead cells using 7-AAD viability dye.

Data from at least three independent experiments are presented as mean ± SEM. Statistical significance was determined using two-tailed Student’s t-tests for comparisons between two groups. For comparisons across multiple groups, one-way ANOVA was performed, followed by Tukey’s honest significant difference (HSD) test for post-hoc analysis. In the survival analysis, the Benjamini-Hochberg procedure was applied to control the false discovery rate (FDR) for multiple comparisons across genes. All statistical tests were two-sided, and a p-value of less than 0.05 was considered statistically significant. Analyses were performed using GraphPad Prism software (version 9.0).

To identify key genes involved in the pathogenesis and progression of breast cancer, we performed differential expression analysis between tumor and normal samples using datasets from TCGA and GSE42568 databases. As shown in Figures 1A, B, thousands of differentially expressed genes (DEGs) were identified in both datasets. To refine these candidate genes, we subsequently employed Weighted Gene Co-expression Network Analysis (WGCNA) to identify prognosis-related key module genes (Figures 1C–F). The results demonstrated that the turquoise module showed the strongest tumor correlation in the TCGA dataset, while the blue module exhibited the highest tumor association in the GSE42568 dataset (Figures 1G, H). Finally, through intersection analysis of DEGs from TCGA, GEO datasets, along with genes from the turquoise and blue modules identified by WGCNA, we identified 102 hub genes (Figure 1I). These findings suggest that these 102 hub genes may play crucial roles in breast cancer progression.

To identify an optimal approach for constructing prognostic models, we expanded our analysis by incorporating additional datasets including GSE9893 and GSE7390 alongside previous cohorts. Machine learning analysis was performed on 51 shared genes across these datasets (Supplementary Figures S1B, C). Comparative evaluation revealed that the StepCox method combined with Random Survival Forest (RSF) demonstrated superior performance in prognostic model construction (Figure 2A). Subsequently, RSF modeling was applied to these 51 genes for prognostic prediction. Through VIMP ranking and minimal depth analysis, we identified a 3-gene prognostic signature (Figure 2B), with corresponding error rate validation (Figure 2C) and VIMP visualization (Figure 2D).

The prognostic model based on these three genes exhibited significant survival discrimination (Figures 2E, F). Patients stratified into the high-risk group showed markedly worse clinical outcomes (Figure 2G). Time-dependent ROC analysis revealed robust predictive accuracy, with AUC values of 0.890, 0.920, and 0.936 for 1-, 3-, and 5-year survival predictions, respectively (Figure 2H). These results demonstrate that the machine learning-optimized model integrating StepCox and RSF methodologies provides enhanced prognostic stratification capacity for breast cancer patients.

The validation cohorts (GSE7390, GSE9893, GSE42568, and METABRIC) were utilized to assess the prognostic performance of the model. Results demonstrated that breast cancer patients with higher prognostic signature scores exhibited significantly worse clinical outcomes (Figures 3A–D). Time-dependent ROC analysis revealed area under the curve (AUC) values of 0.6298 for GSE7390, 0.7023 for GSE9893, 0.7141 for GSE42568, and 0.6028 for METABRIC, further confirming the model’s predictive reliability across independent datasets (Figures 3A–D). Further stratification based on expression levels of the three model-constructing genes (EGR3, RECQL4, and MMP1) revealed two distinct breast cancer subtypes (Figure 3E), with significant inter-subtype expression differences visualized in Figure 3F. The integrated model combining the three-gene signature with key clinical variables (ER, PR, HER2, and Ki67 status) demonstrated a superior prognostic ability compared to the model with clinical variables alone, as evidenced by higher AUC values across 5-year survival predictions (Figure 3G). ESTIMATE algorithm-based analysis demonstrated that patients in the C1 subtype exhibited significantly higher tumor purity (Figure 3H), along with reduced ESTIMATE scores (Figure 3I), lower immune infiltration (Figure 3J), and diminished stromal component (Figure 3K) compared to other subgroups.

To investigate the mechanistic basis of the three-gene prognostic signature (EGR3, RECQL4, MMP1) in breast cancer progression, we performed differential expression analysis followed by GO and KEGG pathway enrichment analyses between the two molecular subtypes. Functional annotation revealed significant enrichment of immune-related pathways, including myeloid cell migration, neutrophil chemotaxis, and IL-17 signaling pathway among differentially expressed genes (Figures 4A–C).

To validate the immunological relevance of the prognostic model, we systematically profiled immune-related molecules across subtypes. The C2 subtype exhibited elevated expression of immune checkpoint inhibitors, including CD274 (PD-L1), LAG3, and TIGIT (Figure 4D), along with upregulated HLA family genes critical for antigen presentation (Figure 4E). Notably, CIBERSORT deconvolution of TCGA breast cancer samples demonstrated an immunosuppressive microenvironment in the C2 subgroup, characterized by reduced CD8⁺ T cell infiltration, enriched M2 macrophage polarization, and expanded neutrophil populations (Figure 4F). These findings collectively establish that our prognostic signature is mechanistically linked to tumor immunomodulation.

To elucidate the functional roles of the three prognostic signature genes (EGR3, RECQL4, and MMP1) in breast cancer pathogenesis, we conducted comprehensive multi-omics analyses across multiple datasets. Pan-cancer expression profiling revealed significant downregulation of EGR3 in tumor tissues compared to normal counterparts, with progressive reduction across advancing tumor stages. Survival analysis demonstrated favorable clinical outcomes in EGR3-high patients, contrasting with the oncogenic patterns of RECQL4 and MMP1, though their survival correlations lacked statistical significance (Figures 5A–C). Immune infiltration analysis revealed EGR3 exhibited inverse correlation with tumor purity, positive association with CD8⁺ T cell infiltration, and negative correlation with myeloid-derived suppressor cells (MDSCs). Conversely, RECQL4 and MMP1 displayed opposing immunomodulatory patterns (Figures 5D–F). Single-cell RNA sequencing across tumor microenvironment compartments demonstrated ubiquitous but cell type-specific expression patterns, with notable enrichment in malignant cells observed in select datasets (Supplementary Figures S2A–C; Figure 6A). Subtype-specific evaluation revealed lowest EGR3 expression in triple-negative breast cancer (TNBC), correlating strongly with immune response signatures. Multi-cohort validation confirmed EGR3-high tumors exhibited superior prognosis (Figures 6B–D; Supplementary Figure S3A). In contrast, RECQL4 and MMP1 showed minimal subtype-specific variation and immune association, though differential survival patterns emerged (Figures 6E–J; Supplementary Figures S3B, C). Cell line validation confirmed TNBC-specific EGR3 suppression (Supplementary Figures S3D–F).

Prioritizing EGR3 for mechanistic validation, clinical specimen analysis confirmed significant EGR3 downregulation in breast carcinomas (Figure 7A). TIDE database interrogation revealed positive correlation between EGR3 expression and cytotoxic T lymphocyte (CTL) infiltration (Figure 7B), with EGR3-high cohorts demonstrating enhanced immunotherapeutic response (Figure 7C). Furthermore, EGR3 expression was significantly correlated with the levels of multiple immune-related molecules, including LAG3, CD274 (PD-L1), TIGIT, and CD69 (Supplementary Figure S4). Functional studies in 4T1 murine models demonstrated EGR3 overexpression significantly inhibited tumor proliferation and suppressed migratory capacity (Figures 7D–F). In vivo xenograft experiments revealed 56.7% tumor volume reduction with EGR3 overexpression (Figures 7G, H). Flow cytometric analysis of tumor-infiltrating lymphocytes demonstrated enhanced CD8⁺ T cell activation and increased TNF-α and IFN-γ production in EGR3-overexpressing tumors (Figures 7I–L). These findings suggest an association between EGR3 expression and tumor-suppressive immunomodulatory effects with therapeutic potential in breast cancer. To determine if EGR3 expression directly augments T cell effector function, we performed a tumor cell-PBMC co-culture assay. We found that EGR3-overexpressing tumor cells significantly boosted the proportion of cytokine-producing CD8⁺ T cells, indicating a potent role for EGR3 in promoting T cell activation (Figures 7M–O).