Genome-wide discovery of circulating cell-free DNA methylation signatures for the differential diagnosis of triple-negative breast cancer

Patients and sample collection

This study recruited 113 primary BC patients diagnosed by pathology at Cancer Hospital of Harbin Medical University between May 2021 and December 2023. The exclusion criteria included: (1) incomplete molecular subtype data; (2) prior diagnosis of other malignancies; (3) any therapy received before surgery. Tumor tissues (n = 14; 5 TNBC and 9 non-TNBC) were obtained through surgical resection and immediately snap-frozen in liquid nitrogen. These samples were used for genome-wide DNA methylation profiling using the Infinium HumanMethylation850K BeadChip to identify candidate markers. Blood samples (1–2 mL, n = 113) were collected using Ardent Cell-Free DNA blood tubes and centrifuged at 800 g at room temperature for 10 min to separate plasma and buffy coat. Plasma was then transferred to a new tube and centrifuged at 16,000 g at 4 °C for 10 min to remove the remaining cells. Plasma and buffy coat were immediately stored at −80 °C until analysis. cfDNA extracted from the plasma of all enrolled patients (33 TNBC and 80 non-TNBC) was used to validate selected methylation markers through mddPCR assay. Demographic and clinical data were collected from medical records of all patients. This study was approved by the Medical Ethics Committee of Harbin Medical University (KY2016-01). Written informed consent was obtained from all patients before enrollment and sample collection.

Public data sources

The Illumina-normalized methylation data of tissues from the TCGA BC cohort, including 83 TNBC and 691 non-TNBC, were downloaded from UCSC Xena (https://xena.ucsc.edu/public/). Additionally, Illumina-normalized methylation data from 342 BC tissues and 233 BC white blood cells (WBC) were downloaded from GEO database (GSE69914, GSE72251, and GSE50132). The methylation level for each CpG was represented as a beta value (β), which is a ratio of intensities between the methylated allele and the sum of M and unmethylated, ranging from 0 (no methylation) to 1 (full methylation).

DNA extraction and bisulfite conversion

Tumor tissue DNA was extracted using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) and quantified using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). Circulating cfDNA was extracted from plasma using the QIAamp® Circulating Nucleic Acid Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions, with concentrations measured using a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Both cfDNA and genomic DNA underwent bisulfite conversion using the EZ DNA Methylation-Gold™ Kit (Zymo Research, Irvine, CA, USA), according to the manufacturer’s guidelines.

Tissue methylation markers discovery and selection

As shown in Fig. 1, we first performed the genome-wide methylation analysis on five in-house TNBC and nine non-TNBC tissues (Table S1) using Infinium MethylationEPIC BeadChip. After data filtering, correction, and normalization, CpG sites with |Δβ| > 0.10 and P < 0.05 were identified as DMCs. Given the limited sample size of the in-house dataset, a more lenient threshold was used to include more potential candidate markers. To improve reliability, these markers were further evaluated in a larger cohort of 83 TNBC and 691 non-TNBC tissues from the TCGA dataset, using a stricter cutoff (|Δβ| > 0.25 and P < 0.05). Subsequently, to minimize false positive detection and increase the signal-to-noise ratio in cfDNA analysis, CpG sites that were almost completely methylated (average β value > 0.90) or unmethylated (average β value < 0.10) in WBC were selected. Pearson’s correlation method was used to calculate correlation for methylation markers. Finally, considering the existence of multicollinearity among them, the LASSO analysis with 10-fold cross-validation was performed to select the optimal markers through the “glmnet” package (Friedman, Hastie & Tibshirani, 2010) in the TCGA dataset. The eight tissue methylation markers with non-zero coefficients were selected for further analysis.

Evaluation and validation the diagnostic effect of DNA methylation markers in tissues

The diagnostic accuracy of individual CpG sites in differentiating TNBC and non-TNBC tissues was assessed using ROC curves. CpG sites with an AUC greater than 0.80 were retained. An eight-CpG methylation diagnostic score (MDS) was then constructed using logistic regression in TCGA dataset and validated in GSE69914 dataset. It was calculated based on the corresponding coefficients of the selected methylation markers. The formula was as follows:

$$\mathrm{M}\mathrm{D}\mathrm{S}=\sum _{i=1}^n\left(m_i\ast Coe{f}_i+b\right).$$

Here, n represents the number of methylation markers; $m_i$ represents the methylation level of each marker; $Coe{f}_i$ represents the coefficients; $b$ represents to the intercept; and MDS represents a weighted sum of the methylation level of each marker.

Evaluation and validation the prognostic effect of DNA methylation markers in TNBC tissues

The primary endpoint for prognostic analysis in the study was overall survival (OS), defined as the time from surgery to death from any cause or the last follow-up visit. In addition, disease-free survival (DFS) was evaluated as a secondary endpoint, defined as the time from surgery to the first local recurrence, distant metastasis, or death from any cause. For each candidate CpG site, patients were divided into low-risk and high-risk groups based on median methylation levels. Univariate Cox regression analysis was performed to screen out CpG sites that significantly associated with OS or DFS in the TCGA-TNBC patients (P < 0.10). Multivariate Cox regression was subsequently performed, incorporating the significant CpG sites as predictor variables and adjusting for age and stage, to identify independent prognostic markers. The stability of these markers was validated in the GSE72251-TNBC cohort.

Development of mddPCR assay

We successfully designed forward and reverse primers, as along with minor groove binder (MGB) Taqman probes (FAM^TM/VIC^TM-reporter dyes) for two of the targeted markers (Table S2), while design attempts for the others were unsuccessful due to unfavorable sequence features. The methylation-insensitive β-actin (ACTB) gene served as an internal control in each PCR well (primers and probes reported elsewhere) (Zhang et al., 2023). A mddPCR assay (including three genes) was then performed in a 21 μL final volume system with 10 μL of ddPCR™ Supermix for Probes (No dUTP), adjusted volumes of primers and probes, and 5–6 μL bisulfite-converted DNA. The QX200™ Droplet Generator (Bio-Rad, Hercules, CA, USA) generated droplets and the PCR conditions were as follows: 10 min at 95 °C, 40 cycles of 94 °C for 30 s, and annealing and extension at 60 °C for 1 min; 10 min hold at 98 °C. The ramp rate was set at 2 °C/s for all steps. The emulsions were analyzed on the QX200^TM Droplet Reader device to count droplets containing amplified DNA targets and empty droplets based on fluorescence. Data were analyzed using the QuantaSoft^TM Analysis Pro software v1.0. Positive droplets (containing an amplified DNA target) were manually identified in a 2D plot based on the distribution of droplets of control samples (methylated DNA control, unmethylated DNA control, RNase-free water control, and non-template control). cfDNA methylation was quantified as the number of copies of methylated alleles per 1 mL of plasma.

Statistical analysis

Categorical variables were described as the number (percentage), and continuous variables as median (interquartile range (IQR)) or mean (standard deviation (SD)). Methylation levels between the two groups were compared using the t-test, Wilcoxon test, or Mann-Whitney U test. Differential methylation analysis between TNBC and non-TNBC tissues were conducted using ChAMP (2.32.0) package (Tian et al., 2017) or “Limma” package (Ritchie et al., 2015). The unsupervised hierarchical clustering heatmap of DMCs methylation was generated using the “pheatmap” package (Kolde, 2019). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (Gene Ontology Consortium, 2021; Kanehisa et al., 2021) were conducted using the “clusterProfiler” package (Wu et al., 2021). Sensitivity, specificity, and AUC were used to evaluate the diagnostic performance of the individual CpG sites and models. The R package “survival” (Therneau, 2024) and “survminer” (Kassambara, Kosinski & Biecek, 2021) were used for survival analysis. The Kaplan-Meier and log-rank tests were utilized to generate the survival curves and compare the differences between the groups. Time-dependent ROC curves were performed to evaluate the prognostic performance of CpG sites by the R package “timeROC” (Blanche, Dartigues & Jacqmin-Gadda, 2013). Hazard ratios (HR) and 95% confidence interval (CI) were calculated with the Cox regression analysis. A minimum sample size of 31 in each group could achieve 80% power with a significance level of 0.05 using a one-sided z-test, which was conducted with PASS 11.0 (UCSS, USA) (Hanley & McNeil, 1983; Obuchowski & McClish, 1997). Statistical analyses were performed in R (version 4.3.1), and P < 0.05 was considered significant.

Basic characteristics of patients

A total of 33 TNBC and 80 non-TNBC patients were enrolled in this study. The age distribution was relatively balanced between TNBC and non-TNBC groups (mean age, 59.5 years vs. 58.7 years, P = 0.800). Detailed characteristics of patients are summarized in Table 1 and Table S3.

Differential methylation markers based on in-house primary tissues

Based on a genome-wide methylation analysis between in-house five TNBC and nine non-TNBC tissues using Infinium HumanMethylationEPIC BeadChip, a total of 32,787 DMCs (|Δβ| > 0.10 and P < 0.05) were identified; with 9,294 (28.35%) DMCs higher methylation levels in TNBC tissues (defined as hypermethylated CpG) and 23,493 (71.65%) with lower methylation levels (defined as hypomethylated CpG, Fig. 2A). An overview of the relative CpG island position and gene location annotation were shown in the Sankey plot (Fig. 2B). Pathway enrichment analysis showed that the DMCs were mainly enriched in signaling-related pathways (Fig. 2C) and the biological process of protein binding (Fig. S1).

Discovery and selection of tissue-based methylation markers

Based on above obtained 32,787 DMCs from the in-house dataset, the methylation data of 83 TNBC and 691 non-TNBC tissues from TCGA were compared and selected 1,130 CpG that were still significantly differentially methylated (|Δβ| > 0.25 and P < 0.05). Since cfDNA can be derived from WBC of cancer patients, we excluded the CpG sites with average methylation β ≥ 0.10 or β ≤ 0.90 in 233 WBC to reduce the false positive possibility, and we obtained 113 DMCs. The 113 methylation markers showed high correlations among them (Fig. S2). To further determine CpG sites that could distinguish TNBC from non-TNBC, we applied the Lasso analysis to select optimal markers in TCGA dataset. We obtained eight markers (cg19758859, cg01095157, cg14534279, cg17588293, cg06268921, cg04016621, cg23247845, and cg02096552) with non-zero coefficients for further analysis, and the optimal λ value = 0.0234, log (λ) = −3.7547 (Table 2 and Fig. 3). The methylation levels for two CpG sites were significantly higher in TNBC than those in non-TNBC tissues, six residual CpG sites were significantly lower in TNBC than those in non-TNBC tissues (Fig. 2D and Table S4).

Diagnostic performance of candidate markers in tissues

We investigated the diagnostic performance of eight CpG sites, with AUCs ranging from 0.809 to 0.893 for differentiating TNBC from non-TNBC tissues (Table S5). These CpG sites (AUC > 0.80) were used to construct the MDS through logistic regression analysis. The MDS was a sum of the methylation levels of eight CpG sites, each multiplied by its corresponding coefficient from the logistic regression analysis: MDS = cg19758859 × (−0.394) + cg01095157 × (−1.894) + cg14534279 × (−1.088) + cg17588293 × 1.317 + cg06268921 × 1.207 + cg04016621 × 0.126 + cg23247845 × (−1.559) + cg02096552 × (−1.489) + 1.189. The MDS was significantly different between TNBC and non-TNBC tissues. (Fig. 4A), with an AUC of 0.922 (95% CI = 0.895–0.950) in the TCGA dataset (Fig. 4B). Using the cutoff value of −2.51, the MDS yielded a sensitivity of 93.8% (61/65) in stage I/II, 82.4% (14/17) in stage III/IV at a specificity of 84.9% (587/691, Table S6). Furthermore, there was no significant difference in MDS between stage I/II and stage III/IV TNBC tissues (Fig. S3A). Notably, the MDS could discriminate the early-stage TNBC from non-TNBC, with an AUC of 0.932 (95% CI [0.908–0.957], Fig. S3B).

The MDS was then applied to the GSE69914 validation dataset, where it also showed a significant difference between TNBC and non-TNBC tissues (Fig. 4C). Using the same cutoff value of −2.51, the MDS achieved an AUC of 0.875 (95% CI [0.789–0.961], Fig. 4D), with a sensitivity of 86.70% (26/30) and a specificity of 90.2% (248/275).

Prognostic effect of candidate markers in tissues

Next, we investigated prognostic potential of the eight CpG sites in the TCGA-TNBC cohort by analyzing their association with OS and DFS. Among them, only cg06268921 (≥0.50 vs. <0.50) was significantly associated with both OS and DFS in univariate Cox regression analysis (P < 0.10; Tables S7 and S9). After adjusting for age and stage in multivariate Cox regression analysis, cg06268921 remained significantly associated with OS (HR: 0.249, 95% CI [0.064–0.966], P = 0.044) and DFS (HR: 0.194, 95% CI [0.052–0.727], P = 0.015), suggesting that it served as an independent prognostic factor for TNBC (Tables S8 and S10). The Kaplan–Meier survival curves further confirmed that patients in the low-methylation group (<0.50) had significantly worse OS and DFS compared with those in the high-methylation group (≥0.50; log-rank P < 0.10, Figs. 5A, 5B). The prognostic performance of cg06268921 was evaluated using time-dependent ROC analysis, yielding AUCs of 0.521 and 0.422 for 3-year OS and DFS, respectively (Figs. S4A, S4B). To validate our findings, we analyzed the prognostic effect of cg06268921 in an independent external cohort (GSE72251). The association was not statistically significant, possibly due to limited sample size and cohort heterogeneity (Tables S8 and S10). Meta-analysis of both cohorts showed non-significant pooled HR for OS and DFS, with moderate heterogeneity (I² > 60%), indicating variability in the prognostic effect of cg06268921 across cohorts (Fig. S5).

Development of cfDNA mddPCR assay and the evaluation of limit of quantification

We constructed an mddPCR assay (cg06268921, cg23247845, and ACTB, Fig. 6A and Table S11) and then optimized its annealing temperature. By comparing the fluorescent of methylated and unmethylated control in the FAM channel, we found that an annealing temperature of 58.8 °C yielded the best performance (Fig. S6).

Next, we compared the limit of quantification (LOQ) of mddPCR assay with multiplex quantitative methylation-specific PCR (mqMSP) assay using 10 ng of mixed methylated DNA control and the unmethylated DNA control at concentrations of methylation at 100%, 10%, 5%, 1%, 0.5%, 0.1%, 0.01%, and 0. Three replicates of each reaction were performed. The mddPCR assay could simultaneously and independently quantify the methylation level of each target gene, while mqMSP could only determine the combined methylation level of multiple genes. In the FAM channel of the mqMSP assay, one methylated allele in a background of 20 unmethylated alleles was detected (LOQ = 5%, R² = 0.982). In contrast, the LOQ of cg06268921 (R² = 0.790) and cg23247845 (R² = 0.990) of mddPCR were 0.01% and 0.1%, respectively, indicating that the LOQ of mddPCR assay was 0.01%, and greater sensitivity than the mqMSP assay (Fig. S7).

Validation the diagnostic effect of methylated markers detected with mddPCR in cfDNA cohort

Moreover, we applied the mddPCR assay to detect methylation levels in cfDNA from 33 TNBC and 80 non-TNBC patients. The TNBC and non-TNBC yielded the median cfDNA concentration of 7.15 ng/mL (IQR: 4.58–9.80), and 6.98 ng/mL (IQR: 2.54–10.4), respectively, with no statistically significant different (P = 0.97, Fig. S8). However, the number of methylated molecular copies of cg06268921 was significantly higher in TNBC than in non-TNBC subtypes with AUC of 0.719 (95% CI [0.614–0.823]), but cg23247845 was not significantly different between two groups (Fig. S9). We incorporated these two markers to construct the cfDNA-based MDS (cf-MDS), which was significantly different between TNBC and non-TNBC groups in cfDNA and achieved a sensitivity of 54.5%, a specificity of 82.5%, and an AUC of 0.728 (95% CI [0.618–0.839], Figs. 6B, 6C). In addition, the cf-MDS could distinguish stage I/II TNBC from non-TNBC, with an AUC of 0.701 (95% CI [0.578–0.825], Fig. S10).