Integrated fragmentomic profile and 5-Hydroxymethylcytosine of capture-based low-pass sequencing data enables pan-cancer detection via cfDNA

Sample collection and study design

In total, low-pass (mean mapped reads: 38.32 ± 8.02 million) 5hmC-seq data from 85 healthy controls and 311 cancer patients were included in this study. We obtained 5hmC sequencing data from 33 LUAD (18), 74 PDAC (37), 132 HCC (38), and 85 control (37) samples from three publications, and 72 GBM samples from an unpublished paper.

Plasma samples from 15 healthy controls and 59 cancer patients were collected and subjected to WGS (mean: 14×, range: 8.6-23.7×) to establish the preferred end marker set, while the low-pass WGS data from 19 healthy controls and 98 cancer patients (mean: 3×, range: 2.3-6.2×) were obtained to cross-check with low-pass 5hmC sequencing data. 13 of the samples with low-pass WGS sequencing had above 10x WGS sequencing data. The demographics and clinical characteristics of the cohort are summarized in Supplementary Tables 1 and 2. The design of the study is shown in Figure 1.

• Download figure
• Open in new tab

Figure 1. Overview of the study design.

cfDNA from cancer patients and controls were sequenced with WGS. Available 5hmC sequencing data from the same four cancer types and controls were combined with WGS data to get the accurate fragmentomic information and upon it a machine learning model was built to distinguish healthy controls and cancer patients.

Plasma sample collection and cfDNA WGS sequencing

We collected plasma samples from 46 GBM patients, 29 PDAC patients and 34 controls from West China Hospital. Plasma samples of 37 HCC patients and 32 LUAD patients were collected from Mengchao Hepatobiliary Hospital of Fujian Medical University and the First Affiliated Hospital of the Second Military Medical University, respectively. For every subject, the QIAseq cfDNA All-in-one Kit (Qiagen) was applied on the 2 ml plasma sample for the cfDNA extraction and library construction. Pair-end 150 bp sequencing of the libraries was performed on the Illumina Novaseq 6000 platform.

cfDNA 5hmC profiling and sequencing

cfDNA extraction and 5hmC library construction was performed as previously described (28). Firstly, cfDNA was extracted from 2 ml plasma sample using QIAamp Circulating Nucleic Acid Kit (QIAGEN Inc., Valencia, CA, USA). Then, cfDNA (5–10 ng) ligated with sequencing adaptors was incubated in a 25 μL reaction solution containing HEPES buffer (50 mM, pH 8.0), 25 mM MgCl2, 60 μM N3-UDP-Glc (ActiveMotif, Carlsbad, CA, USA), and 12.5 U β-glucosyltransferase (NEB, Beverly, MA, USA) for 2 h at 37 °C. Next, 2.5 μL DBCO-PEG4-biotin (Sigma, Carlsbad, CA, USA) was directly added and incubated for 2 h at 37 °C. 10 μg sheared salmon sperm DNA (Life Technologies, USA) was added, then the Micro Bio-Spin 30 Column (Bio-Rad, Hercules, CA, USA) was used to purify the DNA following the instruction, and the final volume was adjusted to 25μL. After that, the purified DNA was incubated with 5 μL C1 streptavidin beads (Life Technologies, USA) in buffer 1 (5 mM Tris pH 7.5, 0.5 mM EDTA, 1 M NaCl and 0.2% Tween 20) for 30 min. The beads were subsequently undergone three 5-min washes each with buffer 1, buffer 2 (buffer 1 without NaCl), buffer 3 (buffer 1 with pH 9) and buffer 4 (buffer 3 without NaCl). Then, the beads were resuspended in water and amplified with 11 cycles of PCR amplification (initial denaturing at 98 °C for 45 s, followed by 11 cycles of denaturing at 98 °C for 15 s, annealing at 60 °C for 30 s, extension at 72 °C for 30 s, and a final extension at 17 °C for 5 min). The amplified products were purified using 0.8 × AMPure XP beads (Beck-man Coulter, Fullerton, CA, USA). Pair-end 150 bp sequencing was performed on the Illumina Novaseq 6000 platform.

Fragmentomic profiling of cfDNA and cancer prediction analysis

The details of the preferred end, 5hmC signatures, size profile and coverage profile calculation are described in the Supplementary Material. We used a random forest model to distinguish healthy people from cancer patients using fragmentomic features. To estimate the prediction error, we used five cross-validations. Near-zero variance features were removed. The training, validation, and test sets account for 60%, 20%, and 20% of the data, respectively. The samples were selected randomly in a balanced way to keep the ratio of the number of cancer to non-cancer samples similar in the training, validation, and test subsets. Feature importance was calculated using the training data in each cross-validation run, and we sorted the features according to the mean value of feature importance. The number of features used in the final model was obtained by the highest AUC value in the validation set. To build an integrative model, we fitted the features selected by every feature-alone model into one random forest model, and performed feature selection, training, testing as described above. Random forest machine learning was implemented using the python package sklearn.ensemble.RandomForestClassifier with the following parameters: n_estimators = 100, criterion=“gini”.

Statistics

All statistical analyses were performed using R version 3.6.1 and python 3.7.0. All tests were two-sided, and P values < 0.05 were considered statistically significant. The ROC-AUC plot and AUC value were implemented using sklearn.metrics.plot_roc_curve.

Ethics

This study was approved by the Ethics Committee of Sichuan Cancer Hospital (SCCHEC-02-2016-005). The written informed consent was obtained from all participants.

Constructing the cfDNA preferred end map

cfDNA fragments ending at specific genomic coordinates that were statistically more abundant than the Poisson distribution was an essential fragmentation signature referred to as the preferred end (10). We utilized WGS data from the four cancer types to build a preferred end map and checked the feature using low-pass data. The basis of utilizing cfDNA preferred end coordinates to indicate the occurrence and location of cancer was its ability to locate nucleosomes, the position of which is related to cell identity (19). To confirm this, we defined upstream (U) ends and downstream (D) ends by their genomic coordinate order on the reference genome calculated via alignment. We surveyed the distance distribution between adjacent U ends, D ends, and nucleosome centers, and the result suggested the preferred ends were enriched on the DNA linker (Supplementary Figure 1a-f). We further checked the distribution of the predicted nucleosome location based on non-malignant tissues at chr12: 34517269-34519122, and the preferred U and D ends of WGS data from controls (n=15) and HCC patients (n=15). As expected, the preferred U and D ends of the controls were enriched at the midpoint between two nucleosomes, which is the location of the DNA linker (Figure 2a), and U ends were found downstream of their nearest D ends. However, in HCC patients, U ends were upstream of their nearest D ends, and only D ends were found at one midpoint between two nucleosomes (Figure 2b), which indicated inconsistent nucleosome positions between HCC patients with controls. We observed the same abnormal preferred end distribution in GBM, LUAD and PDAC patients (Supplementary Figure 1g-i).

• Download figure
• Open in new tab

Figure 2. Genome-wide cfDNA fragment preferred end map construction.

cfDNA preferred end signals in a nucleosome array region (chr12: 34517269-34519122) in healthy (a) and HCC (b) subjects. Brown dots at the bottom represent the predicted nucleosome center loci reported in a previous study (8, 22). (c) Venn diagram showing the intersection of the preferred U and D ends between the cancer cohort and the control cohort. (d) Circle diagram showing the density of preferred ends in 1 Mb windows in four types of cancer and in healthy controls (the order of sample sources from outer ring to inner ring are control, GBM, HCC, LUAD and PDAC). (e) Enrichment of preferred U/D ends in promoter regions. (f) Enrichment of preferred U/D ends in health-specific open chromatin regions. (g) Boxplots showing the ratios of the cancer-specific to health-specific preferred U/D ends in low-pass WGS data for the four cancer types (Wilcoxon rank-sum test).

By using WGS data, we constructed cfDNA fragment preferred end maps for the four cancer types and control (Figure 2c-d). Preferred end enrichment was observed in intron and open chromatin regions (Supplementary Figure 2). Additionally, we noticed that for every cancer type, both U and D preferred ends were more enriched in health-specific open chromatin regions than in cancer-specific open chromatin regions (Supplementary Figure 2). One explanation was the nucleosome pattern in the health-specific open chromatin regions was unstable, as the enrichment of both U and D preferred ends in healthy cfDNA were significantly lower than cancer patients (Mann‒Whitney, P<0.05). Another was the heterogeneous nature of cfDNA, as cfDNA can be derived from various tissue sources. To search for preferred ends specific to each of the four cancer types, we determined genome-wide cancer-specific preferred U/D end coordinates and health-specific preferred U/D end coordinates based on the Youden index. In total, we got 81,083 D and 80,505 U preferred ends. Intriguingly, compared with all repeatable ends (preferred ends appeared in at least two samples in each cancer type), cancer-specific preferred U/D ends were no longer enriched in health-specific open chromatin regions, but they showed enrichment in promoter regions and cancer-specific open chromatin regions (Figure 2e-f, Supplementary Figure 3-4). Furthermore, health-specific preferred U/D ends showed enrichment in intergenic regions and were not found in promoters, as cancer-specific preferred ends were (Supplementary Figure 4). To cross check the cancer specificity of the preferred ends in low-pass WGS data, we calculated a ratio, which was the ratio of detected cancer-specific preferred end numbers to health-specific preferred end numbers for each sample. The results showed significantly higher ratios of cancer-specific to health-specific preferred U/D ends in the three cancer types (except GBM) than in the control (Figure 2g, Mann‒Whitney, p value<0.001). In addition to the difference between cancer and health, statistically significant differences in the ratios between cancer types were observed (Supplementary Figure 5).

Fragmentomic features in 5hmC data

We checked the batch effect of the 5hmC sequencing data based on size profile and coverage profile, and found no clustering of samples according to cancer types (Supplementary Figure 6).

The cfDNA captured by 5hmC marks comprised a subset of DNA in whole plasma. In theory, one sequenced fragment should contain at least one 5hmC. Hence, we investigated the consistency of the size profile between WGS data and 5hmC sequencing data. According to the results, fragment sizes of 5hmC cfDNA in the four cancer samples were significantly shorter than those in the control (Mann‒Whitney test, p value<0.001, pan-cancer Median N50 = 174 bp, control Median N50 = 210 bp), though the difference between cancer and control groups in WGS data was only observed for HCC and LUAD (Figure 3a). Next, we confirmed that the cfDNA size profile of the 5hmC sequencing data was consistent with that of the WGS data in the size range of 0-220 bp (Figure 3b), with both data types displaying a dominant peak at 167 bp (Supplementary Figure 7a). A difference between the cancer and control samples in both short (100-150 bp) and long (151-220 bp) cfDNA fragments (7) was found in the 5hmC data, as in the WGS data (Supplementary Figure 7b-c), though the 5hmC data showed a unique secondary peak at ∼320 bp (Figure 3c, Supplementary Figure 7a). In the WGS data, cfDNA fragments ranging from 0-220 bp accounted for nearly 100% of fragments ranging from 0-800 bp, but the percent value of the same size interval in the 5hmC sequencing data ranged from 31.8% to 99.6% (Figure 3f). Thus, we confirmed the presence of ultra-long fragments (221-500 bp) in the 5hmC cfDNA. We suspected that the existence of ultra-long fragments may have resulted from the capture-based technique applied in 5hmC sequencing, and we consequently assessed the relationship between fragment length and 5hmC peak number. In ultra-long fragments, the percentage of fragments with more than one peak was significantly higher than that in non-ultra-long fragments (0-220 bp) (Mann‒Whitney test, p value<0.001) (Figure 3d). Indeed, enrichment analysis of the ultra-long fragments with more than one peak revealed primary enrichment at CpG islands (Supplementary Figure 7d). This indicated that the capture-based sequencing technique is more likely than WGS to capture ultra-long cfDNA fragments. Another finding was that in all cancer types, the percentage of ultra-long fragments was significantly lower than that in the control (Supplementary Figure 7b). We examined the percentage of ultra-long fragments with more than one 5hmC peak in the four cancer types and the control and found that the value was significantly lower in the former (Figure 3e). As reported by a recent study, the stability of circulating DNA derives mostly from the nucleosome structure (21). We hence evaluated the proportion of 5hmC cfDNA fragments at 146 bp, 166 bp, 312 bp, and 332 bp, which indicates the length of one mono-nucleosome and one di-nucleosome (plus the linker size). The result showed a higher proportion of mono-nucleosome-sized cfDNA fragments and a lower proportion of di-nucleosome-sized cfDNA fragments in the cancer group than in the control (Supplementary Figure 7e).

• Download figure
• Open in new tab

Figure 3. cfDNA fragmentomic information in 5hmC sequencing data.

(a) Violin plots showing fragment size profile comparison between cancer samples and controls in terms of N50 (Mann‒Whitney test). (b) tSNE plot showing 5hmC sequencing data and WGS data by fragments below 220 bp. (c) The cfDNA size profile of 5hmC sequencing data. (d) Comparison of the percentage of fragments with more than one 5hmC peak between ultra-long and non-ultra-long fragments in every sample. The dashed line indicates the diagonal line. (e) Violin plots showing the percentage of ultra-long fragments with more than one 5hmC peak in cancer samples and controls (Mann‒Whitney test). (f) Boxplots showing the percentage of fragments ranging from 0-220 bp and 221-500 bp among fragments ranging from 0-800 bp in WGS data and 5hmC sequencing data (Wilcoxon rank-sum test). (g) Boxplots showing the ratios of the cancer-specific to health-specific preferred U/D ends in each of the four cancer types (Wilcoxon rank-sum test).

We also explored the cancer specificity of the pan-cancer preferred end set built upon WGS data in the 5hmC sequencing data by calculating the ratio of cancer-specific to health-specific preferred U/D ends. The ratios were significantly increased in cfDNA samples in three cancer types (except PDAC) compared with healthy individuals (Figure 3g, Supplementary Figure 5c-d, Mann‒Whitney, p value<0.005). Coverage profiles in 5hmC sequencing data also showed a higher variance than WGS data in samples from control and cancer (Supplementary Figure 8). We explored the coverage patterns of short, long, and ultra-long cfDNA fragments at the 5 Mb bins of the genome. Similar to WGS data (Supplementary Material), cancer patients exhibited multiple unstable genomic regions: coverage of short/long/ultra-long fragments was inconsistent in cancer patients but consistent in controls (Figure 4a, Supplementary Figure 9a). We analyzed the genome-wide coverage profile correlation of cancer patients to controls in the 5hmC data and found that the coverage profiles were consistent among controls and that the correlation value in cancer patients was significantly lower (Wilcoxon rank-sum test, p value<0.005). Importantly, the largest difference between cancer samples and controls was found with respect to the coverage profile of ultra-long fragments (Wilcoxon rank-sum test, p value<0.005) (Figure 4b).

Genomic distribution of 5hmC in cancer and control cohorts

To gain an understanding of the genomic regions associated with hydroxymethylation in cfDNA, we first determined 5hmC-enriched loci in the control and four types of cancer, which were detected as peaks via MACS2 (22). 5hmC-modified regions among samples were compared in 1 kb bins on the reference genome, and thus a consensus list of the absence and presence of 5hmC-modified peaks among the samples was obtained. 5hmC signatures were enriched over genic features, most significantly in the promoter, 5’UTRs, 3’UTRs, exons, and transcription end sites (TESs) (Figure 5a). Comparison of 5hmC signature enrichment in the four types of cancers with the control revealed significant differences, whereby 5hmC peak enrichment was lower in all cancer types than in the control over promoters, 5’UTRs, exons, and CpG islands (Wilcoxon rank-sum test, p value<0.001) (Figure 5a).

Differential analysis of 5hmC-modified peaks between cancer patients and controls by Fisher’s exact test detected 1,010, 6,395 and 773 differentially modified 5hmC peaks in GBM, HCC and PDAC, respectively. Moreover, this pan-cancer 5hmC peak set was able to separate most cancer samples from the control (Figure 5b). KEGG enrichment analysis of genes located within the differentially modified 5hmC peaks revealed associated oncogenesis pathways for each cancer type (Supplementary Figure 10), such as the neurotrophin signaling pathway and platelet activation for GBM (23, 24), the MAPK signaling pathway and focal adhesion for HCC (25, 26), and the FoxO signaling pathway and insulin signaling pathway for PDAC (27).

Cancer detection by combining fragmentomic features and 5hmC signatures using 5hmC sequencing data

As 5hmC sequencing data retain both 5hmC signatures and fragmentomic information, it is expected to theoretically show high sensitivity and specificity in cancer detection. A receiver operator characteristic (ROC) curve was used to evaluate the performance of the classifier. In total, 53 5hmC signatures were selected, and the AUC of the classifier was 0.876 in the validation set (sensitivity = 82.26%, specificity = 82.35%) and 0.872 in the test set (sensitivity = 81.97%, specificity = 82.35%) (Figure 5d-e, Supplementary Figure 11 c).

Overall, the fragmentomic information in the 5hmC sequencing data performed well. Using the size profile as features, 40 were selected to build a pan-cancer prediction model, among which 32 features were of the ultra-long and longer size range (>= 220 bp) (Figure 5c). The AUC value was 0.981 in the validation set (sensitivity = 93.55%, specificity = 94.12%) and 0.882 in the test set (sensitivity = 71.43%, specificity = 88.24%) (Figure 5d-e). Using preferred ends as features, 37 were selected to build a pan-cancer prediction model. ROC analysis showed an AUC value of 0.940 in the validation set (sensitivity = 83.87%, specificity = 94.12%) and 0.899 in the test set (sensitivity = 73.02%, specificity = 94.12%) (Figure 5d-e, Supplementary Figure 11a). Regarding the coverage profile of the 5hmC sequencing data, we explored the genomic distribution of short, long, and ultra-long cfDNA fragments on a 100 kb window of the genome. This model achieved an AUC value of 0.946 in the validation set (sensitivity = 80.65%, specificity = 94.12%) and 0.882 in the test set (sensitivity = 71.43%, specificity = 88.24%) (Figure 5d-e).

Finally, we constructed an integrated cancer screening model by combining fragmentomic features and 5hmC signatures in 5hmC data. Sixty-three features, including 10 5hmC signatures, 24 size profile features, 21 coverage profile features, and 8 preferred end features, were selected to build the pan-cancer prediction model. The AUC value was 0.927 in the validation set (sensitivity = 93.44%, specificity = 88.24%) and 0.920 in the test set (sensitivity = 88.52%, specificity = 82.35%) (Figure 5d-e). The performance of single-cancer detection is depicted in Supplementary Figure 11b. With data from the four cancer types, we further explored tissue-of-origin prediction. The performance of the random forest model with 5hmC fragmentomic features and the integrated model are shown in Supplementary Tables 3-5.