We aggregated data from 11 independent breast cancer cohorts derived from four databases: The Cancer Genome Atlas (TCGA), the Gene Expression Omnibus (GEO), and Metabric (8). For a robust and comprehensive prognostic analysis, we focused on samples with complete survival data. Detailed dataset characteristics, including platform usage and cohort breakdowns, are summarized in Supplementary Table S1 . Epigenetic regulators were sourced from the EpiFactors database (9).

Utilizing the approach of Liu et al. (10), we incorporated 10 computational tools to construct a distinctive telomerase signature for breast cancer: Random Survival Forest (RSF), Least Absolute Shrinkage and Selection Operator (LASSO), Gradient Boosting Machine (GBM), Survival Support Vector Machine (Survival-SVM), Supervised Principal Component (SuperPC), Ridge Regression, Partial Least Squares Cox Regression (plsRcox), CoxBoost, Stepwise Cox Regression, and Elastic Net (Enet). RSF, LASSO, CoxBoost, and Stepwise Cox played crucial roles in dimensionality reduction and variable selection, resulting in 108 unique configurations to create a predictive signature. We assessed all cohorts, including TCGA and eight external datasets, using the average Concordance index (C-index) to determine the most reliable prognostic model, thereby establishing a specific signature for predicting breast cancer outcomes.

To further refine our model and ensure it included only the most predictive genes, we employed an exhaustive search strategy. After identifying 28 candidate genes using RSF and univariate Cox proportional hazards regression, we evaluated all possible combinations of these genes to identify the subset that provided the best model performance. Performance was assessed using predefined criteria, including the C-index and Akaike Information Criterion (AIC). This process reduced the candidate genes to nine, which demonstrated the highest prognostic value across all training datasets.

A risk score for each patient was then calculated using the expression levels of the selected genes weighted by their regression coefficients from the Cox proportional hazards model. The final nine-gene signature was validated across multiple independent cohorts, demonstrating consistent predictive accuracy for breast cancer outcomes.

We analyzed genetic variations between two MLEM groups using the TCGA-BRCA database, focusing on both mutation levels and Copy Number Alterations (CNA). Tumor Mutation Burden (TMB) for high- and low-MLEM breast cancer patients was derived from raw mutation data, with the most frequently mutated genes (mutation rate > 5%) visualized via maftools. Patient-specific mutational signatures were identified using the deconstructSigs tool (11), revealing four significant mutational signatures (SBS2, SBS13, SBS7b, SBS7d) with elevated mutation frequencies. Additionally, we identified the five most common regions of amplification and deletion, particularly in genes located at 17q23.1 and 15q13.1.

To process single-cell RNA sequencing data, we utilized the GSE161529 dataset with Seurat (v4.0) (12). Initially, we excluded genes with zero expression, focusing only on those with detectable expression levels. The expression matrix was normalized using Seurat’s “SCTransform” function, followed by dimensionality reduction via PCA and UMAP techniques. Cellular groupings were identified using the “FindNeighbors” and “FindClusters” functions. The DoubletFinder package was applied to remove potential doublets, ensuring data integrity (13). Cells with over 15% mitochondrial genes or fewer than 500 genes were excluded. Ultimately, 50,214 cells passed quality control and were categorized by manually annotating cell types based on established marker genes.

In our study, we adapted the Single-Cell rEgulatory Network Inference (SCENIC) methodology to construct gene regulatory networks (GRNs) from single-cell RNA sequencing data (14). SCENIC functions through a three-step process: identifying co-expression modules between transcription factors (TFs) and potential targets, pinpointing direct targets using enriched TF motifs, and defining regulons composed of a TF and its direct targets. We calculated the regulatory activity score (RAS) for each cell using the area under the recovery curve. To address SCENIC’s limitations with large datasets and sequencing depth variability, we preprocessed data into metacells, enhancing scalability and robustness, thereby significantly improving data handling and computational efficiency (15).

Our study employs an advanced computational approach to map the regulatory crosstalk among TFs and their target genes, with a focus on TF clustering. The process starts by filtering TF-target interaction data to concentrate on significant pairs (significance threshold > 1), emphasizing the most relevant regulatory interactions. Key regulatory TFs, termed hub genes, are identified by quantifying their regulation of target genes. These interactions are represented using an undirected graph model, spatially refined by a force-directed algorithm to clearly illustrate the network architecture and TF-target interactions. Additional structural insights are obtained through the Leiden algorithm, which detects community structures and groups TFs into clusters based on their regulatory links, enhancing our understanding of the regulatory landscape.

We employed the “CellChat” R package to analyze cell-cell communication (16), creating CellChat objects from UMI count matrices and utilizing the “CellChatDB.human” database for ligand-receptor interactions. The analysis was conducted with default settings, merging objects via the “mergeCellChat” function. Differences in interaction number and intensity across cell types were visualized using “netVisual_diffInteraction.” Changes in signaling pathways were assessed with “rankNet,” while gene expression distributions were illustrated using “netVisual_bubble” and “netVisual_aggregate.” Additionally, the NicheNet package was used to explore ligand activity and the regulated expression of downstream targets (17), providing deeper insights into signaling dynamics and communication pathways within the cellular microenvironment.

To evaluate immune cell infiltration in the tumor microenvironment (TME), we employed multiple algorithms: MCPcounter, EPIC, xCell, CIBERSORT, quanTIseq, and TIMER (18–23). These analyses enabled patient categorization by their MLEM scores, offering a detailed view of the immune landscape. To benchmark their consistency, we calculated Spearman’s correlation coefficients between the outputs of these methods for major immune cell populations. We also assessed the ESTIMATE and TIDE indices to understand immunotherapy potential and prognostic implications for breast cancer (24, 25). Quantification of immune checkpoints was conducted to predict patient responsiveness to immune checkpoint inhibitor (ICI) therapy, thereby supporting personalized medicine and optimizing treatment strategies.

To identify potential therapies for high-MLEM patients, we first filtered out duplicate compounds from the Drug Repurposing Hub, resulting in 6,125 unique compounds. We performed Spearman correlation analysis to select genes associated with breast cancer outcomes, targeting those with correlation coefficients above 0.15 (P < 0.05) and those indicating poor prognosis with coefficients below -0.15 (P < 0.05). Gene significance was further evaluated using CERES scores from the Cancer Cell Line Encyclopedia (CCLE) in relation to breast cancer cell (26).

We also assessed drug responsiveness using data from the Cancer Therapeutics Response Portal (CTRP) and the PRISM project, which involve drug screening across various cancer cell lines. The predictive accuracy of drug responses was enhanced using the pRRophetic package’s ridge regression model, validated through 10-fold cross-validation (27).

Furthermore, we explored potential drugs using Connectivity Map (CMap) analysis by comparing gene expression profiles and identifying compounds inversely related to CMap scores, indicating higher therapeutic potential against breast cancer.

For gene expression analysis in breast cancer samples, RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and cDNA synthesis was performed using GoScript reverse transcriptase and Master Mix (Promega) according to the manufacturer’s instructions. Quantitative expression was measured via qRT-PCR on the CFX96 Touch Real-Time PCR Detection System (BioRad, Hercules, CA, USA) utilizing the 2^-ΔΔCq method with GAPDH as a normalization control. Patients were then stratified based on gene expression levels calculated using a formula derived from the MLEM, aiding in the identification of varying risk profiles and the development of tailored treatment strategies.

Breast cancer tissue samples were collected from 30 patients at Guizhou Provincial People’s Hospital, and Hematoxylin and Eosin (HE) staining was performed following standard protocols, with diagnoses confirmed by two independent pathologists. For immunohistochemistry (IHC) analysis on paraffin-embedded samples, procedures and scoring systems from our previous studies were followed (28, 29). Protein expression levels were independently assessed by the same pathologists, ensuring consistency with our earlier research (29).

To assess the correlation between epigenetic regulation and BC prognosis, we collected epigenetic regulators from the EpiFactors database and applied them across nine independent cohorts. This multi-platform dataset enabled the construction of a robust model using machine learning techniques. The model’s performance was evaluated by calculating the average C-index across the nine cohorts using 10 machine learning algorithms in 108 combinations ( Figure 1A ). The RSF model, which achieved the highest score, was selected as the final model. Genes corresponding to the point with the lowest error rate were chosen as candidate variables ( Figures 1B, C ). Subsequently, a univariate Cox analysis was performed on these genes to assess their impact on BC prognosis ( Figure 1D ). Ultimately, nine genes were identified to construct a machine learning-based epigenetic model (MLEM) ( Figure 1E ). The MLEM effectively stratified BC patients into different risk groups with accurate prognostic predictions ( Supplementary Figure S1 ).

Univariate and multivariate Cox analyses were utilized to compare MLEM with other clinical indicators, revealing that MLEM possessed strong predictive capability ( Supplementary Figure S2A ). A nomogram incorporating MLEM, staging, and age was developed to predict the overall survival (OS) of patients at 1, 3, and 5 years ( Supplementary Figure S2B ). The analysis of calibration curves and Hosmer-Lemeshow tests showed a high level of consistency with the standard curve ( Supplementary Figures S2C, D ). Additionally, decision curve analysis (DCA) suggested that MLEM holds significant potential for clinical application ( Supplementary Figure S2E ). Compared to other clinical factors, MLEM demonstrated superior discrimination ( Supplementary Figure S2F ).

Subsequently, we manually curated 68 published BC models for evaluation against MLEM. Univariate Cox analysis was conducted on these models across 10 BC cohorts. MLEM exhibited the best stability across all 10 datasets ( Figure 2A ). Furthermore, the predictive ability was evaluated using the average C-index, where MLEM consistently ranked high in all 10 cohorts, underscoring its robustness and effectiveness in BC prognostication ( Figure 2B ).

A multi-omics analysis was conducted on genes potentially associated with the BC process in patients with varying MLEM. The findings revealed that patients with high-MLEM group exhibited a significantly higher tumor mutation burden (TMB) compared to those with low-MLEM, which was accompanied by a higher rate of multigene mutations ( Figures 3A, C ). When integrating mutation data from ten oncogenic signaling pathways, it was observed that the mutation rates of classic tumor suppressor genes such as TP53, STK11, and AXIN1/2 were significantly elevated in the high-MLEM group. Conversely, the mutation rates of classic oncogenic genes such as PIK3CA/B and AKT1/2/3 were lower in the high-MLEM group ( Figures 3A, B ). Additionally, a comparison of copy number variations between the two MLEM groups indicated more pronounced amplifications and deletions of chromosome arms in the high-MLEM group. Notable amplifications included 4q13.3, 8q24.21, 17q12, 17q23.1, and 20q13.2, while significant deletions were observed at 9p21.3, 9p23, 15q13.1, 18q23, and 19p13.3 ( Figures 3A, D ). High-MLEM patients exhibited a greater prevalence of SBS13, which is associated with increased tumor mutation burden and aggressive tumor phenotypes. This finding aligns with the poorer prognosis observed in this group. SBS2 and SBS13, indicative of APOBEC activity, have been shown to correlate with heightened sensitivity to immune checkpoint inhibitors. This may explain the differential therapy responses between MLEM groups. Additionally, the presence of SBS7b and SBS7d suggests potential vulnerabilities to therapies targeting nucleotide excision repair pathways.

Single-cell transcriptome analysis was conducted to evaluate the biological signifi-cance of MLEM in 14 BC patients, comprising 6 normal tissues and 8 BC tissues ( Supplementary Figures S3A, B ). Initially, 18 clusters were identified, and 8 cell types were further determined ( Figures 4A, B ). The proportions of these 8 cell types were then analyzed, highlighting the differences between the patient’s body and the 14 tissue samples ( Supplementary Figures S3C, D ). Veri-fication of these cell types using representative markers confirmed consistency with the original results ( Figure 4C ; Supplementary Figure S3E ). Analysis of the distribution of these cell types revealed that epithelial cells, plasma cells, macrophages, T cells, and B cells were more prevalent in BC tissue, whereas pericytes and endothelial cells were more abundant in normal tissue. Fibroblasts were found in similar proportions in both BC and normal tissues ( Figure 4D ).

The previously established MLEM was integrated into the single-cell analysis and categorized into high- and low-MLEM groups based on the peak of epithelial cells ( Figures 4E, F ). Differential gene expression between the MLEM groups was then analyzed across these 8 cell types ( Supplementary Figure S3F ). For instance, in epithelial tumor cells, genes related to rRNA binding, translation regulatory factors, and nucleic acid binding were significantly upregulated in high-MLEM cells, whereas genes associated with histone binding were significantly upregulated in low-MLEM cells ( Supplementary Figure S3G ). Additionally, the copy number alteration of genes in epithelial cells was examined using the CNA-based copyKat algorithm ( Figure 4G ), revealing that aneuploid tumor cells were more prevalent than diploid normal cells ( Figure 4H ).

To comprehensively construct GRNs, we utilized the SCENIC pipeline to analyze single-cell RNA-seq data along with cis-regulatory sequence information. This process transformed the gene expression data into RAS for TFs ( Figures 5A, B ). We then performed variance decomposition analysis using principal component analysis (PCA) to identify specific regulons for MLEM and different cell types. The results indicated that PC1 accounted for cell type-specific TFs, while PC2 was associated with MLEM-specific TFs ( Figures 5C, D ).

In our investigation to identify key regulators of cell identity, we evaluated the activity of each regulon across various cell types. We assigned a regulon specificity score (RSS) based on Jensen-Shannon divergence to measure each regulon’s association with specific cell identities. By concentrating on regulons with the highest RSS values, we analyzed their functional characteristics ( Figure 5E ). This analysis identified MESP1, SPDEF, and GATA3 as key regulons uniquely associated with epithelial cells. UMAP plots further highlighted the specificity of these regulons’ activities within epithelial cells ( Figure 5F ). Additionally, we presented and analyzed regulons specific to other cell types, underscoring the distinct regulatory networks that define cell identities ( Supplementary Figure S4A ).

Recognizing that TFs collaborate to regulate gene expression, we systematically examined the combinatory patterns of these regulatory elements. By assessing the similarity of RAS scores for each regulon pair across the entire atlas using the Leiden algorithm, we organized 343 regulons into 12 distinct modules, unveiling complex regulatory patterns within the cellular landscape ( Figure 5G ; Supplementary Figure S4B ). Notably, modules A and C were predominantly associated with MLEM ( Figure 5H ). We then investigated the specific TFs driving transcriptomic changes in epithelial cells influenced by MLEM. GSEA identified several pathway alterations ( Figure 5I ), including the translation initiation in epithelial cells with high-MLEM ( Figures 5I, J ). Further analysis identified the TFs contributing to this pathway and their roles in MLEM progression ( Figure 5K ). A comprehensive regulatory network is illustrated in Figure 5L .

The relationships of intercellular communication differences between MLEM groups among 8 types of cells were analyzed using CellChat. Although the type of cellular communication was more diverse, the intensity of interactions was slightly lower in high-MLEM cells ( Figure 6A ). Differences in intercellular communication revealed a significant decrease in the number of interactions involving plasma cells with other cells, whereas interactions among the other cells significantly increased in high-MLEM cells ( Figure 6B ). Additionally, the analysis of intercellular signaling pathways showed that seven pathways were highly active in low-MLEM cells, while 18 pathways were highly active in high-MLEM cells ( Figure 6C ).

Further analysis assessed the strengths of outgoing and incoming interactions to gauge cell-cell communication. Compared to low-MLEM cells, macrophages and endothelial cells showed a significant increase in their signal reception ability, whereas plasma cells and T cells exhibited a significant decrease in high-MLEM cells ( Figure 6D ). Several pathways in epithelial cells were specific to MLEM, such as LAMININ ( Figure 6E ).

A deeper analysis was performed to summarize the expression of receptor ligands in the most valuable signaling pathways ( Figure 6F ). It was found that MDK-SDC4 and CCN1-SDC4 ligand-receptor pairs were highly expressed in fibroblasts ( Figure 6G ). Finally, signaling communication pathways involving MDK-SDC4 and CCN1-SDC4 ligand-receptors were analyzed, identifying potential pathways ( Figure 6H ).

Six algorithms were utilized to evaluate the immune microenvironment, revealing that patients in the low-MLEM group exhibited a higher degree of immune cell infiltration, including CD4⁺ and CD8⁺ T cells ( Figure 7A ). The comparison of immune checkpoint inhibitors (ICIs) indicated higher expression levels of ICIs in the low-MLEM group, further suggesting that immunotherapy might be more suitable for these patients ( Figure 7B ). Representative cell markers were evaluated using IHC experiments ( Figure 7C ). Spearman’s correlation analysis revealed a high degree of concordance between MCPcounter, quanTIseq, and TIMER for key immune populations such as CD8⁺ T cells (ρ > 0.75). In contrast, xCell exhibited lower concordance with other methods, particularly for macrophages and B cells (ρ < 0.50). Despite these differences, all methods consistently showed higher immune infiltration levels in the low-MLEM group compared to the high-MLEM group. This robustness across algorithms supports the reliability of our findings.

Furthermore, patients in the low-MLEM group had superior ESTIMATE scores, immune scores, and stromal scores, whereas high-MLEM patients had higher tumor purity, suggesting that the low-MLEM group might be more responsive to immunotherapy ( Figure 8A ). Subsequently, TIDE analysis showed that TIDE, Dysfunction, and Exclusion indicators were higher in the low-MLEM group compared to the high-MLEM group ( Figure 8B ), and patients in the low-MLEM and low-TIDE group demonstrated significantly better outcomes ( Figure 8C ). Further analysis revealed a significant correlation between immune signaling pathways, the immune cycle, and MLEM in BC ( Figure 8D ). Finally, studies from IMvigor210 ( Figures 8E–H ) and GSE78220 ( Figures 8I–L ) indicated that patients in the low-MLEM group benefited more from PD-L1 and PD-1 administration, confirming that immunotherapy was more suitable for patients in the low-MLEM group.

Spearman correlation analysis was applied to investigate the correlation between MLEM and gene expression (positive correlation) as well as CERES (negative correlation). This analysis identified seven potential therapeutic targets for patients with high MLEM: CYCS, SLC25AB, COX7B, NDUFB3, NDUFA4, NDUFB6, and NDUFB9 ( Figure 9A ). These targets were found to be closely related to various pathways of action of breast cancer-related drugs, marking them as key therapeutic targets for BC patients ( Figure 9B ).

Subsequently, six potential compounds were identified from the CTPR (SB-743921, paclitaxel) and PRISM datasets (deforlimus, romidepsin, ispinesib, vincristine). The AUC values for these compounds were compared between the two MLEM groups, revealing lower AUC values in high-MLEM patients, which suggests a better chemotherapy efficacy for these patients ( Figures 9C, D ).

A multi-perspective analysis was performed to identify the most effective therapeutic drugs among the six candidates. This analysis incorporated clinical status, experimental evidence, mRNA expression levels, and CMap scores for each compound. Vincristine emerged as the potential therapeutic drug for high-MLEM patients, based on its CMap score ( Figure 9E ).