The mRNA raw count profiles of the TCGA-BRCA project and relevant clinical information were downloaded from the GDC data portal (https://portal.gdc.cancer.gov). A total of 1097 breast cancer patient samples with corresponding clinical information were available in TCGA. In addition, 888 DNA methylation profiles (Methylation450k, 790 tumor samples, and 98 non-tumor samples) of breast cancer patients were acquired from the UCSC Xena website (http://xena.ucsc.edu). Among 790 TCGA breast cancer DNA methylation samples, 785 samples included both RNA-sequencing data and paired DNA methylation data. In this paper, the gene methylation value was defined as the average DNA methylation value of all CpG sites in the promoters (transcription start sites (TSS) 1500 and TSS200) of a gene. An independent cohort (GSE86166 microarray dataset) with gene expression profiles of 366 breast cancer samples was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo).

The methylation level of each probe was represented by the β value ranging from 0 (unmethylated) to 1 (fully-methylated). To simplify the DNA methylation data, probes with missing data in more than 50% of the samples were removed. The remaining probes with NAs were imputed by the K-nearest neighbor (KNN) imputation procedure. For raw count data from TCGA, the Ensembl IDs were transformed to gene symbols, and protein‐coding genes were selected. Next, we computed the transcripts per million (TPM) values, which showed a greater similarity to those generated from microarray analysis and higher comparability among samples. Then, TPM values were normalized using the log2(TPM+1) formula. For the GEO datasets, we normalized the expression datasets by Robust Multiarray Average with the R package “affy” (16). For multiple probes mapped to one gene, the mean value of expression was taken.

In this study, the methylation analysis following chip analysis methylation pipeline (ChAMP) was performed with the required R package “IlluminaHumanMethylation450-kanno.ilmn12.hg19” as annotation (17, 18). Methylation β matrix was first filtered and imputed and then normalized through the embedded BMIQ method with the champ.filter function was applied to filter methylation probes with default parameters. The β value was normalized by the champ.norm function. Differentially methylated sites (DMSs) were detected by champ.DMP function. A probe was identified as a hypermethylated probe if the probe methylation level was greater than 0.3 in the tumor group but less than 0.2 in the normal group with an adjusted P value less than 0.05, and vice versa for hypomethylated probes. Furthermore, differentially methylated genes (DMGs) were determined according to probe location in corresponding genes. The average DNA methylation value for all CpG sites in the promoters (transcription start sites (TSS) 1500 and TSS200) of a gene was defined as the DNA methylation value for that gene.

The MethylMix package in R was used to analyze integrated DNA methylation data of 785 breast cancer samples, 98 non-tumor samples, and paired gene expression data to identify DNA methylation events that significantly affect the expression of corresponding genes and to show that the gene was a DNA methylation-driven gene (19). For the MethylMix analysis, the correlation between the methylation data and paired gene expression data of DMGs in 785 breast cancer samples was determined to identify DNA methylation events that led to changes in gene expression, and only genes that met the correlation filter were recruited into further analysis. We also used a β mixture model to define a methylation state across a large number of patients, precluding the need for an arbitrary threshold. DNA methylation states between the 785 breast cancer samples and 98 corresponding non-tumor samples and was compared by conducting a Wilcoxon rank sum test. In this paper, significant q values of 0.05 calculated using P-value multiple testing correction with false discovery rate (FDR), and other parameters were set as default parameters.

By using the random forest as a classifier, further feature selection was performed on DNA methylation-driven genes (20). Multivariate cox regression analyses were further conducted to evaluate relationships between the expression of the DNA methylation-driven genes and breast cancer prognosis and to detect the independent DNA methylation-driven genes significantly associated with the cancer prognosis. A risk score prediction model based on DNA methylation-driven genes was established through a linear combination of the expression levels of independent DNA methylation-driven genes using coefficients from multivariate Cox regression as weights. Breast cancer patients were stratified into low-risk and high-risk groups by the prediction model, with the optimal risk score as the cutoff point. Survival differences between high-risk patients and low-risk patients were evaluated by Kaplan-Meier survival plots and then compared by log-rank test. The GSE86166 dataset from the GEO database was used to validate the prognostic model. Whether the predictive power of the predictive model was independent of other clinical features of breast cancer patients was analyzed by univariate and multivariate Cox regression analyses.

In this paper, we calculated the tumor immune signatures from the following two aspects. Firstly, the levels of infiltrating immune and stromal cells were calculated by CIBERSORT (21), TIMER (22), MCP-counter (23), and quanTIseq (24) algorithms. Secondly, the expression of the leukocyte antigen (HLA) family genes and immune checkpoints genes also was calculated.

The information of somatic mutations in TCGA samples was downloaded from Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Concerning different mutation types, non-synonymous mutation variants and synonymous mutation variants are analyzed respectively in this paper. Differential mutated genes between the low- and high-risk groups were identified and the interaction effects of mutated genes also were analyzed by maftools package (25). In our analysis, only genes mutating more than 10 times in at least one group were considered.

Expression profile data of human cancer cell lines were obtained from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) (26). The CERES scores of genome-scale CRISPR knockout screening of 18,333 genes in 739 cell lines were acquired from the dependency map (DepMap) portal (https://depmap.org/portal/) (27, 28). CERES score measures the dependency of the gene of interest in a specific cancer cell line, with a lower score indicating a higher possibility of the essentiality of the gene in cell growth and survival of the given cancer cell line. Drug sensitivity data of cancer cell lines were obtained from the Cancer Therapeutics Response Portal (https://portals.broadinstitute.org/ctrp) and PRISM Repurposing dataset (https://depmap.org/portal/prism/) (29, 30). The CTRP contains the sensitivity data of 481 compounds from over 835 cancer cell lines, and the PRISM contains the sensitivity data of 1497 compounds over from 480 cancer cell lines. Both two datasets provide area under the dose-response curve (AUC, area under the curve) values as a measurement of drug sensitivity, with a lower AUC indicating a higher sensitivity to treatment. In this study, compounds with more than 20% of missing data and cell lines derived from hematopoietic and lymphoid tissue were excluded. The KNN imputation was applied to impute the missing AUC values. Molecular data in CCLE were used for subsequent CTRP and PRISM analyses, as the cancer cell lines in both datasets were obtained from the CCLE project.

Human breast cancer cell line MDA-MB-231 and non-malignant breast epithelial cell line MCF10A were kind gifts from the Cancer Hospital of the University of Chinese Academy of Sciences (Zhejiang Cancer Hospital, Hangzhou China). The human breast cancer cell line MCF7 was purchased from American Type Culture Collection (ATCC). All the cell lines were maintained in RPIM 1640 medium (Gibco, CA, USA) supplemented with 10% fetal bovine serum (FBS, BI, Israel), 100 unit/mL penicillin, and 100 mg/mL streptomycin in a 5% CO₂, 37°C incubator. Oligomycin A (Lot#: S147806), panobinostat (Lot#:120036), (+)-JQ1 (Lot#:S711013), voxtalisib (Lot#:120367), and arcyriaflavin A (Lot#:E0817) were commercially purchased from Selleck (USA).

Cell counting kit-8 assay (CCK-8, Beyotime, China) was performed to detect the sensitivity of the breast cancer cell lines to the drugs (31). In short, when the cells were in good condition and growing exponentially, the culture medium was discarded and the cells were washed with PBS and incubated with trypsin for 3 minutes. The reaction was then terminated with the serum-containing culture medium to prepare cell suspension. Cells (1×10⁴ cells/well) were inoculated into 96-well plates, cultivated overnight to adhere to the wall. All five compounds were dissolved in DMSO at various concentrations and added to the 96-well plates for 24 h. CCK-8 assay was performed to determine cell viability according to the instructions and to assess the sensitivity of breast cells to the drugs. The IC₅₀ values were calculated by GraphPad Prism 8.

The mRNA expression levels of methylated-gene CENPE in the MCF10A, MDA-MB-231, and MCF7 cell lines were examined. The primer sequence used for PCR is F: 5’-GATTCTGCCATACAAGGCTACAA -3’, R: 5’-TGCCCTGGGTATAACTCCCAA -3’. Briefly, the total RNA samples were obtained from all three cell lines using RNAiso Plus (lot# 9108, Takara, Japan) according to the manufacturer’s instruction. 5× Prime Script TM RT Master Mix (Perfect Real Time, lot# AJ21979A, Takara, Japan) was used to reverse-transcribe 20 uL RNA into the cDNA. The relative gene expression was quantified by the Universal SYBR Green Master (lot# 50837000, Roche, Switzerland). The Applied Biosystems (USA) was used for PCR amplification. The data were calculated with the 2^-△△Ct.

When cells increased exponentially and had normal cell morphology, PBS was used washed and trypsin digestion. Next, cell suspensions were seed in 12-well plates. Cell state was observed and incubated with various concentrations of drug solution. Cultivated 24h, single-cell suspensions were prepared and incubated with APC Annexin V apoptosis Detection Kit with 7-AAD (lot#: 640930, BioLegnd, USA) for 15min, followed by flow cytometry analysis (32).

R software was used for statistical analysis of all the data. Survival analysis was performed using the R package survival, and Kaplan-Meier plots and log-rank tests were conducted to assess the difference in overall survival (OS) between the two groups (33). GSEA (gene set enrichment analysis) was performed using the package clusterProfiler in R (34). The Chi-square test was used to verify the association between categorical clinical information and defined groups. The Wilcoxon test (Mann–Whitney test) was used for continuous data. For all statistical analyses, a P value less than 0.05 was considered statistically significant.

The workflow of this study is shown in Figure S1 . Under the previously mentioned thresholds, we detected a total of 2415 upregulated CpG sites and 134 downregulated sites. The detailed differentially methylated CpG sites are shown in Supplementary Table S1 . The circle plot of differentially methylated sites is shown in Figure 1A and the top 100 sites heatmap is shown in Figure 1B . These differentially methylated sites were mapped onto corresponding genes. A gene was defined as differentially methylated if there is a differentially methylated CpG site on its promoter. After integrating the mRNA expression of genes, we finally obtained 1135 differentially methylated genes (DMGs, Supplementary Table S2 ).

As DNA methylation-driven genes play an important role in the initiation and progression of breast cancer, we further identified the DNA methylation-driven genes in this disease. Firstly, the DNA methylation and gene expression data of DMGs (DNA methylation data of 785 breast cancer samples and 98 non-tumor samples and the paired gene expression data of 785 breast cancer samples) were selected for further analysis. Then, using the MethylMix package, 194 DNA methylation-driven genes in total were screened ( Supplementary Table S3 ). To examine the relationship between methylation-driven genes and the prognosis of breast cancer, the univariate cox proportional hazard regression analysis was performed and identified 28 statistically significant (P < 0.05) genes and the Kaplan-Meier survival analysis detected 33 statistically significant (P < 0.05) genes ( Supplementary Tables S4-5 ). Finally, 17 intersect genes were obtained for further analysis ( Figure 1C ). The feature selection was conducted by recursive feature elimination (RFE) with a random forest as the classifier and a 10-fold cross-validation method in R package caret (35). We obtained an accuracy of 0.854 with 13 DNA methylation-driven genes ( Figure 1D ). These 13 DNA methylation-driven genes were AKR1E2, LIMD2, STAC2, SLC7A4, MAP4K1, PLAT, CXCL1, NACAD, RBP1, SFRP1, TOX, OSR1, and UBA7. The correlation analysis showed that the expression of these DNA methylation-driven genes had a significantly negative correlation with the level of DNA methylation ( Figures S2A–M ), showing a negative regulatory relationship between DNA methylation and gene expression.

By using recursive feature elimination (RFE) with the random forest as a classifier and a 10-fold cross-validation method in R package caret, we obtained 13 DNA methylation-driven genes. A prognostic model was established with the regression coefficient from a multivariate Cox proportional hazard model. The risk score was calculated according to the formula: 0.1944 × AKR1E2 expression level + 0.0335 × SFRP1 expression - 0.0930 × LIMD2 expression - 0.0644 × STAC2 expression - 0.0965 × SLC7A4 expression level - 0.0042 × MAP4K1 expression level - 0.0932 × PLAT expression level - 0.1366 × CXCL1 expression - 0.0340 × NACAD expression - 0.0644 × RBP1 expression -0.1161 × TOX expression -0.0767 × OSR1 expression -0.1658 × UBA7 expression. The risk score was normalized to 0 to 1. The breast cancer patients were stratified into high-risk and low-risk groups according to the optimum cutoff point. The high-risk patients showed markedly worse OS than those with low risk ( Figure 2A ). The predictive performance of the prognostic model was further evaluated using 366 breast cancer samples with OS time and survival status in the validation dataset (GSE86166). The patients in the validation dataset were classified into low-risk and high-risk groups utilizing the earlier mentioned formula based on the optimal cutoff value. Consistent with the above findings, in the validation set patients in the high-risk group showed a markedly shorter median OS than those in the high-risk group ( Figure 2B ). To investigate whether the prognostic model was independent of the clinicopathological features, the univariate and multivariate Cox regression analyses were performed on the TCGA dataset. The results of univariate Cox regression demonstrated that the risk score was significantly associated with OS (high-risk of death group vs. low-risk of death group, HR =2.883, 95% CI 2.087-3.983, P < 0.01, Figure 2C ). Additionally, in multivariable Cox regression, the risk score also showed a significant relationship with OS (high-risk of death group v.s. low-risk of death group, HR = 2.921, 95% CI 2.080-4.102, P < 0.01, Figure 2D ). The same analysis was performed on the validation dataset with similar results ( Figures 2E-F ). These results indicated that the prognostic ability of the prognostic model was independent of other clinical features.

Transcriptome-based gene set enrichment analysis (GSEA) was conducted to determine the pathways and biological processes involved in breast cancer patients with different risks. We found that “Cytokine-cytokine receptor interaction”, “Chemokine signaling pathway”, “PI3K-Akt signaling pathway”, “JAK-STAT signaling pathway”, “Th17 cell differentiation”, and “NF-kappa B signaling pathway” were significantly enriched (P < 0.05) ( Figure 3A and Table 1 ). In addition, the GO enrichment results showed that many biological processes associated with immunity were enriched ( Figure 3B and Table 2 ), which was consistent with the KEGG results. These results indicated that the DNA methylation-driven genes-based risk score could influence these pathways and predict the survival of breast cancer patients.

Proteins presenting a high positive correlation with risk score might have potential therapeutic effects on high-risk patients. However, most of the human proteins lack obvious active sites to which small molecule compounds can bind or have been found in cells where are inaccessible for biological agents. With the purpose of identifying potentially druggable therapeutic targets for high-risk patients, we collected data of 2249 drug targets of 4484 compounds and conducted a two-step analysis for screening candidate targets (36) ( Supplementary Table S6 ). Firstly, the correlation coefficient between the expression of druggable proteins and the risk score was calculated, and 17 protein targets were found to have a correlation coefficient higher than 0.30 (P < 0.05) in the TCGA dataset. Similarly, 29 protein targets were identified from GSE86166 and ten genes, including HMMR, CCNB1, CDK1, CCNA2, CDC25C, P4HA1, PAICS, AURKA, CENPE, and RAD51 were identified in both datasets mentioned above. Finally, the correlation analysis between the CERES score and the risk score based on breast cancer cell lines was conducted. Notably, CDK1, CCNA2, P4HA1, PAICS, and RAD51 showed a CERES score higher than zero in most breast cancer cell lines, indicating that they might not be essential to breast cancer. The remaining five genes HMMR, CCNB1, CDC25C, AURKA, and CENPE were considered to be potential therapeutic targets, and inhibiting the functions of these genes in high-risk patients could potentially achieve a favorable treatment efficacy ( Figure 3C ). The correlation figures of these five genes in TCGA and GSE86166 datasets are shown in Figures 3D-M .

In this research, using CIBERSORT, TIMER, MCP-counter, and quanTIseq algorithms, the abundance of infiltrating immune cells between two groups was calculated and the results demonstrated that most of the immune and stromal cells decreased in the high-risk group ( Figure 4A ). We also calculated the estimate score, immune score, stromal score and tumor purity by using estimate package for each sample. According to the result, we found that estimate scores, immune scores and stromal scores were significantly lower in high-risk patients, while tumor purity was higher ( Figure 4B ). We further investigated gene expression of the 35 immune checkpoints and 19 HLA family genes between the high-and low-risk patients. According to the Wilcoxon test, 32 immune checkpoints and all HLA family genes were significantly modulated in the high-risk group ( Figures 4C-D ). In addition, our analysis also showed that the risk score had a strong negative correlation with the expression of 18 immune checkpoints and 16 HLA genes, including TNFRSF14, CD40LG, CD27, TMIGD2, HLA-DPB1, and HLA-E ( Figures 4E-F ).

In this paper, we downloaded somatic mutations from the TCGA database to investigate risk score-related mechanisms in BRCA. In the mutant frequency level, we found that more somatic mutations were observed in the high-risk group, including non-synonymous and synonymous mutations ( Figure 5A ). The waterfall plots of mutation in all samples were shown in Figure 5B . Meanwhile, differential mutant genes also were calculated and 23 genes mutated more frequently in BRCA patients in the high-risk group, such as TP53, CDH9, RYR3, DYNC2H1, and TMEM132D ( Figure 5C ). However, only 9 genes mutated more frequently in BRCA patients in the low-risk group. Moreover, these differential mutant genes also show significant co-occurring or mutually exclusive mutation patterns ( Figure 5D ), which could shed lights on therapeutic strategies for breast cancer patients.

The CTRP and PRISM datasets, which contain the gene expression profiles and drug sensitivity profiles of hundreds of cancer cell lines can be utilized to construct the prediction model of drug response. We excluded compounds with NAs from more than 20% of the samples and cell lines derived from hematopoietic and lymphoid tissues. 680 cancer cell lines with 354 compounds in the CTRP dataset and 480 cancer cell lines with 1285 compounds in the PRISM dataset were recruited into the subsequent analysis. The pRRophetic package that had a built-in ridge regression model was used to predict the drug response of clinical samples based on their expression profiles, and an estimated AUC value of each compound in each clinical sample was obtained (36).

Two different approaches were adopted to identify the candidate with higher drug sensitivity in high-risk patients. The analyses were performed using CTRP- and PRISM-derived drug response data, respectively. First, differential drug response analyses between high-risk score (top decile) and low-risk score (bottom decile) groups were conducted to identify compounds with lower estimated AUC values in the high-risk score group (log2FC > 0.10). Next, Spearman correlation analysis between AUC value and the high-risk score was used to select compounds with a negative correlation coefficient (Spearman’s r < -0.30 for CTRP or -0.35 for PRISM). The analyses detected two CTRP-derived compounds (oligomycin A and panobinostat) and three PRISM-derived compounds ((+)-JQ1, voxtalisib, and arcyriaflavin A). All these compounds showed lower estimated AUC values in the high-risk group and a negative correlation with risk score ( Figures 6A, B ). Therefore, these compounds may have therapeutic effects on high-risk breast cancer patients.

We identified five potential target genes HMMR, CCNB1, CDC25C, AURKA, and CENPE for patients with high-risk scores. Among them, the CENPE gene encodes a centromere binding protein and mitotic kinesin, which has been demonstrated as a promising target for cancer drug development (37). We further examined the mRNA expression level of CENPE in breast epithelial cell line MCF10A and breast cancer cell lines MDA-MB-231 and MCF7. The results showed that CENPE was significantly elevated in breast cancer cell lines ( Figure 7A ).

To explore the sensitivity of breast cancer cells to the selected drugs, we examined the effects of all five compounds on the viability of MCF10A, MCF7, and MDA-MB-231 cell lines. The results showed that only (+)-JQ1 had selective cytotoxicity for breast cancer cells ( Figures 7B–E and Figure S3 ). (+)-JQ1 showed comparable cytotoxicity in MCF7 and MDA-MB-231 cell lines, with IC₅₀ values of 13.57 μM and 13.72 μM, respectively. However, MCF10A cells (IC₅₀ = 30.74 μM) were much less sensitive to (+)-JQ1 than MCF7 and MDA-MB-231 cell lines. We further examined the effects of (+)-JQ1 on breast cancer cell apoptosis. The results showed that with the increased concentration of (+)-JQ1, the compound induced apoptosis in both breast cancer cell lines ( Figure 7F ). Taken together, our results indicated that (+)-JQ1 might exert its anti-breast cancer activity through CENPE.