Fragmentomic liquid biopsy enables early breast cancer detection, molecular subtyping and lymph node assessment

Domain: chemistry
Published: 2026-03-06
Editorial outcome: minor_revision

Abstract
Breast cancer remains a leading global health concern in women, while screening is still limited by imaging accessibility and reduced sensitivity in dense breasts. Here we conduct a multicenter case-control study including 503 breast cancer patients and 289 benign controls to develop TuFEst, a machine learning model based on genome-wide cell-free DNA fragmentomic features. TuFEst achieves high sensitivity (95%) and specificity (78.3%) for early cancer detection and reliably identifies malignancies missed by conventional imaging. Extension of this framework enables non-invasive molecular subtyping (TuFEst-MS) and lymph node status prediction (TuFEst-LN), with strong performance in independent validation cohorts and imaging-pathology discordant cases. Transcriptomic profiling of paired bulk tumor samples (n=79) demonstrates that elevated TuFEst-derived cancer scores reflect tumor aggressiveness and immune-related biological programs. Together, these findings support cfDNA fragmentomics as an integrated liquid biopsy strategy for breast cancer management, enabling concurrent detection, molecular subtyping, and lymph node evaluation with potential clinical utility.

Introduction
Breast cancer (BC) is the most common malignancy and a leading cause of cancer-related mortality in women worldwide, accounting for approximately 30% of female cancers and 15% of cancer deaths1. In China, its incidence has surged, with new cases rising from 3 million in 2015 to 4.2 million in 2020, now ranking as the fourth most prevalent cancer overall2.

Early detection is critical to reduce mortality. The U.S. Preventive Services Task Force strongly endorses population-level BC screening because of its proven survival benefits3. However, global screening practices vary widely. In the U.S., annual or biennial mammography is standard for women aged 40-74, with intervals tailored to individual risks3. In contrast, China’s “Two Cancers” program (BC and cervical cancer) recommends triennial to quinquennial screenings, contingent on local resources. During 2018–2019, only 22.3% of eligible Chinese women underwent BC screening4. This limited coverage stems from both structural and clinical factors. On the infrastructure side, many primary-level hospitals lack mammography equipment, which creates significant access barriers. On the clinical side, while national guidelines favor ultrasound over mammography for better detection of typically dense breast tissue in Asian women4, this preference comes with operational challenges. Diagnostic inconsistencies lead not only to missed cancers (higher interval rates) but also to overtreatment, straining healthcare insurance funds. These disparities directly impact diagnosis timing: 24.1% of Chinese BC patients are diagnosed at Stage III+ versus 11.8% in the U.S., which significantly impairs the prognosis and increases treatment costs2.

These challenges underscore the urgent need for standardized operator-independent screening solutions. Liquid biopsy approaches, particularly cell-free DNA (cfDNA) analyses, offer a promising opportunity to overcome the current limitations. Unlike imaging modalities, cfDNA-based detection requires minimal technical expertise for sample collection and processing, making it particularly suitable for implementation in diverse healthcare settings. Current liquid biopsy approaches for BC screening – including mutation analyses and methylation profiling – face practical clinical limitations for early-stage disease. Mutation-based tests often require large blood volumes and ultra-deep sequencing because variant allele fractions are extremely low in early cancers, which limits sensitivity and inflates cost and complexity5. Meanwhile, methylation-based assays achieve high specificity but only modest stage-I sensitivity in large studies6, and typically rely on computationally intensive pipelines which restrict their widespread implementation7. In contrast, fragmentomic analyses, which profile patterns in cfDNA, including size distributions, end motifs, and fragmentation patterns, provide an alternative method that circumvents these limitations. This approach requires minimal blood volume (≤5 mL), enables high-throughput processing, and preserves tissue-specific signatures, making it ideal for population screening. While fragmentomic methods like DELFI (DNA evaluation of fragments for early interception) show multi-cancer detection potential8, early BC studies remain limited by small sample sizes, inconsistent methodologies, and underrepresentation of early-stage cases. Large-scale multicenter validation is needed to establish standardized protocols and to demonstrate their clinical utility. Recent studies have shown that cfDNA fragmentation patterns are shaped not only by nucleosome positioning but also by transcription factor (TF) binding and repeat element–associated chromatin architecture9. These signals provide complementary layers of biological information, enhancing early cancer detection beyond what can be achieved with size-based fragmentation alone. Motivated by these findings, our study integrates multiple fragmentomic modalities—including TF motif–associated fragmentation, repeat-element instability, and genome-wide nucleosome footprints—into a unified machine learning framework for early BC detection.

Here, we conducted a multicenter case-control study (NCT06016790) across diverse clinical settings, developing TuFEst (Tumor Fraction Estimator), a machine-learning framework that leverages ultra-low-coverage whole-genome cfDNA sequencing to quantify tumor burden through fragmentomic analyses. To the best of our knowledge, this represents the largest and most comprehensive investigation of fragmentomics for early-stage BC, systematically evaluating its clinical utility for screening while extending applications to molecular subtyping and lymph node metastasis prediction. Our integrated approach using a single blood test represents a potential shift toward more comprehensive and personalized BC management via cfDNA fragmentomics.

Results
Study design and clinical characteristics
Our multicenter case-control study enrolled 503 patients with histologically confirmed malignant BC (77.1% early-stage [0-IIa], n=388) from six tertiary Grade-A hospitals and one community healthcare center, along with 289 benign controls (Supplementary Table 1). Following standardized collection protocols, we obtained peripheral blood samples from all participants, isolated approximately 1 mL of plasma per case, and performed low-pass whole-genome sequencing (WGS) at 2× average coverage. All blood samples were collected, processed, and transported in accordance with uniform standard operating procedures across participating centers to ensure pre-analytical consistency. To investigate the clinical utility of cfDNA fragmentomics in early BC characterization, we developed a machine learning framework capable of simultaneously discriminating malignant from benign breast lesions, classifying molecular subtypes, and predicting lymph node metastasis status (Fig. 1).

cfDNA Fragmentation Reflects Transcription Factor Activity in Early BC
Through a comprehensive evaluation of genome-wide cfDNA fragmentation patterns across all chromosomes, we identified distinct chromatin accessibility signatures correlated with disease progression. While individuals with benign breast lesions demonstrate remarkably consistent fragmentation profiles, invasive BC cases exhibited progressively more disorganized fragmentation landscapes with increasing tumor stage. This phenomenon was detectable in Stage I tumors and became more pronounced in advanced stages (Stage II-III) and metastatic disease (Stage Ⅳ), where widespread fragmentation abnormalities were observed across multiple chromosomal regions (Fig. 2A). This heterogeneity reflects the accumulation of genomic alterations during the development and progression of BC.

Because cfDNA fragmentation reflects nucleosome positioning and TF binding, we reasoned that altered TF activity may imprint detectable footprints on plasma DNA. Building on prior fragmentomic studies showing that TF-site coverage enhances detection and biological interpretability9,10, we hypothesized that TF binding alterations contribute to early-stage fragmentation patterns. By combining chromatin immunoprecipitation sequencing (ChIP-seq) data from the ReMap 2020 database with our cfDNA fragmentation profiles, we performed a comprehensive comparison of TF binding site coverage between patients with BC and benign controls. Using stringent filtering criteria (read count > 25, read coverage ≥ 0.9, fold change ≤ 0.5, or ≥ 2), we identified 146 differentially fragmented genomic regions, corresponding to 844 potential TF regulators (Fig. 2B). Subsequent validation in BC cell lines (MCF-7 and MDA-MB-231; ReMap 2020) refined this set to 148 high-confidence TFs with conserved binding patterns in malignant cells (Fig. 2C).

KEGG pathway enrichment analyses revealed that these selected TFs were significantly involved in key BC-related processes, including tumor pathogenesis, endocrine therapy resistance, estrogen signaling, and transcriptional dysregulation – all hallmarks of malignant progression (Fig. 2D). Among the most significantly altered TFs are several with well-characterized roles in BC biology: FOXA1 (Forkhead Box A1), a pioneer factor critical for ERα-positive BC proliferation and endocrine resistance11; AR (Androgen Receptor), which promotes tumor progression and therapeutic resistance particularly in triple negative breast cancer (TNBC)12; HDAC2 (Histone Deacetylase 2), associated with advanced disease and poor outcomes13; RELA (RELA Proto-Oncogene NF-κB Subunit), the central NF-κB subunit driving inflammatory responses, tumor proliferation, and metastasis14; YBX1 (Y-box Binding Protein 1), a known promoter of Epithelial-Mesenchymal Transition (EMT), metastasis, and drug resistance15; ESRRA (Estrogen-Related Receptor Alpha), involved in metabolic reprogramming16; and SMARCE1 (SWI/SNF Chromatin Remodeling Factor), a chromatin remodeler frequently dysregulated in metastatic BC17. Conversely, TFs showing reduced fragmentation signals included tumor suppressors and differentiation factors, such as REST (RE1 Silencing Transcription Factor)18, CREB1 (cAMP Responsive Element Binding Protein 1)19, and CEBPB (CCAAT Enhancer Binding Protein Beta)20, which maintain normal breast tissue homeostasis and oppose malignant transformation (Table 1).

These findings revealed that cfDNA fragmentation patterns capture dual clinical insights: transcriptionally active programs driving early breast carcinogenesis and machine-learning-ready biomarkers, with the identified TF-associated signatures linking epigenetic reprogramming to actionable targets for therapy and liquid biopsy monitoring.

TuFEst model for BC detection
We developed a machine-learning framework to evaluate whether cfDNA fragmentomic features could reliably distinguish malignant from benign breast lesions. We systematically compared six algorithms: Generalized Linear Model (GLM), Support Vector Machine (SVM), eXtreme Gradient Boosting (XGBoost), Multilayer Perceptron (MLP), Gradient Boosting Machine (GBM), and Stacked Ensemble Model (SEM). The TuFEst model employs a comprehensive multi-modal approach to cfDNA fragmentomic analyses, integrating (1) repeat element fragmentation (Alu, Alternating locus units; MER, Medium Efficiency Repeats; LTR, Long Terminal Repeats; families and DNA_TcMar-Tigger) derived from WGS21; (2) quantitative fragmentation indices including the N-index and E-index for DNA end characterization22; and (3) TF binding site coverage analyses capturing chromatin accessibility signatures23,24. Multi-dimensional features – encompassing end motifs (4bp), breakpoint (6bp), fragment-size distributions, and regulatory element coverage – were evaluated across six machine learning architectures. The allocation of cohorts for the training and validation of the TuFEst model is summarized in Fig. 3A.

SEM emerged as the top-performing classifier, demonstrating a robust discriminative ability across all datasets. In the training cohort (Fig. 3B), both SEM and SVM showed comparable classification accuracies (SEM AUC=0.942, 95% CI: 0.918-0.967; SVM AUC=0.944, 95% CI: 0.918-0.969). This robust performance persisted in internal validation (Fig. 3C), with SEM (AUC=0.937, 95% CI: 0.903-0.970) outperforming other models, including GBM (AUC=0.927) and XGBoost (AUC=0.928). Similarly, SEM maintained superior generalizability in external validation (Fig. 3D), achieving an AUC of 0.968 (95% CI: 0.941-0.995). Based on this demonstrated stability and predictive superiority across diverse datasets, we selected SEM as our final model for the subsequent analyses (Supplementary Table 2). We also evaluated the individual performance of each modality, including DELFI profiles, end-motif signatures, and repeat-element features. The results show that each modality demonstrated meaningful discriminatory power when assessed independently, whereas the integrated TuFEst model achieved the highest overall performance (Supplementary Fig. 1). Confusion matrix analyses (Fig. 3E–G) further validated the diagnostic reliability of SEM. The model demonstrated particularly strong performance in external validation, correctly classifying 91 true-positive BC cases while maintaining low misclassification rates (4 false negatives and 8 false positives). These results translated to clinically robust predictive values [positive predictive value (PPV)=0.919, negative predictive value (NPV)=0.930], meeting the key thresholds for practical cancer screening applications.

Moreover, to evaluate the robustness of TuFEst under ultra-low-coverage conditions, we down-sampled the original 2× LP-WGS (Low-pass Whole-genome Sequencing) data to 1×, 0.5×, and 0.1× (Supplementary Table 3). Across all three datasets—training, internal validation, and external validation—the model maintained highly stable performance over the 1× and 0.5× settings. For example, AUCs changed minimally (ΔAUC≤0.02 in all datasets), and sensitivities remained within 1–3% of the 2× baseline. Performance deterioration became more apparent only at 0.1×, where AUCs decreased modestly by 0.03–0.11 and specificity dropped more substantially, particularly in the external validation cohort. These findings indicate that sequencing depths around 1×–2× provide an optimal balance between cost efficiency and diagnostic accuracy, while confirming that TuFEst retains strong discriminative power even under ultra-low-coverage conditions.

Further investigation of the SEM classifier revealed the key fragmentomics signatures driving its diagnostic performance. Notably, we identified significant enrichment of particular 4bp (e.g., AAGC, AGGC) end motifs and 6bp (e.g., CTTGTT, TGGGAC) breakpoint motifs in BC samples, which showed marked differences between malignant and benign lesions (all P<0.001). In addition, repetitive genomic elements with Alu sequences (e.g., AluSx and AluY), long terminal repeats (e.g., LTR12C, LTR37B), and satellite repeats (e.g., Tigger4b) all showed significant frequency alterations in BC patients compared to benign controls (all P<0.001). Hierarchical clustering analyses demonstrated that these combined fragmentomic features achieved robust separation of malignant and non-malignant samples (Supplementary Fig. 2), strongly supporting their value as biomarkers for cfDNA-based diagnostic applications.

Comparative analyses of TuFEst-derived cancer scores revealed significantly elevated values for malignant cases across all cohorts versus benign controls (training set 9.98-fold; internal validation 4.49-fold; external validation 10.4-fold). Notably, the median cancer scores demonstrated remarkable consistency between cohorts (training set 0.8676; internal validation 0.8671; external validation 0.8659), indicating the robust quantitative performance of the TuFEst algorithm (Fig. 3H–J).

To evaluate TuFEst’s potential in addressing imaging false-negatives, we analyzed blood samples from 26 early-stage BC patients from four hospitals (SAHZU, HZFH, ZJTCM, and LHLI) who had been classified as BI-RADS 3 on both ultrasound and mammography during routine screening. Despite their initially benign imaging assessment, all patients subsequently underwent biopsy because of lesion enlargement within six months, with pathological confirmation of invasive carcinoma. In a single-blind validation using blood samples collected at the time of initial imaging, TuFEst correctly identified 25 of the 26 cases (96.2% sensitivity). Notably, all detected cases represented early-stage cancers (DCIS n=13; pT1a n=1; pT1b n=2; pT1c n=10) without lymph node involvement, demonstrating both the model’s high sensitivity for early-stage BC detection and its potential to minimize false-negative imaging (Fig. 3K, L).

Collectively, these results demonstrate the clinical validity of our cfDNA fragmentomic approach, achieving high diagnostic accuracy, sensitivity, and specificity in distinguishing malignant from benign breast lesions, thereby addressing a critical unmet need in current screening paradigms.

cfDNA fragmentomics predicts BC molecular subtypes
BC treatment decisions are fundamentally guided by molecular subtyping based on hormone receptor and HER2 status, as distinct subtypes demonstrate characteristic biological behaviors and tumor microenvironment profiles25. Although liquid biopsies have shown promise in detecting tumor-associated genetic alterations, the ability of cfDNA fragmentomic patterns to discriminate molecular subtypes remains unexplored. To address this gap, we developed a fragmentomics-based machine learning classifier, TuFEst-MS (TuFEst-Molecular Subtype), designed to non-invasively predict BC molecular subtypes using cfDNA fragmentation profiles. The cohort allocations and subtype distributions used for model training and validation are shown in Fig. 4A.

Patients with BC were classified into three molecular subtypes: ER+/PR+ HER2-, HER2+ (triple-positive and HER2-enriched), and TNBC (ER−PR−HER2−). Given the limited size of TNBC and HER2-enriched subgroups, we combined all available samples and used a unified training-validation split rather than partitioning into separate internal and external validation sets. The TuFEst-MS model demonstrated robust predictive performance across all subtypes, with training set AUCs of 0.906 (95% CI: 0.872-0.940) for ER+/PR+ HER2-, 0.925 (95% CI: 0.890-0.961) for HER2+, and 0.891 (95% CI: 0.820-0.963) for TNBC (Fig. 4B). Remarkably, the performance of the validation set showed excellent consistency, achieving AUCs of 0.939 (95% CI: 0.907–0.970), 0.925 (95% CI: 0.887–0.964), and 0.893 (95% CI: 0.813–0.974) for the respective subtypes (Fig. 4C). This minimal performance variance between the training and validation sets confirms the robustness and generalizability of the TuFEst-MS classifier.

Our fragmentomic analyses revealed that subtype-specific patterns in end motifs, breakpoint motifs, CNV-related features, and repetitive element distributions were the key discriminative features in the TuFEst-MS model (Supplementary Fig. 3). These findings demonstrate that cfDNA fragmentation profiles can mirror the transcriptional and genomic heterogeneity underlying different BC subtypes.

Notably, molecular profiles (e.g., hormone receptor/HER2 status) often differ between primary and metastatic lesions owing to tumor evolution and heterogeneity26. While the current ABC5 (5th ESO-ESMO international consensus guidelines for advanced breast cancer) guidelines recommend treatment strategies based on metastatic lesion profiling, biopsy in accessibility frequently limits clinical implementation27. To address this, we evaluated TuFEst-MS in 21 oligometastatic patients from two centers (SAHZU and SRRZU) (Supplementary Table 4). The model achieved an accuracy of 85.7% (18/21) in predicting metastatic lesion subtypes (Fig. 4D), including an 87.5% accuracy (7/8) in primary-metastatic discordant cases. Three classification discrepancies were observed, all involving HER2 status interpretation: one biopsy-confirmed ER+/PR+/HER2+ (IHC 3+) case was misclassified as HER2-negative by TuFEst-MS; two metastatic lesions confirmed as HER2-negative (IHC 2+/FISH-) were incorrectly predicted as HER2-positive. These discordances mirror the well-documented challenges in HER2 assessment28, where spatial intratumoral heterogeneity and temporal expression dynamics during disease progression frequently yield discordant pathological results based on biopsy samples.

cfDNA fragmentomics for lymph node metastasis detection
The accurate detection of metastatic lymph node (mLN) remains a significant clinical challenge in early BC management. While mLN represents the most common initial metastatic event, conventional imaging modalities (ultrasound, mammography, and MRI) demonstrate suboptimal accuracy29, potentially leading to dual clinical consequences: missing the opportunity for neoadjuvant therapy for node-positive patients while subjecting node-negative cases to unnecessary axillary surgeries. Considering that mLN reflects more aggressive tumor biology, likely to release greater amounts of circulating tumor DNA, we sought to explore cfDNA fragmentomic signatures for LN status assessment – an approach that has merely been investigated to date. In our cohort of 115 mLN BC patients, we developed the TuFEst-LN (TuFEst-Lymph Node) model through comprehensive training (n=250) and validation (n=154) analyses (cohort allocation summarized in Fig. 5A). The model demonstrated strong discriminatory performance, achieving AUCs of 0.900 (95% CI: 0.856-0.945) in the training set and 0.874 (95% CI: 0.816-0.931) in the validation (Fig. 5B, C). Detailed evaluation revealed robust NPV (95.2% in training and 90.3% in validation), although PPV was more moderate (51.5% and 54.9%, respectively), likely due to the relatively low axillary tumor burden in our cohort, as 62.8% of all LN-positive patients were classified as pN1 (Fig. 5D, E).

To further assess clinical utility, we applied TuFEst-LN to two challenging patient groups with discordant imaging and pathological findings that were not included in the training or validation sets. Among cases with imaging-negative but pathology-confirmed LN involvement (Group A), the model correctly identified 60% (24/40) as LN-positive. Conversely, in patients with imaging-positive but pathology-negative results (Group B), 97.6% (83/85) were accurately classified as LN-negative (Fig. 5F). These independent validations confirmed the model’s excellent performance in LN assessment, with particularly high specificity in identifying patients without mLN (NPV=97.6%), which could potentially support clinical decisions to omit axillary surgery in cN0 patients with negative TuFEst-LN predictions.

TuFEst cancer scores reflect BC molecular subtype aggressiveness
Although stringent standardization of sample collection and processing protocols was implemented across all participating centers, we observed notable inter-sample variability in cancer scores. Given the well-established differences in clinical aggressiveness among BC molecular subtypes, we performed a stratified analysis comparing cancer score distributions across the ER+/PR+ HER2-, HER2-positive, and TNBC groups. Using multiple threshold definitions for high (top 10–30%) and low (bottom 10-30%) cancer score groups, we consistently found that HER2-positive and TNBC subtypes were significantly enriched in high cancer score populations across all cut-off values (Supplementary Fig. 4A). In contrast, low cancer score distributions remained relatively balanced across molecular subtypes (Supplementary Fig. 4B). These findings corroborate the known biological behavior of TNBC and HER2-positive tumors, which demonstrate greater metastatic potential even at early disease stages compared to the ER+/PR+ HER2- subtype, and validate the TuFEst model’s ability to identify high-risk molecular subtypes at early detection phases, potentially facilitating timely interventions to improve patient prognosis.

Transcriptomic profiling links TuFEst-derived cancer scores to BC biology
Because cfDNA fragmentation is closely linked to chromatin accessibility and transcriptional activity, transcriptomic profiling provides a direct mean to interpret the biological programs underlying TuFEst-derived Cancer Scores and further enables inference of the tumor immune microenvironment through an established deconvolution framework. Thus, to molecularly characterize early-stage BC detected by liquid biopsy, we performed bulk RNA sequencing on 79 matched primary tumors from patients whose blood samples were used to develop the TuFEst model. Focusing on the top and bottom 20% of cases based on cancer scores (n=15 per group), we found no significant differences in clinicopathological features between the high- and low-score groups (Supplementary Table 5), ensuring comparability for further analyses.

Transcriptomic profiling revealed distinct biological patterns between the groups. Both the xCell and ESTIMATE algorithms—widely used and well-validated tools for inferring tumor microenvironment composition — consistently showed significantly higher immune infiltration (xCell: P=0.0453; ESTIMATE: P=0.0408) but comparable stromal scores in high-score tumors (Fig. 6A). CIBERSORT analyses identified plasma cells (P=0.0404) and resting mast cells (P=0.0179) as being enriched in low-score tumors – a cellular profile characteristic of immunologically “cold” microenvironments30. In contrast, high-score tumors exhibited increased infiltration of cytotoxic immune populations (T cells, NK cells) and M0/M1 macrophages, suggesting enhanced immune activation even at early disease stages and potentially indicating favorable responses to immunotherapy regimens (Fig. 6B).

Tumors from the high-score group exhibited significant enrichment across multiple malignancy-associated hallmark pathways31, revealing a coordinated activation of proliferative programs (E2F Targets, G2M Checkpoint, MYC Targets V1/V2), microenvironment remodeling processes (Angiogenesis, EMT), metabolic reprogramming (glycolysis), and diverse immune-inflammatory networks (Inflammatory Response, TNF-α Signaling, IFN-γ Response) (Fig. 6C). The top differentially expressed genes have been functionally implicated in fostering a pro-tumorigenic microenvironment by inducing chronic inflammation (S100A8, S100A9, CXCL8), immune suppression (BCL2A1, OSM, TNFAIP3), and tumor cell survival (BCL2A1, S100A8, S100A9)32–37 (Fig. 6D). These findings were strongly corroborated by gene set enrichment analyses, where EMT (NES=2.647, FDR<0.001), Inflammatory Response (NES=2.800, FDR<0.001), and TNF-α/NF-κB signaling (NES=3.026, FDR<0.001) emerged as the most significantly activated pathways in high-scoring tumors (Fig. 6E). Collectively, our bioinformatic analyses support the notion that our cfDNA fragmentomic pattern-based model effectively captures the critical biological characteristics of BC, especially those linked to great metastatic potential and immune microenvironment remodeling.

Methods
Ethical statements
All procedures involving human participants complied with all relevant ethical regulations and were conducted in accordance with the Declaration of Helsinki. All datasets were collected under institutional review board approval (Approval No. 2023-1048 by the Second Affiliated Hospital of Zhejiang University School of Medicine, Approval No. IIT-20230612-0118-01 by Hangzhou First People’s Hospital, Approval No. 2024-KLS-013-02 by Zhejiang Provincial Hospital of Traditional Chinese Medicine, Approval No. 2024-033 by Ningbo Medical Center Lihuili Hospital, Approval No. IRB-20230357-R by Women’s Hospital School of Medicine, Zhejiang University, Approval No. 2023-0548 by Sir Run Run Shaw Hospital, Approval No. 2024-002 by XiangFu Community Health Service Center). All data were de-identified before model development. Written informed consent was obtained from all participants prior to study participation.

Study population
This study was designed as a prospective multicenter case-control study (NCT06016790) that continuously enrolled 839 participants, including 550 BC patients and 289 benign controls from seven Chinese medical institutions: the Second Affiliated Hospital of Zhejiang University School of Medicine (SAHZU), Hangzhou First People’s Hospital (HZFH), Zhejiang Provincial Hospital of Traditional Chinese Medicine (ZJTCM), Ningbo Medical Center Lihuili Hospital (LHLI), Women’s Hospital School of Medicine Zhejiang University (WHZJU), Sir Run Run Shaw Hospital (SRRSH), and Xiangfu Community Health Service Center (XFCHC). Eligible participants were females aged 18-70 years with pathologically confirmed diagnoses of either: (1) primary invasive breast carcinoma or carcinoma in situ, or (2) benign breast lesions (fibroadenoma/adenosis/fibrocystic hyperplasia). All participants provided their plasma samples. Exclusion criteria included pregnancy/lactation, prior malignancies (5-year window), suspected but unconfirmed malignancies within the preceding year, recent blood transfusions (30-day window), active mastitis, and other investigator-determined contraindications.

Sample collection and processing
Peripheral blood (8–10 mL) was collected in specialized cfDNA tubes (Ardent BioMed Cat. # BY10240301) during routine blood tests or before surgery. Samples were subjected to dual centrifugation (1,600× g → 16,000× g, 10 min each at 4°C) for plasma isolation, followed by immediate storage at −80°C. Frozen plasma was batch-shipped on dry ice to the central laboratory (Omix Science Research Institute, Hangzhou, China) for processing.

cfDNA assays and whole-genome sequencing
cfDNA was extracted from a median plasma volume of 1 mL using the VAMNE Mag Ultra circulating cell-free DNA isolation kit (Vazyme Biotech, Nanjing, China, Cat. # N913) according to the manufacturer’s protocol, followed by precise quantification using a Qubit 4.0 fluorometer (ThermoFisher Scientific, Lenexa, KS, USA) and storage at −80°C. For library preparation, 5–20 ng of the extracted cfDNA was processed with the VAHTS Universal DNA library prep kit for Illumina V3 (Vazyme Biotech, Nanjing, China, Cat. # N610) and barcoded using VAHTS multiplex oligos set 5 for Illumina (Vazyme Biotech, Nanjing, China, Cat. # N322). The quality and integrity of the libraries were evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) to confirm that the libraries met the necessary standards for optimal fragment size distribution. The final libraries were sequenced on the NovaSeq X Plus platform (Illumina, San Diego, CA, USA) using paired-end 150bp reads, with an average coverage depth of 2× for the cfDNA samples.

Raw WGS data underwent rigorous quality control beginning with FastQC (version 0.12.1, https://www.bioinformatics.babraham.ac.uk/projects/fastqc) evaluation, followed by quality filtering and adapter trimming using Cutadapt (version 4.5, https://github.com/marcelm/cutadapt) and Ktrim (version 1.4.1, https://github.com/hellosunking/Ktrim). The processed reads were aligned to the human reference genome (GRCh37) using BWA-MEM (version 0.7.17, https://github.com/lh3/bwa) with default settings, ensuring that only high-quality reads with a Phred score above the threshold were retained. Further analyses were restricted to uniquely mapped reads that were free of PCR duplicates. These mapped reads were subsequently sorted and indexed using Samtools (version 1.9, https://github.com/samtools/samtools), and duplicates were identified and eliminated using Picard (version 2.18.29, https://github.com/broadinstitute/picard). The distribution of fragment sizes and relevant genomic features was inferred based on the coordinates of mapped read pairs, following established protocols using Picard. To identify large-scale epigenetic variations in cfDNA fragmentation across the genome, which can be detected with low-coverage WGS, the hg19 reference genome was segmented into non-overlapping 5 Mb bins.

Fragmentomic feature identification
We systematically quantified the frequencies of 4-bp end motifs and 6-bp breakpoint motifs, normalizing total motif occurrences to unity within each category to enable comparative analyses. For the 6 bp motifs, we characterized the 3 bp sequences flanking cfDNA 5’ breakpoints (both upstream and downstream) and established genome-wide distribution profiles. To maintain consistency, the overall frequency of all 6 bp breakpoint motifs was normalized to one. We further organized the motif features into an I×J matrix, where each row corresponds to a participant identified by their sample ID. This matrix format allows for efficient interpretation and comparison of motif patterns within the cohort, leading to deeper insights into cfDNA motif distributions.

Regional fragment size analyses employed 100 kb genomic bins using Deeptools (version 3.1.2, https://github.com/deeptools/deepTools), optimized to balance resolution with sufficient read depth (~25,000 reads/bin). Each bin contained short (S, 100–150 bp) and long (L, 151-220 bp) fragments. To ensure data accuracy, regions overlapping the ENCODE blacklist or the hg19 gap track from the UCSC Genome Browser were excluded, as these regions typically show poor mappability. This approach yielded robust S/L ratios for the fragmentation pattern analyses. By focusing on fewer, higher-quality bins, we improved the reliability of our fragment-size distribution analyses.

Alu repeats analyses of cfDNA
We categorized all elements from the gene and pseudogene families as transfer RNA (tRNA), signal recognition particle RNA (srpRNA), small nuclear RNA (snRNA), small cytoplasmic RNA (scRNA), and ribosomal RNA (rRNA). Likewise, we grouped elements from the DNA, retroposon, and Rolling Circle (RC) families into the transposable element category. This led to eight primary groups for analysis: LINE, SINE, Medium Repeat (MER), DNA_TcMar-Tigger (DNA Transposable elements), Long Terminal Repeat (LTR), satellites, transposable elements, and RNA elements. In our cfDNA feature analyses and validation, we focused on critical elements, including Alu sequences, long terminal repeats (LTRs), RNA elements, and transposable elements. Through this systematic classification, we adopted a thorough approach to understand cfDNA features, which are essential for investigating their roles in BC. The TuFEst algorithm integrates these repeat elements with the N-index and E-index and incorporates TF features to predict non-cancerous and cancerous status.

Machine learning and deep learning algorithm model
In this study, TuFEst was implemented as a binary classification model50. It outputs a “cancer score” between 0 and 1, representing the probability that a plasma cfDNA sample originates from an individual with cancer. A cutoff of 0.5 was applied to this score to classify samples as “cancer” or “non-cancer.” Although the score is a surrogate measure of the relative burden of tumor-derived DNA and the name suggests its long-term potential for quantification, the current application and validation focus on classification performance.

A GLM with a binomial distribution and a logit link function, equivalent to logistic regression, was developed for binary classification. The model utilized fragmentomics features as predictors, with the outcome variable coded as 1 for malignant and 0 for benign status. To address potential sample dependencies, group information was incorporated using Sample ID. The analytical pipeline included feature standardization (StandardScaler) followed by logistic regression with a maximum of 1000 iterations. Model performance was assessed through 10-fold stratified cross-validation with shuffling (random state = 42) to ensure reproducibility while maintaining class balance across folds.

An SVM with a linear kernel was implemented for classification. The model employed standardized fragmentomics features, with malignant status encoded as 1 and benign as 0. Sample identifiers were retained for group-level analysis. The pipeline consisted of feature scaling followed by SVM classification with probability estimation enabled. Performance was evaluated using 10-fold stratified cross-validation, with metrics including accuracy, AUC, sensitivity, specificity, PPV, and NPV. Probability thresholds at specificity levels of 0.9 and 0.95 were derived from ROC analyses to support clinical decision-making. All predictions and probabilities were systematically recorded for aggregated analysis.

An XGBoost classifier was implemented using a reproducible pipeline (random state = 42). Fragmentomics features were standardized prior to model training, with binary outcomes indicating malignant (1) or benign (0) status. The model was evaluated using 10-fold stratified cross-validation, with comprehensive metrics recorded for each fold. Clinical decision thresholds were established by identifying probability values corresponding to 0.9 and 0.95 specificity levels from ROC curves. All predictions and sample metadata were stored for subsequent analysis, and the final pipeline was saved for potential future use.

An MLP neural network was developed with fixed random state (42) for reproducibility. Features were standardized before being input to the MLP classifier. The model was evaluated using 10-fold stratified cross-validation, with performance metrics including accuracy, AUC, sensitivity, specificity, PPV, and NPV. Probability thresholds for specific specificity levels (0.9 and 0.95) were determined from ROC analyses to define clinical decision boundaries. The trained pipeline, including preprocessing and classification steps, was preserved for potential application, and overall discriminative ability was validated through aggregated ROC analysis.

A GBM was implemented using scikit-learn’s GradientBoostingClassifier. The model was trained on standardized fragmentomics features and evaluated via 10-fold stratified cross-validation (random state = 42). Each fold yielded probability estimates and performance metrics, including accuracy, AUC, sensitivity, specificity, PPV, and NPV. Clinical thresholds were identified by setting specificity levels at 0.9 and 0.95 in ROC analysis. Predictions, true labels, and group identifiers were stored for combined analysis, and the final model was saved for potential deployment.

A SEM integrated predictions from five base learners - SVM, GLM, GBM, XGBoost, and MLP - using a logistic regression meta-learner. The pipeline included feature standardization and was evaluated through 10-fold stratified cross-validation. Performance metrics were calculated for each fold, and probability thresholds were established at specificity levels of 0.9 and 0.95. The ensemble demonstrated strong discriminative ability (AUC=0.943), with all predictions and sample information systematically recorded. The complete pipeline was saved to support potential clinical implementation.

For tumor detection and classification, we screened six distinct machine learning algorithms to build the models: GLM, SVM, XGBoost, MLP, GBM, and SEM. These models were rigorously trained and evaluated using bootstrapping combined with 10-fold cross validation. During the 10-fold cross-validation process, the training set was randomly divided into 10 subsets, with nine subsets used for training and the remaining one used for testing, cycling 10 times to ensure that each subset was used for testing once. We only used training cohort samples from benign individuals and cancer patients, training classifiers in the machine learning algorithms with features such as motif feature frequencies, and generating models to predict cancer scores for each sample. Notably, all validation datasets remained untouched during model training. The cancer score ranged from 0 to 1, with a higher score indicating a greater probability of cancer. After the evaluation, we selected the best-performing model for the downstream analyses.

These models were subsequently applied to the validation datasets to generate cancer prediction scores for each validation sample and assess the model performance. We compared the AUC values of the different models in the validation cohort and the sensitivity/specificity at a fixed specificity threshold in the internal validation cohort.

Three specialized models have been developed for different clinical applications. The TuFEst model was designed to distinguish between benign individuals and cancer patients using a stacked ensemble that incorporates feature information, including TF binding site coverage, to enhance predictive accuracy. The TuFEst-MS model employs a GLM method to classify three BC subtypes: ER+/PR+, HER2+, and TNBC. This classification capability is particularly important for predicting the subtype of metastatic lesions in patients with BC. The TuFEst-LN model utilized GLM to predict whether a patient had mLN, specifically for analyzing patients with radiological evidence of mLN but no surgical confirmation, or those with discordant imaging-pathological results. All models underwent rigorous training using 10-fold cross validation, with prediction scores calculated for each sample. The models were assessed using the AUC-ROC metric to ensure their accuracy and reliability in practical applications.

Selection and analyses of differential features
We utilized the wilcox.test() function to perform the Wilcoxon rank-sum test, assessing the statistical significance of differences in various features across different classifications. The parameter exact = FALSE indicates the use of an approximate computational method. Features with a p-value less than 0.05 were selected to ensure statistical significance and eliminate the influence of random factors. Additionally, we applied a fold change threshold greater than 1, retaining only features that exhibited significant differences between the two groups. Subsequently, the top and bottom 30 features were identified, and the median ratio between the benign and malignant groups was calculated for each feature. Finally, a standardized feature heatmap was generated using the pheatmap() function from the pheatmap package, while the scale() function was applied for Z-score normalization, enhancing visualization and facilitating further analyses.

Genome-wide TF analyses
ChIP-seq peaks from 5620 experiments were downloaded from the ReMap 2020 database (https://remap2020.univ-amu.fr/download_page), selecting the CRMs (Cis-Regulatory Modules) datatype, which includes a total of 1,732,560 peaks.

For each autosomal peak, the peak center was designated as position 0, and the mean coverage was computed for each sample within a ±1,000 bp window relative to the peak center. The processed data were then incorporated as a feature into the model training pipeline to distinguish between benign individuals and patients with cancer. A series of filters were applied to ensure sufficient sequencing depth across benign and cancer samples, retaining only peaks with read counts > 25, read coverage ≥ 0.9, and fold change ≤ 0.5 or ≥ 2.

To investigate the biological processes associated with the identified TFs, we performed KEGG pathway enrichment analyses using the clusterProfiler package in R. TFs were first mapped to Entrez IDs, and the KEGG enrichment analysis was applied to identify significantly enriched pathways. Pathways were ranked based on their adjusted p-values, and the enriched pathways were visualized in a bubble chart using ggplot2, with bubble size representing gene count and color indicating p-value significance.

Estimation of the tumor immune microenvironment
To comprehensively profile the tumor immune microenvironment, we employed two computational deconvolution methods: CIBERSORT51 and xCell52, which were used to quantify the relative proportions of distinct tumor-infiltrating immune cell populations from the bulk transcriptomic data. In addition, the ESTIMATE algorithm53 was applied to calculate immune and stromal scores, reflecting the overall infiltration levels of immune and stromal components within the tumor tissue.

Gene enrichment and pathway analyses
Gene Set Enrichment Analysis (GSEA) was conducted using the “clusterProfiler” R package54 to identify significantly enriched hallmark gene sets. GSVA was performed using the “GSVA” R package55 to compute sample-wise enrichment scores for specific gene sets. Differentially expressed genes between the Cancer Score-High and Cancer Score-Low groups were identified using the “limma” R package, with thresholds set at |log₂(fold change)| > 1 and P < 0.05.

Statistics & reproducibility
For statistical analysis, a ROC curve was generated using the pROC package (v. 1.18.5). Based on the true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN) for cancer prediction, we calculated the sensitivity [TP/(TP+FN)], specificity [TN/(Tn+FP)], positive predictive value (PPV) [TP/(TP+FP)], negative predictive value (NPV) [TN/(Tn+FN)], and accuracy [(TP+TN)/(TP+FP+Tn+FN)]. These metrics were computed using the confusion matrix function of the caret package (v. 6.0.94) in R (v. 4.4.1). The preProcess function (v. 1.2.2) from the caret package was used for feature standardization. Additionally, the confusion matrix algorithm was employed to compare ROC curves, providing a comprehensive evaluation of the diagnostic test performance and ensuring accurate differentiation between benign and malignant cases. No statistical method was used to predetermine sample size. No data were excluded from the analyses. Sample processing for cfDNA library preparation and sequencing was randomized by case/control status to minimize batch effects. The Investigators were not blinded to allocation during the experiments or during outcome assessment.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Zhu et al. describes classification models based on cfDNA fragmentomics for the early detection of breast cancer, molecular subtyping, and identification of disease metastasis, which have potential for improving the clinical management and treatment of breast cancer patients. Generally, the study is well-designed, supported by a large sample size, and includes valuable validation samples, such as blood specimens from patients with benign lesions confirmed by imaging diagnosis. The results appear to be robust and interpretable. However, the manuscript lacks sufficient background explanation regarding the rationale of the study design, and the comparisons with existing competitive methods are not adequately addressed, which should be strengthened to improve the overall clarity and contextual significance of the work.\nPoint 1\nThe authors report that their model achieved high performance, with an AUC > 0.93 for breast cancer detection, demonstrating a sensitivity of 75% at a specificity of approximately 85%. Presenting a table with multiple sensitivity and specificity pairs would enhance clarity. Furthermore, discussing which of these two metrics—sensitivity or specificity—is more critical in the context of real-world breast cancer detection and clinical decision-making would be helpful for clinicians to intuitively understand the value of this technique.\nPoint 2\nAlthough the Discussion section (from line 455 onward) provides relevant interpretations, the Introduction generally lacks sufficient explanation regarding the overall study design and comparisons with other technologies. For example, the use of transcription factor (TF) motifs and repeat elements in the model should be clearly stated in the Introduction, along with appropriate background information. In addition, some references (e.g., references #4 and #5) do not appear to be the most appropriate for supporting the claims made. Given the growing number of liquid biopsy technologies for cancer detection, a more thorough comparison with existing approaches should be presented in the Introduction. Finally, the performance of the model in breast cancer detection and subtyping would be more informative if compared directly with that of other existing methods using a summary table (e.g. https://www.frontiersin.org/journals/oncology/articles/10.3389/fonc.2023.1127086/full,\nhttps://www.nature.com/articles/s41523-021-00316-7)\nPoint 3\nAt line 194, the authors describe filtering transcription factors (TFs) enriched in breast cancer using a breast cancer cell line. It would be helpful to include a discussion on how incorporating TFs enriched in benign conditions might influence model performance.\nPoint 4\nIn Figure 3K, the model correctly identified 25 out of 26 samples. Given the importance of early-stage detection—both in this dataset and in the broader clinical context—it would be valuable to present the model’s performance specifically on early-stage samples, or stratified by stage. Also, while the 25/26 result in Figure 3K suggests outstanding early detection performance, a discussion comparing this result to the model’s performance on the broader dataset would provide important context and help assess the generalizability of the approach.\nPoint 5\nRegarding Figure 3B, it would be helpful to present the individual performance of each modality (DELFI, end-motif, and repeat elements) at least in the Supplementary Materials. This would allow readers to better assess the relative significance of each modality and their respective contributions to the performance of the combined model.\nPoint 6\nThe authors highlight the advantage of their approach as being suitable for ultra-low coverage settings. To support this claim, it would be helpful to present the coverage distribution across samples. As the authors note, lower coverage is indeed advantageous in clinical situations due to constraints in blood volume and cost. Therefore, it would strengthen the manuscript to include an analysis of breast cancer detection performance across different coverage levels or by down-sampling samples to low coverages (0.1~1X). If feasible, similar analyses for the other models presented in the manuscript would also be valuable. Additionally, a discussion explaining why low-pass sequencing seems to be required.\nPoint 7\nThe authors claim that alternative technologies require high costs. To support this point, it would be helpful to include a comparison of the authors method to other existing methods regarding sequencing or assay costs.\nPoint 8\nFigure 2A appears to show less centralized plots across samples compared to DELFI. Even in the healthy group, the patterns seem less centralized compared to the DELFI healthy reference. It would be helpful for the authors to provide a rationale or explanation for this observation.\nPont 9\nAt lines 187–189, the manuscript suddenly shifts to discussing transcription factors (TFs), which may be confusing readers with the rationale for incorporating TF features. Providing additional background or a brief explanation would improve readers understanding. Similarly, at line 398, it is hard to get reasons why RNA-seq was chosen for downstream analysis, especially given the availability of other molecular data types such as methylation, SNVs, and CNVs. Although the rationale becomes somewhat clearer later in the Discussion (starting at line 478), the authors should consider introducing the reasoning earlier in the manuscript. A brief explanation of why xCell and ESTIMATE were selected—and how they support the study objectives—would also be helpful.\nPoint 10\nIn the manuscript showed the detection of lymph node metastasis. It would be helpful for the authors to clarify whether lymph node metastasis could be detected within a unified model along with other types of metastasis (Figure 4). If this is not feasible, a discussion explaining the limitations or reasons for separating these models would strengthen the manuscript.\nPoint 11\nIn the Discussion section (line 436), the authors mention the limitations of bisulfite-based methods. However, more recent approaches such as enzymatic C-to-T conversion and Oxford Nanopore Technologies (ONT)-based methylation profiling offer viable alternatives that may mitigate such limitations. The authors should provide more precise, performance-based evidence to support their claims. Furthermore, performance comparisons with existing studies—ideally in tabular format—would enhance the manuscript. In particular, at line 464, the authors compare their model to studies in other cancer types. It would be more clear to include direct comparisons with recent liquid biopsy studies of breast cancer.\nPoint 12\nAt line 438, the authors refer to their approach as “AI-based.” However, given the methodological context, it may be more appropriate to describe it as “machine learning–based,” as the use of “AI” could be considered too advanced or potentially misleading in this case.\nPoint 13\nAt line 457, ㅅhe authors mention using a tumor score, but what is the corresponding limit of detection in terms of variant allele frequency (VAF)? A brief discussion on this point would clarify the model’s sensitivity.", "reviewer_2": "The manuscript introduces a new multi-feature cfDNA-based Breast cancer detection method. The method is trained and tested on a cohort of 503 cases and 289 controls and shows good classification performance. Additionally, they train a ML model that can perform molecular subtyping using the same features. The results are pretty good and of interest to people in the field, but the methods section is lacking details about the features and how the method was trained.\nYou use three different kinds of fragmentomics features. How many features do you end up with of each kind?\nHave you tested which kind of fragmentomics features contribute the most to the performance of the method?\nDescribing the motif analysis, you write that the 6-bp motifs are centered around the breakpoint. Is this also the case for the 4bp motifs, or are they only the end-motif?\nIf you count motifs in both ends of the fragments, do you use the reverse complement of one of the ends so that the motifs are comparable?\nThere are no details about the different ML methods that are trained provided in the manuscript. I looked at the linked GitHub repo, but as far as I can see, that repo only contains code for plotting. You should both add more details about the parameters and implementations of the methods to the manuscript and add the code to GitHub to make the details regarding the training and testing of the models available.\nFor instance, is it a linear SVM you test, or do you use a kernel? Do you use cross-validation to fit any meta-parameters for any of the ML methods, or do you use default parameters? You end up using a Stacked Ensemble ( I would write out the abbreviation the first time you use it), but that is just a method for combining other ML models, and you do not say anything about which models you are combining (is it the 5 other mentioned models?).\nThe method is called TuFEst (Tumor Fraction Estimator), but the manuscript only reports on classification performance; there are no plots showing an estimated tumor fraction compared to some independent estimate of the tumor fraction in the cfDNA samples.\n(Remarks on code availability)\nAs far as I can see, that repo only contains code for plotting. I could not find code for training the different models mentioned in the manuscript.", "reviewer_4": "This is a well written manuscript demonstrating the power of cfDNA fragments to detect breast cancer as well as determine subtype and lymph node status.\nIt was a well conducted study using both benign and malignant breast tissue.\nI think that this study will be of interest to researchers and clinicians involved in early breast cancer detections and treatment.\nI have very minor concerns:\n1. I could not locate Table 1 cited on line 209.\n2. there is no discussion of tumors which are triple positive, whether they were included or excluded in the analyses."}, "author_response": "Point-by-point response\nReviewer #1 (Remarks to the Author):\nThe manuscript by Zhu et al. presents classification models utilizing cfDNA fragmentomics for the early detection of breast cancer, molecular subtyping, and identification of metastasis disease, showing promise for enhancing clinical management and treatment of breast cancer patients. Overall, the study is well-designed, with a large sample size and valuable validation cohorts—including blood samples from patients with imaging-confirmed benign lesions. The findings appear robust and interpretable.\nHowever, the manuscript lacks sufficient background explanation regarding the rationale of the study design, and the comparisons with existing competitive methods are not adequately addressed, which should be strengthened to improve the overall clarity and contextual significance of the work.\nPoint 1\nThe authors report that their model achieved high performance, with an AUC > 0.93 for breast cancer detection, demonstrating a sensitivity of 75% at a specificity of approximately 85%. Presenting a table with multiple sensitivity and specificity pairs would enhance clarity. Furthermore, discussing which of these two metrics—sensitivity or specificity—is more critical in the context of real-world breast cancer detection and clinical decision-making would be helpful for clinicians to intuitively understand the value of this technique.\nResponse：We sincerely thank the reviewer for this constructive suggestion. In response, we have added a dedicated table summarizing multiple sensitivity–specificity operating points for TuFEst across both the training and validation cohorts, derived from different machine learning models. This table has been included as Supplementary Table 2 in the revised manuscript.\nRegarding the relative importance of sensitivity versus specificity in real-world breast cancer early detection, we have expanded the discussion as follows (Lines 484–494 in the revised version):\n“In real-world early detection, the optimal balance between sensitivity and specificity must be tailored to the clinical context. For population-level screening, particularly among Asian populations with dense breast tissue and suboptimal mammographic performance, higher sensitivity is essential to minimize the risk of missed early-stage cancers. Since our cfDNA fragmentomics approach is intended to complement routine ultrasound at the initial diagnostic stage, adopting a higher-sensitivity threshold supports timely referral to MRI or biopsy and reduces the likelihood of false-negative outcomes. At the same time, maintaining reasonable specificity is crucial to avoid unnecessary imaging or invasive procedures. In this regard, the TuFEst model achieves a clinically practical balance (specificity: 0.783 in training database; 0.869 in external validation database), enabling flexible integration across diverse healthcare settings and diagnostic workflows.”\nPoint 2\nAlthough the Discussion section (from line 455 onward) provides relevant interpretations, the Introduction generally lacks sufficient explanation regarding the overall study design and comparisons with other technologies. For example, the use of transcription factor (TF) motifs and repeat elements in the model should be clearly stated in the Introduction, along with appropriate background information. In addition, some references (e.g., references #4 and #5) do not appear to be the most appropriate for supporting the claims made. Given the growing number of liquid biopsy technologies for cancer detection, a more thorough comparison with existing approaches should be presented in the Introduction. Finally, the performance of the model in breast cancer detection and subtyping would be more informative if compared directly with that of other existing methods using a summary table (e.g.\nhttps://www.frontiersin.org/journals/oncology/articles/10.3389/fonc.2023.1127086/full,\nhttps://www.nature.com/articles/s41523-021-00316-7)\nResponse: We thank the reviewer for these excellent suggestions. Firstly, we expanded the Introduction to briefly outline our overall study design and to clarify the rationale for incorporating multiple fragmentomic modalities, as described in Lines 144–150 of the revised manuscript. Specifically, we added concise background on transcription factor–associated fragmentation and repeat-element features to justify their inclusion and to improve the conceptual flow of the Introduction.\nSecondly, we carefully reviewed and updated the references throughout the Introduction. The original References #4 and #5 were replaced with three more appropriate sources (Lines 131–135)—de Abreu et al. (The Oncologist, PMID: 40680231), Liu et al. (Ann Oncol, PMID: 32165213), and Jamshidi et al. (Cancer Cell, PMID: 36372334)—as these studies more accurately support the limitations of mutation-based assays (low VAF requiring deep sequencing) and methylation-based assays (modest Stage-I sensitivity and substantial computational burden).\nThirdly, we added a comprehensive comparison of breast cancer early-detection approaches in the revised manuscript. Supplementary Table 6 now summarizes major diagnostic modalities—including mammography and ultrasound; mutation-based cfDNA assays; bisulfite-based, enzymatic, and ONT methylation assays; and leading fragmentomic platforms such as DELFI, CRAG, SPOT-MAS, cfFrag, GALYFRE, and TAPS—detailing their blood-volume requirements, workflow complexity, assay cost, and diagnostic performance. This table provides a breast-cancer-focused benchmark for evaluating the performance and practical advantages of TuFEst.\nWe believe our revisions substantially enhance the introductory framing and provide readers with a clearer understanding of the motivation, innovation, and clinical value of our work.\nPoint 3\nAt line 194, the authors describe filtering transcription factors (TFs) enriched in breast cancer using a breast cancer cell line. It would be helpful to include a discussion on how incorporating TFs enriched in benign conditions might influence model performance.\nResponse: We thank the reviewer for raising this important point. The TFs analyzed in our study are specifically associated with breast cancer malignancy. Our selection criteria were intentionally designed to identify TFs linked to malignant biological processes, as the primary objective of the model is to distinguish malignant from benign states (PMID: 36399356). As a result, the malignant-associated and benign-associated TF groups are mutually exclusive by definition.\nIn addition, we conducted a dedicated ablation study comparing the full TuFEst model with a “Non-TF” model in which all TF-related coverage features were removed. The results, presented in Rebuttal Figure 1, demonstrate that incorporating TF features substantially enhances discriminative performance across the training, internal validation, and external validation cohorts.\nRebuttal Figure 1. Ablation analysis evaluating the contribution of TF-associated fragmentomic features to TuFEst classification performance.\nROC curves evaluating the discriminative performance of SEM with (A) or without (B) TF-related coverage features across the training, internal validation, and external validation databases.\nPoint 4\nIn Figure 3K, the model correctly identified 25 out of 26 samples. Given the importance of early-stage detection—both in this dataset and in the broader clinical context—it would be valuable to present the model’s performance specifically on early-stage samples, or stratified by stage. Also, while the 25/26 result in Figure 3K suggests outstanding early detection performance, a discussion comparing this result to the model’s performance on the broader dataset would provide important context and help assess the generalizability of the approach.\nResponse: We appreciate the reviewer’s emphasis on early-stage detection. We have performed an expanded stage-stratified analysis of TuFEst performance across Stage 0, I, IIa, and IIb–III samples in the external validation dataset. As shown in Rebuttal Figure 2, the model achieved consistently high detection accuracy across all early-stage categories: 96% in Stage I (46/48 correctly identified), 100% in Stage IIa (21/21). These stage-specific results place the 25/26 detection shown in Figure 3K into broader context and demonstrate that the strong performance observed in that subset is not an isolated finding but is consistent with robust detection across the full staging spectrum.\nRebuttal Figure 2. Stage-stratified performance of the TuFEst model in the external validation cohort.\nPoint 5\nRegarding Figure 3B, it would be helpful to present the individual performance of each modality (DELFI, end-motif, and repeat elements) at least in the Supplementary Materials. This would allow readers to better assess the relative significance of each modality and their respective contributions to the performance of the combined model.\nResponse: We thank the reviewer for this insightful suggestion. Clarifying the relative contributions of individual fragmentomic feature classes is indeed essential for both biological interpretability and future assay optimization. To address this point, we conducted additional analyses to systematically evaluate the predictive performance of each feature category in isolation. Specifically, we assessed models constructed using only DELFI fragmentation profiles, end-motif signatures, or repeat-element features.\nAs shown in the newly added Supplementary Figure 1, each modality exhibited meaningful discriminatory power when evaluated on its own, and we now provide these results alongside the combined model shown in Figure 3B-D (Lines 264-268). Notably, the integrated TuFEst model consistently outperformed all single-modality models across cohorts.\nModel interpretability analysis using SHapley Additive exPlanations (SHAP) further revealed that end-motif features contributed the largest mean absolute SHAP values, followed by repeat elements, TF signals, and DELFI profiles (Rebuttal Figure 3). These results indicate that all feature classes provide complementary and nonredundant biological information, and their integration is responsible for the superior performance of the combined model.\nRebuttal Figure 3. SHAP-based interpretability analysis of different fragmentomic features.\nPoint 6\nThe authors highlight the advantage of their approach as being suitable for ultra-low coverage settings. To support this claim, it would be helpful to present the coverage distribution across samples. As the authors note, lower coverage is indeed advantageous in clinical situations due to constraints in blood volume and cost. Therefore, it would strengthen the manuscript to include an analysis of breast cancer detection performance across different coverage levels or by down-sampling samples to low coverages (0.1~1X). If feasible, similar analyses for the other models presented in the manuscript would also be valuable. Additionally, a discussion explaining why low-pass sequencing seems to be required.\nResponse: We thank the reviewer for this valuable suggestion. We have now incorporated the coverage distribution for all samples and performed the requested down-sampling analysis from 2× to 1×, 0.5×, and 0.1×. As shown in Rebuttal Table 1, Supplementary Table 3, and Rebuttal Figure 4, TuFEst exhibited highly stable performance at 1× and 0.5×, with only minimal changes in AUC (ΔAUC ≤ 0.02) and sensitivity (≤3%). A marked decline in performance was observed only at 0.1× coverage, indicating that sequencing depths in the range of 1–2× offer an optimal balance between assay accuracy and clinical feasibility. These results have been incorporated into the revised Results section (Lines 273–283).\nWe have also added a brief explanation in the Discussion regarding the rationale for low-pass sequencing in fragmentomics as follows (Lines 516–521): “Low-pass sequencing is inherently well suited for fragmentomics because its biological signal arises from genome-wide coverage variability, nucleosome-positioning patterns, and sequence-context features, rather than from high-depth variant detection such as single-nucleotide variants or variant allele frequencies. Importantly, increasing sequencing depth yields diminishing returns, as fragmentomic features tend to plateau once sufficient coverage uniformity has been achieved.”\nRebuttal Figure 4. Robustness of TuFEst performance across reduced sequencing depths.\nROC curves evaluating the diagnostic performance when reducing LP-WGS coverage from 2× (A) to 1× (B), 0.5× (C), and 0.1× (D) across the training, internal validation, and external validation databases.\nPoint 7\nThe authors claim that alternative technologies require high costs. To support this point, it would be helpful to include a comparison of the authors method to other existing methods regarding sequencing or assay costs.\nResponse：We thank the reviewer for this constructive suggestion. In the revised manuscript, we have added a comprehensive comparison of sequencing and assay costs across major early breast cancer detection technologies, including mammography, ultrasound, mutation-based cfDNA assays, methylation-based assays (bisulfite, enzymatic, and nanopore-based), and fragmentomic approaches. These data are summarized in the newly added Supplementary Table 6, which outlines differences in required blood volume, workflow complexity, and estimated per-sample cost.\nAs highlighted in the Results and Discussion, TuFEst operates at the lowest sequencing depth (approximately 2×) and requires minimal sample input, resulting in substantially lower overall assay cost compared to other liquid biopsy methods. These advantages support the feasibility of our approach for large-scale population screening and implementation in resource-constrained clinical settings.\nPoint 8\nFigure 2A appears to show less centralized plots across samples compared to DELFI. Even in the healthy group, the patterns seem less centralized compared to the DELFI healthy reference. It would be helpful for the authors to provide a rationale or explanation for this observation.\nResponse: We thank the reviewer for this thoughtful observation. To clarify this point, we performed an additional comparison and adjusted the visualization accordingly. The apparent “less centralized” pattern in Figure 2A arises primarily from differences in scaling and normalization between TuFEst and DELFI. In original Figure 2A, the visualization used z-score–standardized values, which expanded the dynamic range of signals across genomic bins and made sample-to-sample variation more visually pronounced—even among healthy controls. To address this, we re-generated the plots using raw coverage values instead of z-scores. As shown in therevised Figure 2A, this adjustment produces more centralized patterns without altering the underlying biological signal or affecting model performance.\nPont 9\nAt lines 187–189, the manuscript suddenly shifts to discussing transcription factors (TFs), which may be confusing readers with the rationale for incorporating TF features. Providing additional background or a brief explanation would improve readers understanding. Similarly, at line 398, it is hard to get reasons why RNA-seq was chosen for downstream analysis, especially given the availability of other molecular data types such as methylation, SNVs, and CNVs. Although the rationale becomes somewhat clearer later in the Discussion (starting at line 478), the authors should consider introducing the reasoning earlier in the manuscript. A brief explanation of why xCell and ESTIMATE were selected—and how they support the study objectives—would also be helpful.\nResponse: We thank the reviewer for highlighting the need for clearer conceptual transitions. We have now revised the manuscript to provide explicit background and rationale at each of the indicated points.\nFirst, we added a bridging explanation to clarify the rationale for incorporating transcription factor (TF)–associated features. Because cfDNA fragmentation patterns reflect nucleosome positioning and transcription factor occupancy, alterations in TF binding can leave characteristic footprints on plasma DNA. Prior fragmentomic studies have demonstrated that TF-site coverage is informative of regulatory activity and can enhance cancer detection performance (PMID: 31604930; 41043439). We now introduce this mechanistic rationale earlier in the manuscript to improve continuity when TF-related features are first presented (Lines 196–200 in the revised version).\nSecond, we clarified the reasoning behind selecting RNA-seq for downstream biological interpretation, despite the existence of alternative molecular data types such as methylation profiles, SNVs, and CNVs. Fragmentation patterns are closely linked to chromatin accessibility and transcriptional output; thus, RNA-seq provides a mechanistically grounded, orthogonal perspective for interpreting the biology underlying different TuFEst Cancer Scores. We have now introduced this rationale where RNA-seq analyses first appear, rather than reserving it for the Discussion (Lines 426–429 in the revised version).\nFinally, we added a brief explanation of why xCell and ESTIMATE were chosen for immune and stromal deconvolution (Lines 436–437 in the revised version). These tools are widely used and well-validated for inferring tumor microenvironment composition from bulk RNA-seq, and they directly support our objective of connecting cfDNA fragmentation signals with underlying immune and stromal states.\nPoint 10\nIn the manuscript showed the detection of lymph node metastasis. It would be helpful for the authors to clarify whether lymph node metastasis could be detected within a unified model along with other types of metastasis (Figure 4). If this is not feasible, a discussion explaining the limitations or reasons for separating these models would strengthen the manuscript.\nResponse: We sincerely thank the reviewer for this thoughtful comment and appreciate the opportunity to clarify this point. We apologize for any ambiguity in the original presentation that may have contributed to misunderstanding. Figure 4 does not present a metastasis-prediction model. Instead, it reports the performance of TuFEst in predicting the molecular subtypes of metastatic lesions as defined by ER, PR, and HER2 status. Consequently, Figure 4 is unrelated to lymph-node or distant metastasis.\nRegarding LN metastasis, our current analysis focuses specifically on identifying patients who are clinically node-negative (cN0) in early-stage breast cancer. The TuFEst-LN model demonstrates a high NPV, which is clinically meaningful because reliably excluding LN involvement may help reduce unnecessary invasive axillary procedures. The lower PPV, however, indicates that needle biopsy remains necessary when metastasis is suspected, consistent with standard practice.\nThe reviewer also raises an important point regarding whether LN metastasis and other metastatic sites could be predicted within a unified framework. At present, this is not feasible due to the limited number of patients with visceral or distant metastasis in our dataset, which constrains model development. Additionally, lymph-node metastasis and visceral/distant metastasis may generate distinct fragmentomic signatures due to differences in tumor burden and microenvironmental properties, further justifying separate modeling frameworks at the current stage. We agree that constructing an integrated metastasis-prediction model is an important and promising direction, and this will be explored in future studies as larger and more comprehensive datasets become available.\nPoint 11\nIn the Discussion section (line 436), the authors mention the limitations of bisulfite-based methods. However, more recent approaches such as enzymatic C-to-T conversion and Oxford Nanopore Technologies (ONT)-based methylation profiling offer viable alternatives that may mitigate such limitations. The authors should provide more precise, performance-based evidence to support their claims. Furthermore, performance comparisons with existing studies—ideally in tabular format—would enhance the manuscript. In particular, at line 464, the authors compare their model to studies in other cancer types. It would be more clear to include direct comparisons with recent liquid biopsy studies of breast cancer.\nResponse: We thank the reviewer for this valuable suggestion. In the revised Discussion (Lines 508–516), we now provide a more precise, performance-based summary of recent alternatives to bisulfite conversion, including EM-seq, TAPS, and ONT-based direct methylation profiling. In addition, we have created a new Supplementary Table 6 that systematically compares early breast-cancer detection performance and per-sample assay cost across major liquid biopsy platforms and imaging modalities. This table enables direct comparison between TuFEst and recently published breast-cancer liquid biopsy studies, addressing the reviewer’s request for clearer benchmarking. We have also updated the Discussion to contextualize TuFEst’s cost–performance characteristics relative to these existing methods, while taking care to emphasize that each modality has distinct advantages and trade-offs.\nPoint 12\nAt line 438, the authors refer to their approach as “AI-based.” However, given the methodological context, it may be more appropriate to describe it as “machine learning–based,” as the use of “AI” could be considered too advanced or potentially misleading in this case.\nResponse: We thank the reviewer for this important clarification. We agree that the term “AI-based” may overstate the methodological scope of our work and that “machine learning–based” more accurately reflects the approaches employed in this study. Accordingly, we have replaced “AI-based” with “machine learning–based” at the indicated location and have carefully reviewed the entire manuscript to ensure consistent and precise use of terminology throughout.\nPoint 13\nAt line 457, the authors mention using a tumor score, but what is the corresponding limit of detection in terms of variant allele frequency (VAF)? A brief discussion on this point would clarify the model’s sensitivity.\nResponse: We thank the reviewer for raising this important question regarding the model’s detection sensitivity. We agree that the limit of detection (LOD) in terms of variant allele frequency (VAF) is an essential parameter when interpreting the performance of sequencing-based assays.\nExisting literature provides useful benchmarks for VAF LOD in targeted sequencing panels. For example, Petrackova et al. (PMID: 31552176), using a binomial probability model, recommend a minimum depth of ≥1,650× and at least 30 mutant reads to reliably detect variants at ≥3% VAF. Jennings et al. (PMID: 28341590; 33486075) similarly suggest coverage levels of 250–500× for clinical-grade targeted NGS.\nHowever, these standards apply specifically to deep targeted sequencing. The low-pass whole-genome sequencing (LP-WGS) approach used in our study (average coverage ~2×) operates under a fundamentally different technical principle. At this coverage, the “tumor score” does not rely on detecting low-frequency single nucleotide variants (SNVs) or directly computing a VAF. Instead, it reflects tumor burden by capturing larger-scale genomic features—such as copy-number alterations (CNAs) and fragmentation-based chromatin signals—that can be robustly inferred even at ultra-low coverage.\nFor this reason, conventional VAF-based LOD benchmarks (e.g., 3% VAF thresholds) are not applicable to LP-WGS fragmentomic analyses, and we have added a brief clarification in the revised manuscript accordingly (Lines 518–519).\nReviewer #2 (Remarks to the Author):\nThe manuscript introduces a new multi-feature cfDNA-based Breast cancer detection method. The method is trained and tested on a cohort of 503 cases and 289 controls and shows good classification performance. Additionally, they train a ML model that can perform molecular subtyping using the same features. The results are pretty good and of interest to people in the field, but the methods section is lacking details about the features and how the method was trained.\nYou use three different kinds of fragmentomics features. How many features do you end up with of each kind?\nResponse: Thank you for your comment. We ended up with the following number of features for each fragmentomics feature type, including TF (148 features) (PMID: 36399356), DELFI (515 features) (PMID: 31142840), End-motif (256 features for 4BP End-motif and 6 features for 6BP End-motif) (PMID: 32111602) and Repeat elements (1,279 features) (PMID: 38478628).\nHave you tested which kind of fragmentomics features contribute the most to the performance of the method?\nResponse: We thank the reviewer for this insightful question. Clarifying the relative contributions of different fragmentomic feature classes is indeed critical for both biological interpretability and future assay optimization. In response, we conducted additional analyses to systematically evaluate the predictive performance of each modality independently. Specifically, we assessed models built using DELFI fragmentation profiles, end-motif signatures, and repeat-element features. As shown in the newly added Supplementary Figure 1, each modality exhibited detectable discriminatory power on its own.\nModel interpretability using SHapley Additive exPlanations (SHAP) further demonstrated that end-motif features contributed the largest mean absolute SHAP values, followed by repeat elements, TF features, and DELFI profiles (Rebuttal Figure 3). These results indicate that all feature classes provide complementary and nonredundant biological information, and their integration is responsible for the superior performance of the combined model.\nRebuttal Figure 3. SHAP-based interpretability analysis of different fragmentomic features.\nDescribing the motif analysis, you write that the 6-bp motifs are centered around the breakpoint. Is this also the case for the 4bp motifs, or are they only the end-motif?\nResponse: We thank the reviewer for this helpful clarification request. To clarify, the 4-bp motifs in our analysis are used exclusively as end-motifs (EDMs), representing the nucleotide sequences located precisely at the fragment termini. In contrast, the 6-bp motifs are defined as breakpoint-motifs (BPMs), which are centered around the predicted fragmentation breakpoints, as described in previous studies [PMID: 32111602; 35690859; 35780566]. Thus, the two motif classes reflect distinct biological signals: EDMs capture terminal sequence preferences at cleavage sites, whereas BPMs characterize the broader genomic context surrounding fragmentation breakpoints.\nIf you count motifs in both ends of the fragments, do you use the reverse complement of one of the ends so that the motifs are comparable?\nResponse: We thank the reviewer for raising this important methodological question. In standard end-motif analysis, reverse-complement transformation is often applied to one end to ensure direct comparability between the two fragment termini.\nIn our study, however, we intentionally analyzed the 5′ and 3′ end motifs separately, without applying reverse-complement merging. This design choice was motivated by the possibility that strand-specific cleavage preferences may exist, reflecting asymmetries in nuclease activity or DNA damage/repair processes. Treating the two ends independently preserves such potential biological signals, which would be lost if the ends were collapsed through reverse complementation.\nAccordingly, our 256 end-motif (EDM) features represent the full set of 4-bp motifs derived independently from both fragment ends, providing a more comprehensive view of possible asymmetric fragmentation patterns, as also suggested by recent fragmentomic studies (PMID: 32111602; 35690859; 35780566).\nThere are no details about the different ML methods that are trained provided in the manuscript. I looked at the linked GitHub repo, but as far as I can see, that repo only contains code for plotting. You should both add more details about the parameters and implementations of the methods to the manuscript and add the code to GitHub to make the details regarding the training and testing of the models available.\nResponse: We sincerely thank the reviewer for this valuable comment and fully agree that improving methodological transparency is essential. In the revised manuscript, we have expanded the Methods section (Lines 678–719 in the revised manuscript) to provide detailed descriptions of the machine-learning workflow. Specifically, we now clarify that: “A SEM integrated predictions from five base learners-SVM, GLM, GBM, XGBoost, and MLP-using a logistic regression meta-learner. The pipeline included feature standardization and was evaluated through 10-fold stratified cross-validation. Performance metrics were calculated for each fold, and probability thresholds were established at specificity levels of 0.9 and 0.95. The ensemble demonstrated strong discriminative ability (AUC = 0.943), with all predictions and sample information systematically recorded. The complete pipeline was saved to support potential clinical implementation.”\nTo further ensure full transparency and reproducibility, we have uploaded the complete machine-learning codebase—including training, evaluation, and prediction scripts—to a dedicated GitHub repository:\nhttps://github.com/yanwenfei0926/BC_CFDNA/tree/main/BC_Model/\nThis repository now contains all scripts necessary to reproduce the results presented in the manuscript, addressing the reviewer’s concerns.\nFor instance, is it a linear SVM you test, or do you use a kernel? Do you use cross-validation to fit any meta-parameters for any of the ML methods, or do you use default parameters? You end up using a Stacked Ensemble ( I would write out the abbreviation the first time you use it), but that is just a method for combining other ML models, and you do not say anything about which models you are combining (is it the 5 other mentioned models?).\nResponse: We sincerely thank the reviewer for these important clarifying questions. We have now revised the Methods section to explicitly provide the requested technical details.\nFirst, regarding the SVM implementation: we used a linear SVM, consistent with prior work demonstrating improved stability in high-dimensional settings and to reduce the risk of overfitting with limited samples. We now explicitly report this choice in the revised manuscript (Lines 685–692).\nSecond, for all machine-learning methods—including SVM, GLM, XGBoost, MLP, and GBM—we employed a nested cross-validation strategy to tune hyperparameters. Specifically, a 10-fold stratified outer loop was used for performance evaluation, while an inner loop performed grid-search–based hyperparameter optimization on each training fold. This procedure reduces optimistic bias and ensures that the reported metrics reflect generalizable performance. All hyperparameter grids and training scripts are available in the GitHub repository:\nhttps://github.com/yanwenfei0926/BC_CFDNA/tree/main/BC_Model/\nThird, we appreciate the reviewer’s suggestion to clarify the composition of the Stacked Ensemble Model (SEM) and to write out the abbreviation at first use. As now stated in the revised text, the SEM combines predictions from the five base learners evaluated in this study (GLM, linear SVM, XGBoost, MLP, and GBM). A logistic regression meta-learner was then trained on the out-of-fold predictions from these base models to generate the final ensemble prediction. We have updated the manuscript to clearly describe this structure (Lines 713–719).\nWe thank the reviewer again for these insightful questions, which have helped us substantially improve the clarity and transparency of our methodological description.\nThe method is called TuFEst (Tumor Fraction Estimator), but the manuscript only reports on classification performance; there are no plots showing an estimated tumor fraction compared to some independent estimate of the tumor fraction in the cfDNA samples.\nResponse: We thank the reviewer for this insightful observation. We fully acknowledge that the name “Tumor Fraction Estimator (TuFEst)” may suggest direct quantification of tumor fraction. In the current study, however, TuFEst is used primarily as a classification model, and its output—referred to as the cancer score (Figure 3H–J)—represents a probability value (cutoff = 0.5) used for malignancy prediction rather than an explicit estimate of the percentage of tumor-derived cfDNA.\nAt present, we do not have an independent ground-truth measure of tumor fraction for these clinical samples (e.g., through deep targeted sequencing or orthogonal cfDNA assays), and therefore cannot provide direct tumor fraction–to–tumor score comparisons. For this reason, we have focused on reporting classification performance. We retained the name “TuFEst” because the cancer score reflects a surrogate measure of tumor-derived DNA burden, derived from fragmentomic features known to correlate with tumor fraction in prior studies. This framework provides the conceptual basis for extending TuFEst into a true quantitative tumor fraction estimator in future work. We have revised the manuscript to clarify this distinction and avoid ambiguity (Lines 672–677in the revised version).\nReviewer #2 (Remarks on code availability):\nAs far as I can see, that repo only contains code for plotting. I could not find code for training the different models mentioned in the manuscript.\nResponse: Thank you for carefully examining our code repository and for raising this important point. We apologize for the earlier omission. We have now uploaded the complete code for training all models described in the manuscript to the following directory in the GitHub repository:\nhttps://github.com/yanwenfei0926/BC_CFDNA/tree/main/BC_Model/\nReviewer #4 (Remarks to the Author):\nThis is a well written manuscript demonstrating the power of cfDNA fragments to detect breast cancer as well as determine subtype and lymph node status.\nIt was a well conducted study using both benign and malignant breast tissue.\nI think that this study will be of interest to researchers and clinicians involved in early breast cancer detections and treatment.\nI have very minor concerns:\n1. I could not locate Table 1 cited on line 209.\nResponse: First, we sincerely thank the reviewer for the positive and encouraging comments. We truly appreciate the reviewer’s recognition of our study design. In terms of this issue, we apologize for the oversight regarding Table 1, which has now been moved to the correct location (Line 1093). We have carefully checked all table and figure placements to ensure consistency between the text and the manuscript layout.\n2. there is no discussion of tumors which are triple positive, whether they were included or excluded in the analyses.\nResponse: We thank the reviewer for this helpful comment. In our study, molecular subtyping was based on ER, PR, and HER2 expression status, following standard therapeutic categories (HR⁺/HER2⁻, HER2⁺, and TNBC). Under these commonly used clinical definitions, triple-positive breast cancers (TPBC; ER⁺/PR⁺/HER2⁺) are not treated as a separate subtype but are included within the HER2⁺ group, as HER2-targeted therapy primarily determines treatment stratification.\nIn our cohort, the number of TPBC cases was relatively small, which further limited the statistical power to train and evaluate a dedicated TPBC class. For these reasons, TPBC tumors were grouped within the HER2⁺ subtype during model development. We have now explicitly clarified this in the revised manuscript.\nWe also acknowledge this as a limitation (Lines 576–578in the revised version) and note that future studies with larger cohorts will allow us to assess whether cfDNA fragmentomic signatures can further distinguish TPBC as a biologically or clinically distinct subgroup. We appreciate the reviewer’s suggestion, which has helped improve the clarity of our study design and highlighted an important direction for future research.\nRESPONSE TO REVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nThe inclusion of details about the model training as well as the feature importance analysis is a clear improvement.\nBut there are two minor points I still feel need to be addressed:\n1) The explanation of the 4 bp end motifs and 6bp breakpoint motifs in the manuscript remains cryptic. The clarification you wrote in the response to my comment is not reflected in the manuscript. In the graphical abstract, the 4bp and 6bp motifs are drawn the same way, and in the text, you do not distinguish between end motifs and breakpoint motifs.\nResponse: Thank you for the reviewer’s suggestion. We have carefully revised the graphical abstract to clarify the difference between the 4-bp end motifs and 6-bp breakpoint motifs, and have also updated the corresponding description in the main text. We appreciate the valuable feedback.\nIn line 252 you refer to them the same way: “encompassing end motifs (4/6bp breakpoint signatures)”.\nResponse: Thank you for pointing this out. We have revised the phrasing to \"encompassing end motifs (4 bp) and breakpoint motifs (6 bp)\" to clearly distinguish between the two. We appreciate your careful reading.\nAnd line 285: “particular 4 bp (e.g.,AAGC, AGGC) and 6 bp (e.g., CTTGTT, TGGGAC) end motifs”.\nResponse: Thank you for pointing this out. We have revised the text to \"particular 4 bp (e.g., AAGC, AGGC) end motifs and 6 bp (e.g., CTTGTT, TGGGAC) breakpoint motifs\" to ensure clarity and distinction. We appreciate your careful reading.\nAnd line 348: “patterns in breakpoint end motifs”.\nResponse: Thank you for your careful review. We have revised the text to \"patterns in end motifs, breakpoint motifs, \" to accurately and clearly distinguish between the two concepts.\nAnd line 641: “4 bp and 6 bp end motif frequencies”.\nResponse: Thank you for pointing this out. We have revised the phrasing to \"4 bp end motif and 6 bp breakpoint motif frequencies\" to maintain clear and consistent distinction between the two terms throughout the manuscript.\n2) I am surprised to see that you include so many Repetitive Element features, not just the main classes, but many subclasses. I can understand why Alu elements could be a relevant feature if they are actively transposing in the tumor tissue, but it is only a subset of the Alu subfamilies that can still be active. Furthermore, the results in Supp. Figure 2D shows that the main difference between benign and malignant samples is that there are more Alu reads in a subset of the benign samples. What could explain this? The figure text conflicts with the legend. Line 1065: “breast malignancy (blue) and benign controls (pink)” but the legends specify the colors the other way around. Which is correct?\nResponse: Thank you for your insightful comment. Our ARTEMIS pipeline analyzes a comprehensive set of repetitive elements, including not only Alu subfamilies but also other classes such as MERs, satellite repeats, LTR elements, and Tigger subfamilies, allowing for a broader feature space.\nRegarding the observation that a subset of benign samples shows a higher Alu read count, we agree that this is an interesting finding. One plausible explanation is the natural inter-individual variation in repetitive element content within the general population, where some individuals (and thus some benign samples) may inherently have higher baseline levels of certain repeats. The specific biological mechanisms driving this observation would indeed require further, more targeted investigation. We appreciate your valuable perspective on this point.\nWe have corrected the description in line 1065 to \"cfDNA fragmentomic features between breast malignancy (pink) and benign controls (blue)\" to ensure consistency with the figure colors.\nReviewer #2 (Remarks on code availability):\nIt is not reproducible.\nThere is no README and no description of the input files (not included). There are some comments in the code but they are all in chinese.\nResponse: Thank you for your feedback. We sincerely apologize for the oversight. We have updated the code repository to include a comprehensive README file in English. We believe that researchers with basic coding skills will be able to set up and run the code accordingly.\nRegarding the full input data, the proprietary, individual patient-level training data cannot be made publicly available due to commercial confidentiality and pending patents. This is clarified in the Data Availability statement. The provided code, example data, and documentation ensure the reproducibility of the analytical methodology."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Dear Editor,\nI have received and carefully reviewed the authors’ responses to reviewers comments by thoroughly addressing the points with appropriate justification and supporting results. Based on these revisions, the manuscript is appropriate for publication in Nature Communications and will provide valuable insights to researchers in the relevant field.", "reviewer_2": "The inclusion of details about the model training as well as the feature importance analysis is a clear improvement.\nBut there are two minor points I still feel need to be addressed:\n1) The explanation of the 4 bp end motifs and 6bp breakpoint motifs in the manuscript remains cryptic. The clarification you wrote in the response to my comment is not reflected in the manuscript. In the graphical abstract, the 4bp and 6bp motifs are drawn the same way, and in the text, you do not distinguish between end motifs and breakpoint motifs.\nIn line 252 you refer to them the same way: “encompassing end motifs (4/6 bp breakpoint signatures)”.\nAnd line 285: “particular 4 bp (e.g.,AAGC, AGGC) and 6 bp (e.g., CTTGTT, TGGGAC) end motifs”.\nAnd line 348: “patterns in breakpoint end motifs” .\nAnd line 641: “4 bp and 6 bp end motif frequencies”.\n2) I am surprised to see that you include so many Repetitive Element features, not just the main classes, but many subclasses. I can understand why Alu elements could be a relevant feature if they are actively transposing in the tumor tissue, but it is only a subset of the Alu subfamilies that can still be active. Furthermore, the results in Supp. Figure 2D shows that the main difference between benign and malignant samples is that there are more Alu reads in a subset of the benign samples. What could explain this? The figure text conflicts with the legend. Line 1065: “breast malignancy (blue) and benign controls (pink)” but the legends specify the colors the other way around. Which is correct?", "reviewer_4": "the authors have addressed all of my comments in this revised version"}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nThe inclusion of details about the model training as well as the feature importance analysis is a clear improvement.\nBut there are two minor points I still feel need to be addressed:\n1) The explanation of the 4 bp end motifs and 6bp breakpoint motifs in the manuscript remains cryptic. The clarification you wrote in the response to my comment is not reflected in the manuscript. In the graphical abstract, the 4bp and 6bp motifs are drawn the same way, and in the text, you do not distinguish between end motifs and breakpoint motifs.\nResponse: Thank you for the reviewer’s suggestion. We have carefully revised the graphical abstract to clarify the difference between the 4-bp end motifs and 6-bp breakpoint motifs, and have also updated the corresponding description in the main text. We appreciate the valuable feedback.\nIn line 252 you refer to them the same way: “encompassing end motifs (4/6bp breakpoint signatures)”.\nResponse: Thank you for pointing this out. We have revised the phrasing to \"encompassing end motifs (4 bp) and breakpoint motifs (6 bp)\" to clearly distinguish between the two. We appreciate your careful reading.\nAnd line 285: “particular 4 bp (e.g.,AAGC, AGGC) and 6 bp (e.g., CTTGTT, TGGGAC) end motifs”.\nResponse: Thank you for pointing this out. We have revised the text to \"particular 4 bp (e.g., AAGC, AGGC) end motifs and 6 bp (e.g., CTTGTT, TGGGAC) breakpoint motifs\" to ensure clarity and distinction. We appreciate your careful reading.\nAnd line 348: “patterns in breakpoint end motifs”.\nResponse: Thank you for your careful review. We have revised the text to \"patterns in end motifs, breakpoint motifs, \" to accurately and clearly distinguish between the two concepts.\nAnd line 641: “4 bp and 6 bp end motif frequencies”.\nResponse: Thank you for pointing this out. We have revised the phrasing to \"4 bp end motif and 6 bp breakpoint motif frequencies\" to maintain clear and consistent distinction between the two terms throughout the manuscript.\n2) I am surprised to see that you include so many Repetitive Element features, not just the main classes, but many subclasses. I can understand why Alu elements could be a relevant feature if they are actively transposing in the tumor tissue, but it is only a subset of the Alu subfamilies that can still be active. Furthermore, the results in Supp. Figure 2D shows that the main difference between benign and malignant samples is that there are more Alu reads in a subset of the benign samples. What could explain this? The figure text conflicts with the legend. Line 1065: “breast malignancy (blue) and benign controls (pink)” but the legends specify the colors the other way around. Which is correct?\nResponse: Thank you for your insightful comment. Our ARTEMIS pipeline analyzes a comprehensive set of repetitive elements, including not only Alu subfamilies but also other classes such as MERs, satellite repeats, LTR elements, and Tigger subfamilies, allowing for a broader feature space.\nRegarding the observation that a subset of benign samples shows a higher Alu read count, we agree that this is an interesting finding. One plausible explanation is the natural inter-individual variation in repetitive element content within the general population, where some individuals (and thus some benign samples) may inherently have higher baseline levels of certain repeats. The specific biological mechanisms driving this observation would indeed require further, more targeted investigation. We appreciate your valuable perspective on this point.\nWe have corrected the description in line 1065 to \"cfDNA fragmentomic features between breast malignancy (pink) and benign controls (blue)\" to ensure consistency with the figure colors.\nReviewer #2 (Remarks on code availability):\nIt is not reproducible.\nThere is no README and no description of the input files (not included). There are some comments in the code but they are all in chinese.\nResponse: Thank you for your feedback. We sincerely apologize for the oversight. We have updated the code repository to include a comprehensive README file in English. We believe that researchers with basic coding skills will be able to set up and run the code accordingly.\nRegarding the full input data, the proprietary, individual patient-level training data cannot be made publicly available due to commercial confidentiality and pending patents. This is clarified in the Data Availability statement. The provided code, example data, and documentation ensure the reproducibility of the analytical methodology."}]

Decision letter