RNA-seq data in fragments per kilobase of transcript per million mapped reads (FPKM) form and clinical information of BC were downloaded from the Cancer Genome Atlas (TCGA: https://portal.gdc.cancer.gov/) database. Illumina Human Methylation450 BeadChip array (450k array) and Illumina Human Methylation27 BeadChip array (27k array) data of TCGA database were downloaded from UCSC xena (https://xenabrowser.net/) (Goldman et al., 2020). DNA Methylation levels were evaluated by the β value, which ranged from 0 to 1 (0 means unmethylated and 1 means fully methylated). Probes with over 70% of missing values and probes located at chromosomes X and Y were removed. The missing values of the remaining probes were imputed using the k-nearest neighbours (knn) imputation algorithm of the impute R package. Since DNA methylation in promoter regions would strongly influence gene expression, we focused on the methylation probes in promoter regions defined as 2.0 kb upstream to 0.5 kb downstream from transcription start sites (TSS). Batch effects were removed by the ComBat algorithm of the sva R package (Leek et al., 2012). Ultimately, 560 patients including methylation data (from 450k array) and corresponding clinical data, 986 patients (from TCGA database) containing both gene expression data and corresponding clinical data were used for mainly analysis. And we obtained 557 overlapping patients (from the above two datasets) with complete gene expression data, methylation data and clinical data. Moreover, RNA-seq data and clinical information of 106 samples from GSE146558 were downloaded from NCBI (GEO: https://www.ncbi.nlm.nih.gov/) as external validation dataset. The mRNA expression profile from GEO dataset was normalized by the Robust Multichip Average (RMA) algorithm with background adjustment, quantile normalization, and final summarization. The workflow of our study was illustrated in Figure 1.

Univariate Cox regression analysis was performed to screen out the prognosis-related CpGs of the 560 BC patients. In our study, we used the OS as clinical parameter of prognosis. Next, all these CpGs were subjected to multivariate Cox regression analysis, with age, pathological stage, and TNM stage as covariates to identify independent prognostic CpGs.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/) (Kanehisa et al., 2017) analysis were performed using the R ClusterProfiler package (Yu et al., 2012). A p value < 0.05 was set as the cut-off value for both GO and KEGG analyses in our study.

Consensus clustering (Wilkerson and Hayes, 2010) was performed to determine subgroups of different methylation characteristics of the 560 BC patients based on the independent prognostic CpG loci using the ConsensusClusterPlus package of R. The criteria to determine the number of clusters were: (Freddie et al., 2020) The consistency within the cluster was relatively high; (DeSantis et al., 2019) There was no significant increase in the area under the CDF curve; (Allemani et al., 2018) The relative change in area under CDF curve tended to be stable. We then generated a consensus matrix to better visualize and help determine the number of clusters.

Differentially methylated independent prognostic CpG loci between the different prognosis clusters were screened out using Wilcoxon test. The filtering conditions were false discovery rate (FDR) < 0.05 and | log2 (fold change) | >0.585. On the basis of these differentially methylated CpG loci, a methylation prognostic model of 560 BC patients was constructed using multivariate Cox analysis. The formula of the risk score was as follows: R i s k   s c o r e = ∑ c o e f i β i where coef_i was the multivariate Cox regression coefficient, and βi was the corresponding methylation β value. According to the median risk score (Chen et al., 2020b), (Shen et al., 2019), patients were divided into methylation low-risk (n = 280) and high-risk (n = 280) groups. Survival curves were employed to compare the OS of the two groups. Risk curve was plotted to visualize the relationship of the risk score, survival status and the methylated levels of the six signature CpG loci. Univariate and multivariate Cox regression analysis were performed to explore whether the risk score could be an independent predictor of OS. The sensitivity and specificity of the methylation prognostic model were evaluated by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curve.

Differentially expressed genes (DEGs) were identified from methylation high-risk and low-risk groups (of 557 BC patients’ dataset) with the filter conditions of adjusted p < 0.05 and | log2 (fold change) | >1. And then we further extracted the above-mentioned DEGs from the transcriptomic profiles of 986 BC patients’ dataset for subsequent analysis. Kaplan-Meier analysis and univariate Cox regression analysis were employed to investigate prognosis-related DEGs of the 986 samples. Similarly, a gene prognostic signature was constructed by multivariate Cox regression analysis. And the risk formula was: R i s k   s c o r e = ∑ i = 1 n ( e x p r e s s i o n i ∗ c o e f i )

Patients were also categorized into low-risk (n = 493) and high-risk (n = 493) groups based on the median risk score. Kaplan-Meier survival curve, risk curve, ROC curve, univariate and multivariate Cox regression analysis were also used for evaluating and validating the prognostic signature. Besides, two subgroups from 986 BC patients, the 557 patients’ dataset and GSE146558 dataset were used to validate the prognostic value of the gene signature.

CellMiner (https://discover.nci.nih.gov/cellminer/home.do) is a robust, user-friendly online database that integrates drug sensitivity and genomic data (Reinhold et al., 2017), (Wang et al., 2016). Anti-tumor activity data obtained from NCI-60 tumor cell line panel of the developmental therapeutics program (DTP) and RNA-seq data for the 60 cell lines of the NCI DTP drug screen were downloaded from this website. Subsequently, correlation between the sensitivity of anti-tumor drugs and the signature genes was analyzed.

R 3.6.3 (version 3.6.3, https://pan.baidu.com/s/1sufVf2lmoj9GYG_j5_fJKQ) was used for statistical analysis and plotting. Consensus clustering was performed using the ConsensusClusterPlus package of R; COX regression analysis was performed with the coxph function in survival package of R (Zhang, 2016); Kaplan-Meier curve was plotted using the survival and survminer packages of R; Pheatmap was plotted using the pheatmap package of R; The forest plots were plotted by the forestplot package of R; ROC curve was plotted by the survival ROC package of R. GO and KEGG analyses were performed using the ClusterProfiler package of R.

Mann-Whitney test was used to estimate the statistical significance of two groups of skewed distributed continuous variables, and Kruskal-Wallis test was used to evaluate the statistical significance of multiple groups of skewed distributed continuous variables (with Bonferroni correction for pairwise comparisons among multiple groups). All tests were two-sided and for all statistical tests, p < 0.05 was considered to be statistically significant unless otherwise specified.

In our study, 450k array dataset was defined as train group and 27k array dataset was defined as test group (Table 1). Firstly, 144 CpG loci with p < 0.001 by univariate Cox analysis were screened out and identified as prognosis-related CpG loci. Using age, pathological staging and TNM staging as covariates, 66 CpG loci (49 CpG loci associated with favorable prognosis, 17 CpG loci associated with poor prognosis) with p < 0.001 by multivariate Cox analysis were further selected and used as the methylation classification features (Figure 2).

Then, the 560 patients were categorized into clusters of different methylation characteristics with consensus clustering based on the methylation of the 66 independent prognostic CpG loci. When the patients were assigned to 5 categories, the consistency within the clusters was high, the area under the cumulative distribution function (CDF) curve began to stabilize, and the relative change in area under CDF curve tended to be stable (Figures 3A,B). A consensus matrix representing the consensus for k = 5 also displayed a well-defined 5-block structure (Figure 3C). Accordingly, the optimal number of clusters was determined to be 5.

Subsequently, we conducted a subgroup analysis for the 5 clusters. Firstly, we compared the methylation levels of these 66 independent prognostic methylation loci among the five clusters. As illustrated in Figures 3D,E, Cluster 5 had the lowest methylation levels, followed by Clusters 2 and 4, Clusters 1 and 3. And the methylated difference between Cluster 5 and each of the remaining clusters was statistically significant. To explore the prognostic significance of the five clusters, Kaplan-Meier survival analysis was performed. We found that the prognosis was statistically significantly different among the 5 clusters, where Cluster 1 and Cluster 3 had the best OS, while Cluster 5 had the worst (Figure 3F).

The five clusters were significantly prognosis-associated, and therefore were used to identify potential prognostic biomarkers. Given that Cluster 5 had the lowest methylation level and the worst OS, it was reasonable to be selected as the reference cluster. Next, 20 differentially methylated independent prognostic CpG loci were identified from Cluster 5 and the rest clusters (Table 2). Ultimately, a methylation prognostic model was constructed which included six CpG loci (cg00945507, cg05406101, cg10092957, cg13060154, cg14992108, cg18678121) determined by multivariate Cox analysis (Table 3). Kaplan-Meier analysis showed that cg00945507, cg05406101, cg10092957, cg14992108, cg18678121 were associated with improved survival, and cg13060154 was associated with poor survival (Figure 4).

Then we explored the mechanisms by which these signature CpG loci might act on BC. The six CpG loci, cg00945507, cg05406101, cg10092957, cg14992108, cg18678121, and cg13060154, were located at gene promoter regions of SEC61G, RWDD2B, NCCRP1, SNTB1, SEC61A2, DAB2IP, respectively. We firstly analyzed the correlation of these CpG loci and their corresponding target genes. The methylation of cg10092957, cg05406101, cg18678121, cg00945507 were moderately negatively correlated with the expression of their target genes. Whereas cg13060154 was weakly positively correlated with its corresponding gene (Supplementary Figures S1A–F). Consistent with the above results, the increased expression of NCCRP1, RWDD2B, SEC61A2, and SEC61G was associated with the decreased β values of cg10092957, cg05406101, cg18678121, and cg00945507, respectively. To the opposite of them, DAB2IP had higher expression in the presence of hypermethylated cg13060154 (Supplementary Figures S1G–L). However, there was no relationship between the cg14992108 and its target gene SNTB1. In addition, we took advantage of TCGA Wanderer, an interactive viewer exploring DNA methylation and gene expression data in human cancer (http://maplab.imppc.org/wanderer/) (Díez-Villanueva et al., 2015), to explore the methylated difference of the six CpG loci between breast cancer and normal tissues. We found that the methylation levels of cg05406101, cg18678121, cg14992108 were higher in normal tissues. However, the methylation levels of cg13060154 and cg10092957 were higher in breast cancer (Supplementary Figures S1M–R). We also explored the prognostic roles of the six CpG loci in breast cancer through the public database MethSurv (Modhukur et al., 2018) (https://biit.cs.ut.ee/methsurv/). High methylation levels of cg00945507, cg05406101, cg10092957, cg14992108, cg18678121 were associated with favorable prognosis. On the contrary, high methylation level of cg13060154 was associated with poor survival (Supplementary Figures S1S–X).

On the basis of multivariate Cox regression, we developed the following risk model:

Risk score = −4.016 × cg00945507 + 1.117 × cg05406101−1.78 × cg10092957 + 1.655 × cg13060154−1.283 × cg14992108−1.136 × cg18678121 (Table 3).

According to the formula, we computed the risk score of each BC patient in the train group (n = 560) and assigned them into high-risk (n = 280) and low-risk groups (n = 280) with reference to the median risk score. The methylation levels of the six CpG loci between the high-risk group and low-risk group were showed in Supplementary Figures S2A. Kaplan-Meier curve indicated that the high-risk group had significantly poorer OS than the low-risk group (Figure 5A). To confirm the methylation risk score could effectively predict the BC patients’ prognosis, we plotted the ROC curve. Notably, we observed that the risk score had the highest prediction performance of prognosis compared with the conventional clinical features, with the 3-years and 5-years AUC values being 0.739 and 0.744 (Figures 5B,C). The relationship of methylation risk score, survival status and methylation levels of the six signature CpG loci was shown in Figures 5D–F. Univariate Cox analysis indicated that age, stage, T stage, N stage and risk score were significantly associated with OS. However, when they were introduced into multivariate Cox analysis, only age [hazard ratio, 1.026 (95% CI, 1.008–1.045), p = 0.005], and risk score [hazard ratio, 2.823 (95% CI, 2.131–3.741), p < 0.001] remained as independent prognostic predictors (Figures 5G,H). Similar prognostic significance was observed in the validation cohort (27k array dataset), with AUC values of 0.603 and 0.657 for 3 and 5 years, respectively (Supplementary Figures S2B–D).

Using the thresholds of adjusted p < 0.05 and | log2 (fold change) | >1, a total of 413 differentially expressed genes (DEGs) between methylation high-risk and low-risk groups were obtained. The volcano and heatmap visually displayed the DEGs (Figures 6A,B). To further investigate the biological characteristics of the DEGs, function and pathway annotations were performed. GO analysis indicated that these genes were involved in the regulation of cell cycle processes and mitotic cell cycle phase transition. KEGG analysis showed that these DEGs were mainly enriched in p53 and TGF-β signaling pathways (Figures 6C–F).

To explore the prognostic values of the identified 413 DEGs, we extracted the expression of these DEGs from the transcriptomic profiles of 986 BC patients in the TCGA-BRCA database for subsequent analysis. Firstly, 50 DEGs significantly related with OS were selected by Kaplan-Meier analysis (p < 0.05). All of these prognostic genes were subjected to univariate Cox regression analysis and 13 of them with p < 0.05 were further identified as prognosis-associated genes (Supplementary Table S1). Ultimately, six prognosis-associated genes, namely IRF2, KCNJ11, ZDHHC9, LRP11, PCMT1, and TMEM70, were included in developing a gene prognostic signature by multivariate Cox regression analysis. Among them, four genes (ZDHHC9, LRP11, PCMT1, TMEM70) were related with poor survival and two genes (IRF2, KCNJ11) were associated with good survival (Supplementary Figures S3). The risk formula was as follows:

Risk score = −(0.06128 × IRF2) − (0.0342 × KCNJ11) + (0.0228 × ZDHHC9) + (0.01257 × LRP11) + (0.01082 × PCMT1) + (0.02917 × TMEM70).

Gene expression analysis of the six signature genes in the Oncomine database (https://www.oncomine.org) revealed that ZDHHC9, LRP11, PCMT1, TMEM70 were highly expressed in breast cancer, and IRF2, KCNJ11 were highly expressed in normal tissues (Figure 7A). Then we explored the protein levels of these six genes between breast cancer and normal tissues in the Human Protein Atlas database (Uhlen et al., 2010) (HPA: https://www.proteinatlas.org/humanproteome/pathology). In accordance with the gene expression levels, the protein levels of ZDHHC9, PCMT1, TMEM70 were significantly higher in breast cancer, and the protein levels of IRF2, KCNJ11 were higher in normal tissues (Figure 7B). Moreover, we further checked the prognostic values of our six genes in the public database TCGA portal (version 1.0) (http://tcgaportal.org/TCGA/Breast_TCGA_BRCA/process.php), and we found that ZDHHC9, LRP11, PCMT1, TMEM70 were associated with poor prognosis, while IRF2 and KCNJ11 were related with good prognosis of BC patients (Figure 7C).

Perhaps not surprisingly, a positive and significant correlation was observed between the two prognostic signatures, which was mainly reflected at the following aspects: firstly, a moderate correlation was found between the six signature CpG loci and six signature genes. Specifically, the expression of IFR2 was positively related with the methylation of cg00945507, cg05406101, cg14992108, and cg18678121, above of which were all good prognostic factors, while IFR2 expression was negatively related with the methylation of poor prognostic CpG locus: cg13060154; And KCNJ11 expression was positively correlated with the methylation of favorable prognostic CpG loci: cg05406101 and cg10092957, and negatively related with the methylation of cg13060154; To the opposite of these two genes, ZDHHC9 and TMEM70 expression were negatively related with the methylation of cg05406101 and cg18678121, and positively associated with the methylation of cg13060154; Similarly, the expression of PCMT1 was negatively correlated with the methylation of cg05406101, cg14992108, cg18678121; And LRP11 expression was negatively related with the methylation of cg00945507, cg05406101, cg18678121 (Supplementary Figures S4A, Supplementary Figures S4I).

Subsequently, we examined the correlation between the six CpG loci and the gene risk score, and we observed that except for cg13060154 having the trend being positively correlated with the risk score, cg05406101 (R2 = −0.34, p < 0.001), cg18678121 (R2 = −0.28, p < 0.001), cg14992108 (R2 = −0.26, p < 0.001), cg00945507 (R2 = −0.26, p < 0.001), and cg10092957 (R2 = −0.11, p < 0.01) were all negatively correlated with the risk score (Supplementary Figures S4B–G, Supplementary Figures S4I ). Besides, we also explored the relationship between methylation risk score and gene risk score. Interestingly, we found that these two established risk scores were positively correlated with each other (R2 = 0.34, p < 0.001) (Supplementary Figures S4H,I).

The expression of the six signature genes between the gene high-risk and low-risk groups was shown in Supplementary Figures S5A. Kaplan-Meier analysis showed that survival probability in the low-risk group was higher (Figure 8A). The AUC values of 3-years and 5-years OS were 0.725 and 0.715 (Figures 8B,C). The risk score distribution, the survival status, and the expression of the six genes of 986 BC patients were visualized in Figures 8D–F. Univariate and multivariate Cox regression analyses indicated that the risk score was associated with OS and could be an independent prognostic predictor, with univariate hazard ratio, 1.187 (95% CI, 1.082–1.302, p < 0.001), multivariate hazard ratio, 1.199 (95% CI, 1.094–1.314), p < 0.001, respectively (Figures 8G,H).

To confirm the prognostic value of the six-gene signature, we tested it with the validation subgroup comprising of 557 BC patients. The results were consistent with our previous findings. Specifically, the OS in the high-risk group was poorer (Supplementary Figures S5B), and the AUC values of 3 and 5 years were 0.735 and 0.696 (Supplementary Figures S5C). Univariate and multivariate COX regression analysis also showed that the risk score was an independent prognostic predictor of BC (Supplementary Figures S5D,E). Subsequently, the 986 BC patients’ dataset was randomly assigned into two test subgroups [test group one (n = 494) and test group two (n = 492)] which were with balanced baseline characteristics (Table 4), and both of them were used for validating the gene signature. The two subgroups could also distinguish the favorable OS patients from the poor OS patients (Supplementary Figures S5F,G). AUC values of 1 year, 3 years, 5 years were 0.711, 0.757, 0.721 in test group one, and 0.864, 0.733, 0.702 in test group two, respectively (Supplementary Figures S5I,J). Moreover, the external dataset GSE146558 further confirmed our gene prognostic signature. In high-risk group of GSE146558 dataset, patients were with poorer OS (Supplementary Figures S5H). And AUC value of 3 years was 0.634 (Supplementary Figures S5K).

In addition, we confirmed the prognostic value of the six-gene prognostic signature in the subgroups of BC patients presented with different clinical features (age (<65 and≥65), T staging (T1-2 and T3-4), N staging [N0 and N1-3) and stage (stage I-II and stage III-IV)] (Supplementary Figure S6).

Considering the important roles of BRCA1, BRCA2, CDH1, PTEN, TP53, PIK3CA in BC, we also evaluated these gene expression between gene high-risk and low-risk groups, and observed that the expression of oncogenes such as BRCA1, BRCA2 and CDH1 were significantly higher in high-risk group. On the other hand, the expression of the tumor suppressor gene PTEN was significantly higher in low-risk group (Figure 9).

To investigate potential functions and signaling pathways related to the six-prognosis signature, we performed Gene Set Enrichment Analysis (GSEA: http://www.gsea-msigdb.org/gsea/index.jsp). Notably, we found that more tumor-related GO terms and KEGG pathways were associated with low-risk group (Figure 10A). In detail, low-risk group was mainly associated with the function of regulating epithelial and endothelial cell migration, and high-risk group was related with nuclear chromosome condensing, protein folding. Pathway enrichment analysis indicated that JAK/STAT signaling pathway, cell adhesion molecule signaling pathway, VEGF signaling pathway, and MAPK signaling pathway were active in the low-risk group. On the other hand, P53 signaling pathway was active in the high-risk group (Figure 10B).

Correlation analysis between the expression of the six prognosis genes and the sensitivity of anti-tumor drugs was performed based on the CellMiner database (https://discover.nci.nih.gov/cellminer/), and the results indicated that our signature genes were moderately correlated with the response of some common anti-tumor drugs such as PARP inhibitor (Olaparlib), chemotherapy drugs (Fluorouracil, Decitabine, Oxaliplatin), which might imply potential value in anti-tumor therapy (Figure 11).