A Computational Framework for Analysis of cfDNA Fragmentation Profiles

Circulating cell-free DNA (cfDNA) has emerged as a promising non -invasive medium for studying tumor molecular profiles. Non -random fragmentation patterns in plasma cfDNA, particularly around nucleosome -depleted regions (NDRs) near transcription start sites (TSS), have been shown to reflect epigenetic regulation and gene expression. In this study, coverage profiles of the NDR were utilized to derive an NDR score, which was subsequently used as a proxy for inferring gene expression. To reduce transcript-to-transcript variability and enhance the clarity of these expression - associated signals, we implement a method for GC-bias correction of cfDNA samples. A computational framework (NDRDiff) was then developed to enable comparative analyses of NDR score profiles across different sample groups. The GC-bias correction preserved the overall trend of the NDR signal while improving the separation of gene expression levels, as demonstrated by comparisons of healthy donor cfDNA samples with matched blood RNA -seq data. Validation on a simulated dataset showed that NDRDiff achieved an area under the precision –recall curve (AUPRC) of 0.916, outperforming a standard t-test (AUPRC of 0.777). When applied to a comparison of healthy donor cfDNA and metastatic colorectal cancer (mCRC) cfDNA, NDRDiff identified 531 differential NDR score (DNS) genes that facilitated clear separation between the two groups. These DNS genes were found to correlate with tumor fraction estimates (down -regulated DNS genes: Pearson R = 0.89, p < 0.05; up -regulated DNS genes: Pearson R = –0.88, p < 0.05) and included CLDN4, BIN2, and IRAG2, which exhibit strong associations with colorectal cancer or blood cell expression signatures. Gene set enrichment analysis further revealed enrichment of colon and other gastrointestinal tissue signatures. Collectively, these findings underscore the potential of NDR -based cfDNA analysis as a minimally invasive tool for monitoring tumor-related molecular features in cancer.

cfDNA, anti-EGFR, Colorectal Cancer, Resistance .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 3 Contents Abstract................................ ................................ ................................ ................................ . 2

................................ ................................ ................................ ........................... 4

................................ ................................ ................................ ................................ .. 6 Inference of gene expression with NDR Score ................................ ................................ .. 6 GC-bias correction improves NDR signals ................................ ................................ ......... 8 Identification of differential NDR score (DNS) ................................ ................................ .... 9 Comparison of Healthy vs CRC samples with NDR score ................................ ............... 11

................................ ................................ ................................ .......................... 15 Methodology ................................ ................................ ................................ ....................... 17 Sample Collection ................................ ................................ ................................ ........... 17 Library Preparation ................................ ................................ ................................ .......... 17 Nucleosome-Depleted Region (NDR) Score Estimation ................................ .................. 18 Differential Comparison of NDR Scores ................................ ................................ ........... 19 Tissue-Specific Gene Set Selection ................................ ................................ ................. 20 Gene Set Enrichment Analysis ................................ ................................ ........................ 21 Statistics ................................ ................................ ................................ .......................... 21 Supplementary Materials ................................ ................................ ................................ .... 22 Supplementary Figure 1 ................................ ................................ ................................ .. 22 Supplementary Figure 2 ................................ ................................ ................................ .. 23 Supplementary Figure 3 ................................ ................................ ................................ .. 24 Acknowledgment ................................ ................................ ................................ ................. 25

................................ ................................ ................................ ........................... 26 .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 4

Colorectal cancer (CRC) is one of the most prevalent cancers globally, accounting for about 10% of all diagnosed cancers and is the second leading cause of cancer-related deaths1. Approximately 22% of newly diagnosed CRC patients present with distant metastases2, and up to 70% of CRC patients develops metastatic disease3. Although advances in treatment have improved outcomes for selected patients, survival rates for metastatic CRC (mCRC) remains limited, with a 5-year relative overall survival (OS) of approximately 15%4,5. This underscores the importance of comprehensive management strategies, including tailored treatments based on molecular characteristics, to improve clinical outcomes for patients. Understanding the evolution of the tumor molecular profile throughout the treatment period is critical for understanding the mechanisms of treatment resistance, enabling the early detection of resistance, and identifying opportunities for alternative therapies. However, t his necessitates longitudinal sampling at multiple time points during treatment. Traditional methods for detecting these alterations often require invasive tissue biopsies, which are not only poses significant risk for patients but also may not capture the full heterogeneity of the tumor 6. A liquid biopsy approach, such as the study of plasma cell-free DNA (cfDNA), presents an alternative to tissue biopsy due to the minimally invasive nature, allowing longitudinal sampling to track tumor evolution. Plasma cfDNAs are DNA released from cells into the blood and are released through diverse cell -death modalities, such as apoptosis and necrosis and active excretion mechanisms7. While predominantly originating from hematopoietic cells 8, elevated plasma cfDNA is detected in cancer patients in 19779 and since then, numerous studies have confirmed the contribution of tumor-derived DNA to circulating cfDNA10– .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 5 12, referred to as circulating tumor DNA (ctDNA). Plasma ctDNA can retain genetic and epigenetic profile of tumor and has been used to detect genetic mutations9,13,14 and epigenetic alterations, such as DNA methylation11,12,15 and histone modification16. Plasma cfDNA was found to show non -random fragmentation pattern that reflects epigenetic regulation17. Subsequently, differences in nucleosome protection patterns relative to the transcription start site (TSS) are observed between highly and lowly expressed genes. This pattern, or nucleosome footprint , around the TSS reflects chromatin accessibility possibly arising from the presence of the transcription preinitiation complex in transcriptionally active genes 18. Later studies confirmed the finding and showed that the relative coverage of the nucleosome -depleted region (NDR) positioned from -150 to +50 bp relative to the TSS (henceforth referred to as NDR score) can be employed for the classification of expressed and silent genes 19, and for the estimation of plasma cfDNA tumor fraction20. Considering the NDR score's association with gene expression, we investigated the potential of using NDR score as a metric to infer gene expression profiles of the cellular populations contributing to plasma cfDNA. In this study, we adapted existing GC-bias estimation approach to correct GC-bias in cfDNA samples . Building on these enhanced NDR estimates, we develop a computational framework for comparative analysis of NDR score profiles across different sample groups. We evaluate the performance of this framework against standard statistical test by validating it on a simulated dataset, and further demonstrate its utility by comparing healthy donor cfDNA samples with those from mCRC patients. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 6

Deep whole-genome sequencing (dWGS) of plasma cfDNA was performed, achieving a sequencing coverage of 60x for healthy donors and 120x for mCRC samples (Methods). In total, t he study included 10 mCRC patients, with 31 mCRC d WGS samples, and 17 dWGS samples from healthy donors (Table 1). Table 1 Clinical characteristics of the dWGS cohort, including healthy and mCRC samples. Healthy (n=17) mCRC (n=31) Age Median (Min - Max) 39.5 (23 - 52) 57.2 (40.2 - 68) Gender Male 2 7 Female 10 3 Microsatellite Stability Clinical MSS NA 7 Not Available NA 3 Stage III NA 1 IV NA 9 Tumor Location Left Colon NA 4 Right Colon NA 2 Rectum NA 4 Tumor Fraction Median (Min - Max) NA 0.225 (0 - 0.506) Inference of gene expression with NDR Score As hematopoietic cells are the main contributors of plasma cfDNA in healthy individuals21, we investigated the association between NDR score and gene expression by comparing the relative coverage of TSS in our healthy samples (n=17) with corresponding gene expression from RNA-Seq data of blood samples. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 7 Figure 1 Association of NDR score with gene expression and improvement of NDR signals with GC -bias correction. (a) Schematic of GC-bias correction procedure for NDR score estimation. (b, e) Mean relative coverage around the TSS grouped for (b) original and (e) GC-bias corrected sequencing coverage according to gene expression (TPM) groups based on GTEx Whole Blood RNA -Seq for healthy cfDNA samples (n=17). The ribbon indicates the mean relative coverage ± 2 standard deviation (SD). The NDR site region is indicated b y the dotted lines (-150 to +50 bp). (c, f) Comparison of mean NDR score across blood gene expression groups for (c) original and (f) GC-bias corrected NDR scores. (d) GC-bias corrected (y-axis) compared to original mean NDR scores (x- axis). (g) Comparison of mean absolute deviation (MAD) of original vs GC -bias corrected NDR score for silenced (TPM ≤ 0.1) and expressed genes (TPM > 10). .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 8 Transcripts were categorized based on transcripts per million (TPM) values extracted from GTEx Whole Blood RNA -Seq. Comparing the mean relative coverage profiles across each TPM group , we demonstrated a reduction in relative coverage surrounding the TSS region associated with increased gene expression levels (Figure 1b), consistent with previous studies19. Specifically, we observed a decrease in mean relative coverage within the central region (-150 to +50 bp with respect to TSS) in actively transcribed genes compared to silenced genes. Subsequently, we computed the NDR Score (Methods) and explored its association with TPM levels using data from the GTEX Whole Blood RNA -Seq dataset (Supplementary Figure 1a). We identified a moderate inverse correlation between the NDR Score and gene expression levels (Spearman’s ρ = -0.614, p < 0.05). Notably, we observed distinct clusters of points representing silenced and expressed genes, particularly around NDR Scores of 1.0 and 0.5. A clear separation between expressed and silent genes was evident. However, we observed considerable overlap in mean NDR score across TPM groups of varying expression levels, suggesting limitations in the NDR score's ability to finely discriminate gene expression levels (Figure 1c). GC-bias correction improves NDR signals GC-bias refers to the correlation between read coverage and GC content of the sequenced region. This bias arises during library preparation, notably in PCR amplification and binding-based purification steps, leading to an underrepresentation of regions with extreme GC content22. GC-bias can obscure the desired signals in the analysis of sequencing data related to quantifying of coverage within a genome , and has been applied in cfDNA analysis to improve signals for detection of copy number .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 9 alterations23 and assessment of chromatin accessibility 24. To account for GC bias in the analysis of NDR scores, we adapted an existing GC-bias estimation method24 for the correction of sequencing coverage at the base -pair level. This approach enables the calculation of NDR scores using corrected sequencing coverage (Methods, Figure 1a). Following GC -bias correction, a robust correlation between original and GC -bias corrected NDR scores ( Spearman’s correlation, R = 0.91; Figure 1d), indicating that the correction process preserved the underlying biological signal. In addition, consistent trends in relative coverage at NDR sites across varying gene expression levels were maintained , accompanied by a significant reduction in the coverage variability across gene expression groups ( Figure 1e, f; Supplementary Figure 1b). Comparing mean NDR scores before and after GC -bias correction, an improved separation of different gene expression groups was observed, attributing to the reduction in variation of NDR scores within each group (Figure 1f). This is reflected by a significant decrease in mean absolute deviation (MAD) of the NDR scores observed post-correction (Figure 1g, p < 0.01). Comparison of NDR scores between pre-defined blood- and CRC-specific gene sets (Methods) in CRC and healthy samples shows a slight improvement in the separation of mean NDR scores between the two gene sets in healthy cfDNA samples (Supplementary Figure 1c, d). Identification of differential NDR score (DNS) To detect tissue -specific and treatment -associated changes in NDR scores across samples, we evaluated differences in NDR scores between groups based on sample type and treatment status. Given that NDR score differences follow a normal .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 10 distribution (Supplementary Figure 2b), a Z -test was implemented for comparative analysis. To improve the accuracy of differential NDR score (DNS) analysis, standard deviation shrinkage was applied to address heteroskedasticity before performing the Z-test (Methods, Figure 2a), a method hereafter referred to as NDRDiff. Table 2 Performance metrics of NDRDiff and t-test for differential NDR score (DNS) analysis To evaluate the performance of NDRDiff compared to standard t-test, a simulated dataset was generated for validation . Based on the observation that the overall distribution of NDR scores appear to follow a mixed-effect normal distribution (Figure 2b), we generated 10 simulated datasets using bimodal parameters (Model 1: μ=0.459, σ=0.098, λ=0.354 ; Model 2 : μ=0.886, σ=0.188, λ=0.646 ; Methods), producing a distribution in the mixed-effect model that mirrors expected profiles (Figure 2c). Each dataset included 500 simulated DNS genes. Comparison of the paired groups in the simulated datasets with NDRDiff identified 532 (range 520-542) DNS genes (FDR < 0.05, |ΔNDR| ≥ 0.2) and demonstrated strong predictive accuracy, with a Area Under Precision-Recall Curve (AUPRC) of 0.916 (range 0.896 – 0.925) (Figure 2e, Table 2). In contrast, the t-test identified 506 (range 486-526) DNS genes (p < 0.01, |ΔNDR| ≥ 0.2), with a lower AUPRC of 0.777 (0.721 – 0.800) (Figure 2d and Table 2). Overall, standard deviation-shrinkage prior to Z-test enhanced the sensitivity for detecting DNS genes relative to t-test while maintaining high specificity. NDRDiff FDR < 0.05 FDR < 0.05 p < 0.01 p < 0.05 Positive Predictive Value (PPV) 0.805 (0.793-0.816) 0.973 (0.952-0.993) 0.7 (0.68-0.719) 0.348 (0.343-0.354) Negative Predictive Value (NPV) 0.996 (0.996-0.997) 0.977 (0.977-0.978) 0.993 (0.992-0.993) 0.999 (0.998-0.999) Sensitivity 0.856 (0.849-0.863) 0.0776 (0.0573-0.0979) 0.708 (0.692-0.724) 0.945 (0.94-0.951) Specificity 0.995 (0.994-0.995) 1 (1-1) 0.992 (0.992-0.993) 0.956 (0.955-0.957) Significant Genes Detected 532 (520-542) 40 (29-51) 506 (486-526) 1357 (1334-1380) Area Under Precision-Recall Curve (AUPRC) 0.916 (0.896 - 0.925) Performance Metrics t-test 0.777 (0.721 - 0.8) .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 11 Figure 2 Evaluation of statistical tests for detection of differential NDR score (DNS) genes. (a) Schematic of procedure for NDRDiff (Z-test with shrinkage of standard deviation). Created in BioRender. Skanderup, A. (2025) https://BioRender.com/7zlj88e. (b, c) Distribution of NDR scores in the (b) healthy cfDNA samples (n=17) and (c) simulated dataset (n=10). (d, e) Precision-Recall curve for the detection of significant DNS genes in (d) NDRDiff and (e) t-test. Comparison of Healthy vs CRC samples with NDR score To assess the feasibility of using NDR scores to identify differential molecular profiles between sample groups, our approach was initially applied to compare healthy (n=17) and mCRC (n=31) samples. A total of 53 1 DNS genes were identified (FDR< 0.05, |ΔNDR| ≥ 0.2), with 419 genes showing a decrease in NDR scores in mCRC samples .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 12 (up-regulated in mCRC) and 112 genes showing an increase in NDR scores in mCRC samples (down-regulated in mCRC) (Figure 3a, Supplementary Table 1). The identified DNS genes effectively distinguished healthy donors from mCRC samples, and the degree of separation appeared to be associated with the tumor fraction of the cancer samples (Figure 3b, c). Furthermore, a comparison of the overlapping transcripts from the identified DNS genes with Blood and CRC transcriptomic datasets confirm tissue-specific gene expression of the identified DNS genes (Supplementary Figure 3c). Among the DNS genes identified are CRC-related genes, such as CLDN4, and blood- related genes, such as BIN2 and IRAG2. CLDN4 is known to be overexpressed in various cancers, including CRC, with a study suggesting that targeting CLDN4 can enhance the anti-tumoral effects of chemotherapeutic agents 25. Consistent with this, CLDN4 was found to be activated in mCRC samples ( ΔNDR=-0.33, Supplementary Figure 3b) compared to blood samples. Conversely, an increase in NDR scores was observed in blood-related genes, such as BIN2 and IRAG2 (Supplementary Figure 3b), in mCRC samples (ΔNDR=0.382, 0.301) . BIN2 and IRAG2 are both expressed in platelet cells and involved in platelet activation and aggregation26,27. Gene set enrichment analysis (GSEA) on MSigDB hallmark gene sets28 demonstrated suppression of interferon alpha and gamma responses, inflammatory responses, and JAK/STAT signaling pathways in mCRC samples, which are critical for immune activity and regulation in the blood (Supplementary Figure 3d). Conversely, we detected activation of E2F targets, which play a key role in cell cycle progression and are .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 13 frequently associated with cancer proliferation. GSEA of Human Protein Atlas (HPA) Tissue Gene Expression Profiles gene sets29 revealed significant activation of Figure 3 Validation of NDRDiff with comparison of healthy (n=17) and mCRC (n=31) samples . (a) Volcano plot of -log10 statistical significance (y-axis) against difference in NDR scores (x -axis) for NDRDiff comparison of healthy and mCRC samples. Significant DNS genes are indicated in red. (b) Heatmap of the NDR scores of the identified DNS genes across healthy and mCRC samples. The top annotation indicates the sample type and tumor fraction of the samples. Tumor fraction is not applicable to healthy samples and is greyed out correspondingly. (c) PCA plot of healthy (red) and mCRC (green) samples across identified DNS genes. (d) GSEA plot depicting the normalized enrichment scores (NES) of top activated and suppressed gene sets identified for HPA tissue mRNA gene sets. Activated gene sets show decrease in NDR scores and suppressed gene sets show increase in NDR scores. (e) Scatter plots of mean NDR scores of identified DNS genes against tumor fraction in healthy (red) and mCRC (green) samples. Correlation between mean NDR scores and tumor fraction was examined with Pearson .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 14 correlation. Top panel: DNS genes down -regulated in mCRC samples, bottom panel: DNS genes up -regulated in mCRC samples. gastrointestinal tissue gene sets in mCRC samples, including those from the colon, small intestine, duodenum, and rectum (Figure 3e, Supplementary Figure 3a). In summary, the identification of DNS genes and enriched gene sets related to CRC and blood cells validates that the comparison of NDR scores can distinguish between cancer and normal samples and retrieve tissue-specific expression profiles. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 15

By analyzing cfDNA NDR coverage profiles, gene expression patterns were inferred. To overcome the limitations of existing methods, GC -bias correction and enhanced statistical approaches were implemented to enrich relevant NDR score signals. The approach was validated using simulated datasets and comparisons between healthy and mCRC samples. This work represents a novel application of cfDNA analysis for the comparative analysis of NDR profiles across sample groups. The study of cfDNA fragmentation profiles and their association with gene expression has been previously investigated18,19. However, earlier analyses were largely confined to inference of gene expression states (expressed or silenced) and estimation of cell type contribution to cfDNA. In this study, we demonstrate that incorporating GC -bias correction and statistical methods accounting for the heteroskedasticity of NDR scores enhances the detection of differential NDR score genes. The genes identified using this approach are shown to be relevant across different contexts, including simulated datasets (Figure 2) and comparison between healthy and cancer samples (Figure 3). Prior studies on cfDNA analysis have employed gene signatures derived from orthogonal approaches, such as DNA methylation sequencing or RNA sequencing, to evaluate tissue-specific origins of cfDNA 21,30 and to estimate tumor burden 20. These studies have demonstrated the capability to specifically identify tissue-specific and tumor-specific cfDNA fragments, marking a notable advancement in the clinical utility of cfDNA. Consistent with these findings, our analysis comparing healthy and cancer samples identified DNS genes that distinguish between the two groups and shows a correlation with tumor fraction in cancer samples (Figure 3). Our approach enables the .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 16 identification of gene signatures specific to cfDNA NDR profiles, which may improve tracking of tumor contributions across different cell types and cancer types. The use of cfDNA for inferring gene expression remains in its early stages, and several considerations can be addressed to enhance the accuracy of NDR signals. First, NDR scores were estimated based on relative coverage within a fixed region spanning -50 bp to +150 bp relative to the TSS. However, the actual region of interest, characterized by a dip in the cfDNA profile, may not consistently align with this fixed region. This potential misalignment highlights the need for methods to identify and evaluate optimal NDR sites for each transcript based on their correlation with gene expression. Second, the predominance of cfDNA from blood-derived sources poses a significant challenge, as cancer-specific NDR signals can be obscured by contributions from blood signals. Developing robust approaches to quantify and acc ount for the influence of blood - derived cfDNA will be essential for improving the detection of cancer-specific signals. Lastly, cfDNA analysis primarily requires sequencing of a ±2 kb region around the TSS for each gene, making it feasible to integrate this approach into existing targeted sequencing panels or whole -exome sequencing platforms. Such integration could enable cost-effective evaluation of inferred gene expression and facilitate the broader adoption of cfDNA-based methods in clinical applications. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 17 Methodology Sample Collection This study was approved by the institutional review boards of SingHealth (2018/2795, 2019/2401), and all procedures were conducted in accordance with ethical guidelines. Written informed consent was obtained from all participants prior to sample collection. Venous blood samples were collected in EDTA tubes to prevent coagulation. Samples were either stored at -80°C for later analysis or processed within two hours of collection to maintain integrity. Whole blood was centrifuged at 300 g x 10 min and 9730 g x 10 min to isolate the plasma layer, which was subsequently stored at -80°C until further analysis. Library Preparation Plasma DNA was extracted using the QIAamp Circulating Nucleic Acid Kit ( QIAGEN, 55114) according to the manufacturer’s protocol . Library preparation was conducted with KAPA HyperPrep Kit (Roche, KK8504), using up to 100 ng of cfDNA each sample. Following end repair and A -tailing, the cfDNA was ligated with custom adapters containing a random 8-mer adjacent to the library index site, synthesized by Integrated DNA Technologies (IDT). Post-ligation cleanup was performed using 0.8X Agencourt AMPure XP Beads, and the adapter -ligated cfDNA was eluted in 20 µL of nuclease - free water. Library amplification was carried out using KAPA HiFi HotStart Ready Mix (Roche, KK2602) and 1X Agencourt AMPure XP Beads (Beckman Coulter, A63882), followed by elution of the amplified library in 20 µL of Elution Buffer (10 mM Tris-HCl, pH 8.0). Quantification of the plasma DNA libraries was performed using the KAPA Universal Library Quantification Kit (Roche, KK4824), and the library quality was assessed with Agilent High Sensitivity DNA Kit (Thermo Fisher Scientific, NC1738319). .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 18 Sequencing was carried out on an Illumina NovaSeq 6000 System, generating 2 × 150 bp paired-end reads. dWGS cfDNA samples from healthy donors were sequenced to ~60X coverage, whereas those from mCRC patients were sequenced to ~120X coverage. Nucleosome-Depleted Region (NDR) Score Estimation 𝑁𝐷𝑅 𝑆𝑐𝑜𝑟𝑒 = 𝐶𝑜𝑣̅̅̅̅̅𝑐𝑒𝑛𝑡𝑟𝑎𝑙 𝐶𝑜𝑣̅̅̅̅̅𝑓𝑙𝑎𝑛𝑘𝑖𝑛𝑔 Equation 1 where 𝐶𝑜𝑣̅̅̅̅̅𝑐𝑒𝑛𝑡𝑟𝑎𝑙 denotes the mean sequencing coverage in the central region surrounding the TSS ( -50 to +150 bp relative to the TSS), and 𝐶𝑜𝑣̅̅̅̅̅𝑓𝑙𝑎𝑛𝑘𝑖𝑛𝑔 represents the mean coverage in the upstream and downstream flanking regions (-2000 to -1000 bp and +1000 to +2000 bp relative to the TSS). The NDR score was estimated as described in Equation 1. The mean sequencing coverage of the central region surrounding the transcriptional start site (TSS , -50 to +150 bp relative to the TSS) was normalized by the mean coverage of the flanking regions (-1000 to -2000 bp upstream and +1000 to +2000 bp downstream). Only reads with a mapping quality score of ≥20 was included for the coverage calculation. NDR scores were estimated for selected transcripts in the Matched Annotation from NCBI and EMBL -EBI (MANE) Select set, where transcript at each genomic locus was selected to reflect the biology of the locus31. GC-bias correction was performed prior to NDR score estimation to account for biases arising from variations in GC content across different transcripts. Fragment -length- specific GC -bias estimates were generated using Griffin, with the provided Snakemake workflow file using default parameters (Doebley et al., 2022). Each sequencing read was assigned a GC-bias estimate based on its fragment length and GC content. Bias-corrected coverage was then calculated by applying the reciprocal .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 19 of the GC -bias factor (1/GC -bias) to ensure normalization of coverage across transcripts. Additionally, NDR sites were excluded if more than half of the samples in either the healthy or colorectal cancer cohorts exhibited an average flank coverage of zero or extreme NDR scores ( 2). Differential Comparison of NDR Scores To compare nucleosome -depleted region (NDR) scores across groups, such as healthy versus colorectal cancer (CRC) samples, a robust statistical approach was developed to account for observed heteroskedasticity . Specifically, we observed that higher NDR scores were associated with increased standard deviations (SD) (Supplementary Figure 2a), violating the assumption of constant dispersion required for standard t-tests. This violation adversely impacted performance, leading to lower positive predictive value (PPV) and reduced sensitivity across different significance thresholds (Table 2). Additionally, the differences in NDR scores across comparison groups generally followed a normal distribution ( Supplementary Figure 2b). Drawing inspiration from RNA-Seq methodologies, such as DESeq2 32, a shrinkage approach for variance estimation was implemented to address these challenges. The proposed method, termed NDRDiff, incorporates a Z -test with empirical Bayes shrinkage of the SD to adjust for heteroskedasticity. Initially, SD values were estimated for each gene within the respective sample groups. Observations revealed that the log-transformed SD (logSD) of NDR scores approximated a normal distribution (Supplementary Figure 2c). LOESS regression was performed to model the relationship between the logSD and the mean NDR score. Subsequently, empirical Bayes shrinkage was applied to the logSD values using mean and variance from .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 20 LOESS regression as prior mean and variance respectively, yielding posterior estimates of the logSD: 𝜇𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 = 𝜇𝑝𝑟𝑖𝑜𝑟 𝜎𝑝𝑟𝑖𝑜𝑟 2 + 𝜇𝑜𝑏𝑠 𝜎𝑜𝑏𝑠 2 1 𝜎𝑝𝑟𝑖𝑜𝑟 2 + 1 𝜎𝑜𝑏𝑠 2 Equation 2 where: • 𝜇𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 is the posterior (shrunk) estimate of the logSD. • 𝜇𝑝𝑟𝑖𝑜𝑟 and 𝜇𝑜𝑏𝑠 are the prior and observed means of logSD, respectively. • 𝜎 𝑝𝑟𝑖𝑜𝑟 2 and 𝜎 𝑜𝑏𝑠 2 are the corresponding variances. The shrunk SDs were derived from the posterior estimates of logSD. For genes in which the shrunk SD deviated by more than two SDs above the LOESS regression fit, the original unshrunk SD was retained to avoid false positives due to over -shrinkage. These shrunk SDs were subsequently used in a two -sample Z-test to identify genes with significant differences in NDR scores between sample groups. Tissue-Specific Gene Set Selection Tissue-specific gene sets were identified by comparing gene expression profiles from GTEx Whole Blood and TCGA Colorectal RNA -Seq datasets. Genes were classified as blood-specific if their TPM values were ≥5 in blood and <0.2 in CRCs, and as CRC- specific if their TPM values were ≥ 5 in CRC and <0. 01 in blood. Both datasets were sourced from the UCSC Toil RNA -Seq Recompute Compendium33 and accessed via the Xena Functional Genomics Explorer platform ( https://xenabrowser.net). The datasets were processed using a standardized workflow to eliminate computational batch effects, ensuring consistency in the analysis. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 21 Gene Set Enrichment Analysis Gene Set Enrichment Analysis (GSEA) was performed to explore enriched pathways. GSEA was performed using the fgsea 34 package with Z-statistics as input, analyzing gene sets from MSigDB Hallmarks28 and HPA Tissue mRNA datasets29. Statistics Statistical analyses were conducted as follows: group comparisons were performed using the non-parametric Wilcoxon test, which does not assume an underlying data distribution. Correlation analyses were carried out using Pearson correlation when a linear relationship was expected; otherwise, Spearman correlation was applied. For differential NDR score comparisons across sample groups, standard t-tests and Z- tests incorporating dispersion shrinkage (NDRDiff) were employed. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 22 Supplementary Materials Supplementary Figure 1 Supplementary Figure 1 Characterization of GC-bias corrected NDR scores. (a, b) Correlation of original (a) and GC-bias corrected (b) NDR scores of healthy cfDNA samples with log2 gene expression from GTEX Whole Blood RNA-Seq. (c, d) Comparison of mean NDR score (y-axis) of Blood and CRC-specific gene sets (x-axis) for CRC (left) and healthy samples (right) for original (c) and GC-bias corrected NDR scores (d). .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 23 Supplementary Figure 2 Supplementary Figure 2 Assessment of NDR score characteristics and distribution for differential comparisons. (a) Scatter plots illustrating the distribution of the standard deviation (SD) of NDR scores relative to mean NDR scores. (b) Distribution of differences in NDR scores between two sample groups (NDR Diff), with a fitted normal distribution (red line). (c) Quantile–quantile plot of log10 -transformed SD values. (d) Scatter plot showing the relationship between mean NDR scores and SD of NDR scores, with a LOESS regression curve (red line). The blue and light blue dots represent post-shrinkage SD and original SD of transcripts, respectively, with the blue dots indicating the final SD used for the Z-test. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 24 Supplementary Figure 3 Supplementary Figure 3 Further validation of NDRDiff results from the comparison of healthy and mCRC samples. (a) Enrichment plot for the colon gene set from the HPA Tissue Gene Expression dataset. (a) log2(TPM + 0.01) values for DNS genes detected in the GTEx Whole Blood RNA-Seq dataset (left) and TCGA COAD-READ RNA-Seq data (right). (b) Comparison of mean NDR scores between healthy and CRC samples for selected DNS genes (CLDN4, IRAG2, and BIN2) . (c) Comparison of mean NDR scores for identified DNS genes that are downregulated (top) and upregulated (bottom) in CRC samples, along with their associatio n with ichorCNA tumor fraction estimates. Healthy samples (assumed to have a tumor fraction of zero) are indicated in red, while cancer samples are shown in green. Pearson correlation was used to assess correlation between mean NDR score and tumor fraction estimates. (d) GSEA results showing the normalized enrichment scores (NES) for the significant gene sets from the MSigDB Hallmark database. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 25 Acknowledgment Public RNA -Seq from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) were accessed from the UCSC Toil RNA -seq Recompute Compendium33. TCGA RNA-Seq data used were generated by the TCGA Research Network: https://www.cancer.gov/tcga. The GTEx Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 26

1. Colorectal cancer. https://www.who.int/news-room/fact-sheets/detail/colorectal-cancer. 2. Siegel, R. L., Wagle, N. S., Cercek, A., Smith, R. A. & Jemal, A. Colorectal cancer statistics, 2023. CA: A Cancer Journal for Clinicians 73, 233–254 (2023). 3. Ilyas, M. I. M. Epidemiology of Stage IV Colorectal Cancer: Trends in the Incidence, Prevalence, Age Distribution, and Impact on Life Span. Clin Colon Rectal Surg 37, 57–61 (2023). 4. Colorectal cancer statistics, 2023 - PubMed. https://pubmed.ncbi.nlm.nih.gov/36856579/. 5. Cancer of the Colon and Rectum - Cancer Stat Facts. SEER https://seer.cancer.gov/statfacts/html/colorect.html. 6. Mannelli, C. Tissue vs Liquid Biopsies for Cancer Detection: Ethical Issues. J Bioeth Inq 16, 551–557 (2019). 7. Loy, C., Ahmann, L., Vlaminck, I. D. & Gu, W. Liquid Biopsy Based on Cell-Free DNA and RNA. (2024) doi:10.1146/annurev-bioeng-110222-111259. 8. Lui, Y. Y. N. et al. Predominant hematopoietic origin of cell-free DNA in plasma and serum after sex-mismatched bone marrow transplantation. Clin Chem 48, 421–427 (2002). 9. Leon, S. A., Shapiro, B., Sklaroff, D. M. & Yaros, M. J. Free DNA in the serum of cancer patients and the effect of therapy. Cancer Res 37, 646–650 (1977). 10. Sorenson, G. D. et al. Soluble normal and mutated DNA sequences from single-copy genes in human blood. Cancer Epidemiol Biomarkers Prev 3, 67–71 (1994). 11. Katsman, E. et al. Detecting cell-of-origin and cancer-specific methylation features of cell-free DNA from Nanopore sequencing. Genome Biol 23, 158 (2022). 12. Shen, S. Y. et al. Sensitive tumour detection and classification using plasma cell-free DNA methylomes. Nature 563, 579–583 (2018). 13. Phallen, J. et al. Direct detection of early-stage cancers using circulating tumor DNA. Science Translational Medicine 9, eaan2415 (2017). .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 27 14. Strickler, J. H. et al. Genomic Landscape of Cell-Free DNA in Patients with Colorectal Cancer. Cancer Discovery 8, 164–173 (2018). 15. Stackpole, M. L. et al. Cost-effective methylome sequencing of cell-free DNA for accurately detecting and locating cancer. Nat Commun 13, 5566 (2022). 16. Sadeh, R. et al. ChIP-seq of plasma cell-free nucleosomes identifies gene expression programs of the cells of origin. Nat Biotechnol 39, 586–598 (2021). 17. Ivanov, M., Baranova, A., Butler, T., Spellman, P. & Mileyko, V. Non-random fragmentation patterns in circulating cell-free DNA reflect epigenetic regulation. BMC Genomics 16 Suppl 13, S1 (2015). 18. Snyder, M. W., Kircher, M., Hill, A. J., Daza, R. M. & Shendure, J. Cell-free DNA Comprises an In Vivo Nucleosome Footprint that Informs Its Tissues-Of-Origin. Cell 164, 57–68 (2016). 19. Ulz, P. et al. Inferring expressed genes by whole-genome sequencing of plasma DNA. Nat Genet 48, 1273–1278 (2016). 20. Zhu, G. et al. Tissue-specific cell-free DNA degradation quantifies circulating tumor DNA burden. Nat Commun 12, 2229 (2021). 21. Moss, J. et al. Comprehensive human cell-type methylation atlas reveals origins of circulating cell-free DNA in health and disease. Nat Commun 9, 5068 (2018). 22. Benjamini, Y. & Speed, T. P. Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Res 40, e72 (2012). 23. Adalsteinsson, V. A. et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nat Commun 8, 1324 (2017). 24. Doebley, A. L. et al. A framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. Nature communications 13, 7475 (2022). 25. Fujiwara-Tani, R. et al. Anti-claudin-4 extracellular domain antibody enhances the antitumoral effects of chemotherapeutic and antibody drugs in colorectal cancer. Oncotarget 9, 37367–37378 (2018). .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint 28 26. Prüschenk, S. & Schlossmann, J. Function of IRAG2 Is Modulated by NO/cGMP in Murine Platelets. Int J Mol Sci 23, 6695 (2022). 27. Volz, J. et al. BIN2 orchestrates platelet calcium signaling in thrombosis and thrombo-inflammation. J Clin Invest 130, 6064–6079 (2020). 28. Liberzon, A. et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 1, 417–425 (2015). 29. Uhlén, M. et al. Proteomics. Tissue-based map of the human proteome. Science 347, 1260419 (2015). 30. Li, S. et al. Comprehensive tissue deconvolution of cell-free DNA by deep learning for disease diagnosis and monitoring. Proceedings of the National Academy of Sciences 120, e2305236120 (2023). 31. Morales, J. et al. A joint NCBI and EMBL-EBI transcript set for clinical genomics and research. Nature 604, 310–315 (2022). 32. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15, 550 (2014). 33. Vivian, J. et al. Toil enables reproducible, open source, big biomedical data analyses. Nat Biotechnol 35, 314–316 (2017). 34. An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation | bioRxiv. https://www.biorxiv.org/content/10.1101/060012v1. .CC-BY-NC-ND 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 18, 2025. ; https://doi.org/10.1101/2025.04.17.649313doi: bioRxiv preprint