Paired DNA Methylation chip data (Illumina Infinium 450 K methylation microarray data) and gene expression of BRCA patients were downloaded from TCGA database, comprising 447 tumor and 52 paired normal samples. The corresponding clinical data including cancer stages, survival time, and tumor purity, was also downloaded. Another DNA methylation data based on peripheral blood were obtained from GEO database under accession number (GSE214344).

On the basis of our previously constructed method for identifying tumor-specific methylation markers (Qin et al., 2024), we screened differential methylation sites (DMPs). Firstly, we performed probe quality control by: 1) removing probes with more than 90% missing values in samples, 2) excluding sites located on the X and Y chromosomes, and 3) eliminating non-probe sites. Ultimately, we obtained 383,561 sites. Subsequently, we selected tumor tissues with paired samples and with the tumor purity greater than 0.6 (estimated using the ESTIMATE method). This resulted in a final cohort of 57 BRCA patients with paired tumor and normal samples. Finally, we employed a nonparametric permutation test to compare the paired samples in this cohort. The model is as follows: β ij = β 0 + β 1 × TY i + TO i + ε ij β ij is the methylation level (0,1)of the j th probe of the i th sample, TY i is the tissue type of the i th sample (cancer or paracancerous tissue), β 0 is the intercept, β 1 is the coefficient related to tissue type, ε ij is Gaussian error.

According to the aforementioned model, we performed 1,000 permutations for each methylation site to obtain the corresponding p-value for the methylation level between cancer tissues and adjacent normal tissues. Probes with a p-value less than 0.001 were defined as BRCA associated DMPs. Then, we calculated the difference value of methylation level (∆β = β_cancer - β_normal) for each sample in the cohort of 57 paired samples. Probes with ∆β > 0 in 80% of the samples were defined as putative hypermethylated positions (Hyper-DMPs), while probes with ∆β < 0 in 80% of the samples were defined as putative hypomethylated positions (Hypo-DMPs). This resulted in a total of 51,736 putative Hyper-DMPs and 27,405 putative Hypo-DMPs.

To further determine the methylation level of each probe, we calculated the median ∆β value and ( median ∆ β − ∆ β value for the 57 samples. Moreover, to better control unspecific variation of the ∆ β values, we set up a second-step threshold. Detailly, based on the distribution of median − ∆ β values, we selected the top 5% of the 51,736 Hyper-DMPs, threshold of 0.4124, as the cutoff for identifying final Hyper-DMPs, and the bottom 5% of Hypo-DMPs (n = 27,405), the threshold of −0.4125, to identify final Hypo-DMPs (Supplementary Figure S1).

We performed methylation signature profiling in BRCA patients according to previously published methods (Qin et al., 2024). Firstly, we selected a BRCA population from the TCGA database with tumor tissue purity greater than 0.6 (N = 447). Then, we applied non-negative matrix factorization (NMF) (Lee and Seung, 2000) on the Hyper-DMPs and Hypo-DMPs from the 447 samples. We used two different algorithms, nsNMF (using Kullback-Leibler divergence for multiplicative updates) and Lee (based on Euclidean distance), to perform matrix factorization. Each matrix was iterated 100 times. NMF is an unsupervised learning algorithm that decomposes a non-negative matrix into two matrices, show as following. M p × n = E p × k × S k × n

M is a p × n matrix of methylation, where p = 2 , 587   Hyper − DMPs ;   p = 1 , 370 Hypo − DMPs , and n = 447 . S is a matrix of k methylation signatures for n tumor samples, and matrix E is a weight matrix of p probes for k methylation signature.

The optimal number of labels was selected based on the cophenetic coefficient and average silhouette widths. Finally, we obtained three high methylation labels (Hyper-S1, Hyper-S2, and Hyper-S3) and five low methylation labels (Hypo-S1, Hypo-S2, Hypo-S3, Hypo-S4, and Hypo-S5). NMF classified the probes based on their contribution to the labels’ weights, including: Hyper-DMPs:890 probes in Hyper-S1, 855 in Hyper-S2, 842 in Hyper-S3; Hypo-DMPs: 335 in Hypo-S1, 287 in Hypo-S2, 316 in Hypo-S3, 188 in Hypo-S4, and 244 in Hypo-S5.

Several evaluation indices were used to assess the classification performance of methylation probes for immunotherapy response, including sensitivity, specificity, precision, and accuracy. Sensitivity refers to the ability of the probes to identify positive events (in this case, immune therapy response). Specificity refers to the ability to identify negative events (non-response). Precision considers the false positive rate, defined as the number of incorrectly predicted responses divided by the total number of predicted responses. Accuracy measures the proportion of true results (both true positives and true negatives) among the total number of cases examined. Four Sensitivity = TP TP + FN Specificity = TN TN + FP Precision = TP TP + FP Accuracy = TP + TN TP + TN + FP + FN

TP: True positive rate; TN: True negative rate; FP: False positive rate; FN: False negative rate.

Furthermore, the predictive abilities of the DNA methylation probes were estimated using under the curve (ROC AUC). The calculation formula of the AUC is: A U C = ∫ 0 1 T R P F P R − 1 t d t

In this formula, TRP is the true positive rate (sensitivity); FPR-1(t) is the inverse function of the false positive rate (1-specificity), representing the false positive rate corresponding to the true positive rate t; the integral range is from 0 to 1, covering the entire range of the ROC curve.s.

To evaluate the correlation between DMP signatures and clinical factors, we compared signature activity differences between primary (stages I and II) and advanced tumors (stages III and IV) using the Wilcoxon signed rank-sum test. Patients were assigned to high- and low-signature groups based on the median signature activity as a cutoff. Differences in overall survival between high- and low-signature groups were evaluated using the Cox proportional hazards regression model via the “survival” R package and visualized with the “ggsurvfit” package.

Programmed death-ligand 1 (PD-L1) regulates T cell exhaustion by binding programmed death-1 (PD-1) on T cells. Cancer cells with high amounts of PD-L1 can turn T cells off and inhibit T cells attacks, against which immunotherapy medicines (immune checkpoint inhibitors, ICIs) may be effective. PD-L1 is a widely used predictive biomarker for immunotherapy response (Teng et al., 2018). Higher expression of PD-L1 typically correlate with greater therapeutic benefit from ICIs (Jiang et al., 2019). To validate the immunotherapy predictive power of identified methylation signatures, we assessed the relationship between signature activities and PD-L1 gene expression (Transcript per million, TPM) based on paired BRCA DNA methylation data and RNA transcriptome dataset from the TCGA database using cor.test.

Genes containing DMPs in the methylation signatures were extracted. Gene ontology (GO) and KEGG pathway enrichment analysis were performed on these genes using Enrichr (Chen et al., 2013).

Statistical analysis and plotting were performed using R software. All statistical significances were set as p or adjusted p-value <0.05.

Based on the 450K DNA methylation chip data from TCGA-BRCA, a total of 92 paired tumor samples and 383,561 methylation sites were selected after quality control. The permutation test method screened 148,494 (38.71%) differential methylation sites (DMPs). According to the 5% cutoff values of Hyper-DMPs (0.4124, n = 2,587) and Hypo-DMPs (−0.4125, n = 1,370) (Supplementary Figure S1), we finally obtained 2,587 Hyper-DMPs and 1,370 Hypo-DMPs. The amount of Hyper-DMPs is almost twice the number of Hypo-DMPs.

Given the non-negativity of DNA methylation levels, the NMF method was subsequently used to extract DNA methylation signatures. Based on the cophenetic coefficient and average silhouette widths of NMF method (Figure 1A), we obtained three high methylation signatures (Hyper-S1, Hyper-S2, and Hyper-S3) and five low methylation signatures (Hypo-S1, Hypo-S2, Hypo-S3, Hypo-S4, and Hypo-S5) from 447 BRCA patients (Figure 1B).

To analyze the genomic characteristics of methylation signatures, we analyzed the genomic distribution features of DMPs associated with each methylation signature. We found that those DMPs associated with high methylation signatures (Hyper-) were mostly located in promoter regions, while DMPs associated with low methylation labels (Hypo-) were mainly located in gene body regions (Figure 2). In terms of CGI position, the Hyper-DMPs were mostly located in CpG islands, the Hypo-DMPs were mostly located in Open Sea areas (Figure 2).

Since the DNA methylation signatures were tumor-related, whether they were related to tumor progression remained to be answered. To clarify this question, we performed correlation analysis between DNA methylation signature and tumor stage including primary (stage I-II) and advanced (stage III-IV) stages. The Wilcoxon’s rank sum text showed that four out of eight (50%) DNA methylation signatures exhibited markedly differences between the two stages. Detailly, Hyper-S3 and Hypo-S4 were extremely higher in advanced stages. However, Hyper-S1 and Hypo-S3 were significantly upregulated in the primary stages of tumors. The distinct pattern between these four signatures and tumor progression indicates that these four signatures might act vital roles in tumor development (Figure 3A).

DNA methylation alterations can influence tumor development and patient outcomes. To assess the impact of the methylation signatures on BRCA patient’s survival, we compared overall survival times between groups stratified by high and low signature activity levels. Our results showed that patients with low activities of Hyper-S2 signature exhibited significantly longer survival times than those with high Hyper-S2 activities (Figure 3B), indicating its prognostic value. Unlike the tumor stage, only a small proportion (1/8) of signatures affected overall survival. This implies DNA methylation may have a greater impact on tumor progression than direct effects on the survival of BRCA. We therefore focused subsequent analyses on the signatures associated with tumor stage, as these are likely to provide further mechanistic insights into BRCA pathogenesis.

To clarify whether the probes of methylation signatures in tumor tissues are consistent with the methylation level in circulating cell-free DNA of BRCA patients, we further mined the four methylation signatures (Hyper-S1 = 890, Hyper-S3 = 842; Hypo-S3 = 316, Hypo-S4 = 188) related to tumor stage in the methylation data of peripheral blood from five healthy individuals and eight breast cancer patients. We analyzed the probes with abnormal methylation in BRCR peripheral blood within these four methylation labels, and found that 510 probes were associated with tumor blood abnormal methylation (p < 0.05), including 137 sites in Hyper-S1, 307 in Hyper-S3, 42 in Hypo-S3, and 24 in Hypo-S4. Notably, we found that the 24 significantly different probes in Hypo-S4 all showed significantly lower methylation in BRCA patients than in healthy individuals (Figure 4), suggesting that these probes could serve as molecular markers for early breast cancer diagnosis.

To further validate the relationship between these 510 probes with abnormal methylation and immune therapy efficacy, we performed AUC analysis. Our results showed that 23 probes were associated with patients’ immune therapy response, including 11 probes in Hyper-S1, 7 in Hyper-S3, 2 in Hypo-S3, and 3 in Hypo-S4 (Figure 5A). Their highest predictive power for immune therapy response was 86% (AUC, Table 1). Among them, cg00870269 had a sensitivity of 93%, precision of 86%, accuracy of 74%, and specificity of 66%. cg10017626 had a sensitivity of 75%, precision of 86%, accuracy of 86%, and specificity up to 92%. KEGG pathway enrichment analysis revealed that genes associated with these 23 probes are mainly involved in glycan degradation, Type I diabetes mellitus, cholesterol metabolism, notch signaling pathway, and glycerolipid metabolism pathways (Figure 5B, p < 0.05). GO biological process enrichment showed that the immunotherapy-related methylated genes are primarily enriched in processes including regulation of Type B pancreatic cell proliferation, transmembrane receptor protein tyrosine phosphatase signaling pathway, cleavage furrow formation, insulin secretion involved in cellular response to glucose stimulus, protein deglycosylation, and very-low-density lipoprotein particle remodeling (top five, Figure 5C).

The PD-L1 is a biomarker used to assess the effectiveness of ICIs. Therefore, to validate the predictive utility of these signatures and their 23 probes for immunotherapy efficacy, we evaluated the correlation between PD-L1 gene expression and methylation levels from the TCGA-BRCA dataset. Notably, results showed that the activity of the Hypo-S4 signature was positively correlated with PD-L1 expression (Figure 6A). Consistently, three probes in Hypo-S4 associated with immunotherapy response were markedly and positively correlated with PD-L1 gene expression (Figure 6B). While activities of the other three stage-related signatures showed no significant relationship with PD-L1 expression overall (Figure 6A), several individual probes from their corresponding signatures did demonstrate marked correlations, including cg07517893 and cg18419175 from the Hyper-S1 and cg10881128 and cg18530645 from the Hyper-S3 (Figure 6B).