Global and genomic distribution of 5hmC modifications

We evaluated the global 5hmC level variation in cancer by using ultra-sensitive capillary electrophoresis-electrospray ionization-mass spectrometry (CE-ESI-MS)³⁰. Global 5hmC levels of the tumor gDNA were markedly decreased compared to the adjacent healthy tissue gDNA, with an average of 85% and 64% reduction in colorectal and gastric tumor samples, respectively. Global 5hmC levels of the cancer patients' plasma cfDNA showed a more limited decrease compared to control plasma cfDNA, consistent with low proportions of tumor-derived DNA in the total cfDNA pool (Supplementary information, Figure S2).

In plasma cfDNA, 5hmC was enriched within gene bodies and DNase I sensitive peaks whereas it was depleted at transcription start sites, CpG islands and transcription factor (TF) binding peaks relative to the flanking areas (Supplementary information, Figure S3A-S3F). This suggests an accumulation of 5hmC surrounding TFs at active transcription sites. 5hmC was also enriched in regions bearing permissive histone marks such as H3K27ac, H3K4me1 and H3K9me1, whereas it was underrepresented in regions bearing repressive markers such as H3K9me3 (Supplementary information, Figure S3G-S3R). The genomic enrichment pattern of 5hmC was consistently observed in tissues gDNA (Supplementary information, Figure S3) and was similar between disease and healthy samples (Supplementary information, Figure S4).

Differential 5hmC loci associated with colorectal cancer

The average 5hmC profiles of plasma cfDNA were distinct from those of tissue and blood cell gDNA (Figure 2A), which could be due to their distinct cell origins and/or the different DNA degradation properties in cell-free circulation. Among gDNA profiles, variations attributable to tissue identity (colon and stomach tissues, white blood cells) were dominant over variations attributable to disease status (healthy individual vs cancer patient, tumor vs adjacent tissue). Across plasma cfDNA profiles, colorectal and gastric cancer samples were more closely related with each other than with healthy controls (Figure 2A).

Open in a separate window

Figure 2

Differential 5hmC loci associated with cancer. (A) Average 5hmC levels in gene bodies in healthy controls (health) and cancer patients (colon, stomach), estimated for plasma cfDNA (plasma cf), white blood cell genomic DNA (WBC) and tissue genomic DNA (tumor, adjacent), were clustered by correlation distance. (B) Counts per million reads at SULF1 gene (plus ± 20 kb region) in plasma cfDNA of the 15 healthy controls and 18 colorectal cancer patients. The moving averages at 0.01 smoother span are shown. (C) The distribution of colorectal cancer-associated 5hmC loci detected at 5% false discovery rate in plasma cfDNA. Each vertical bar denotes a differential locus (a histone modification peak or a gene body). The color key indicates the relative magnitude of log₂ fold change in cancer patients vs controls. (D) Pearson's correlation of log₂ fold changes between all analyzed genes and their neighboring genes (points) was plotted against the null distribution of correlation with their first neighboring genes (curves), generated by shuffling gene positions for 1 000 times. Blue and orange points denote data from plasma cfDNA and tissue gDNA, respectively, for colorectal cancer. In C and D, chromosome 1 is shown as an example. (E) Cancer plasma cfDNA and tumor gDNA exhibit correlation in average 5hmC density (library size and feature length normalized log₂ counts, black bars). However, there is no correlation in the log₂ fold changes between differential 5hmC loci detected (between cancer vs health (plasma cfDNA)) and tumor vs adjacent tissue (tissue gDNA), (orange bars). (F) Genes with a 5hmC level elevated in cancer plasma cfDNA (cancer cf) are enriched in genes with high 5hmC level in tissue gDNA (tumor high, adjacent high). To estimate fold enrichment, the 1st, 5th and 10th percentile genes in descending order of the log₂ fold change in cancer cfDNA were compared against the corresponding percentile genes in descending order of the average 5hmC level in tissue gDNA. Similarly, genes with a 5hmC level decreased in cancer plasma, and cfDNA are enriched in genes with low 5hmC levels in tissue gDNA (tumor low, adjacent low). In contrast, differentially marked genes detected in tumor gDNA (tumor) show no such enrichment pattern. Dashed line denotes no enrichment.

We compared 5hmC profiles from plasma cfDNA between 15 colon cancer patients and 18 healthy controls in the discovery batch to identify differentially modified 5hmC loci. The profiles were separated into 18 feature categories: gene bodies, promoters, CpG islands and cis-regulatory elements delineated by the Encyclopedia of DNA Elements (ENCODE)³¹. A parallel analysis compared 5hmC profiles from gDNA between colorectal tumors and adjacent tissues in 30 patients in the tissue discovery batch. All feature categories showed enrichment of differentially modified 5hmC loci (Supplementary information, Table S4). Figure 2B shows a differential locus detected in plasma cfDNA at the SULF1 (sulfatase 1) gene. In cancer patients' plasma cfDNA, the 5hmC levels in SULF1 were elevated in both exons and introns, with a peak pattern similar to that of tissue gDNA (Supplementary information, Figure S5). Neighboring differential 5hmC loci showed regionally elevated or decreased 5hmC levels (Figure 2C). For example, the correlation of cancer-associated 5hmC changes between neighboring genes was significantly higher than a null distribution generated by shuffling gene positions within a chromosome (Figure 2D). This may suggest that 5hmC modifications occur and change in a relatively long-range, region-wise pattern.

The 5hmC densities across genomic features showed a correlation between cancer plasma cfDNA and tumor gDNA samples (Spearman's ρ 0.34-0.84, Figure 2E). This locus-specific correlation of the 5hmC density is expected because of biological constraints. In contrast, we found no correlation between cancer patients' plasma cfDNA and tumor gDNA for cancer-associated 5hmC changes (Figure 2E). This can be explained by the greater variation in 5hmC levels among different tissues than that between disease status (Figure 2A). When gDNA from tumor tissue is released into plasma and mixed with the vast amount of background cfDNA derived from a variety of different tissues, the additional tumor signal observed at a given locus would be determined by the order of locus-, tissue- and disease-specific variations. Consistent with this expectation, we found that genes with 5hmC level elevated in cancer patients' plasma cfDNA were enriched in genes with high 5hmC levels in the tumor tissue gDNA. The top 1% of genes with the most elevated 5hmC levels in cancer patients' plasma cfDNA were enriched by over five-fold in the top 1% of genes with the greatest 5hmC levels in tumor and adjacent tissues (Fisher's exact tests P < 0.001). Similarly, genes with 5hmC levels decreased in cancer patients' plasma cfDNA were enriched in genes with low 5hmC levels in tumor and adjacent tissues (Figure 2F). In contrast, no such enrichment pattern was observed for the differentially modified 5hmC loci detected in tumor gDNA (Figure 2F).

Classification of colorectal cancer by 5hmC markers derived from plasma cfDNA

Unsupervised hierarchical clustering using differentially modified 5hmC loci derived from plasma cfDNA generally separated colorectal cancer patients from healthy individuals in the validation batch (Figure 3A). Across the various feature categories, the log₂ fold change of 5hmC levels in gene bodies showed the greatest correlation between the discovery and validation batches (Spearman's ρ = 0.79, Figure 3B), indicating that 5hmC loci in gene bodies are potentially more stable cancer biomarkers. We selected 989 differential loci in gene bodies detected at a 5% false discovery rate (FDR) with a greater than 1.2-fold change (increase or decrease in cancer; Supplementary information, Table S5) for cancer classification. A model-based classifier that applies elastic net regularization on logistic regression was trained using the discovery samples (15 patients vs 18 controls) and then tested in the validation samples (24 patients vs 35 controls). Receiver operating characteristic (ROC) curves were generated to evaluate the performance using the area under the curve (AUC). The prediction algorithm achieved 83% sensitivity and 94% specificity (AUC = 0.95, Figure 3C) for patient classification.

Open in a separate window

Figure 3

Performance of 5hmC biomarkers for colorectal cancer. (A) The heat map shows clustering of cfDNA samples from both the discovery and validation batches, using the 989 differential gene bodies detected in plasma cfDNA from the discovery batch. (B) Correlation of 5hmC changes in cancer between the discovery and validation batches is higher in plasma cfDNA (cancer patients vs healthy individuals) than in tissue gDNA (tumors vs adjacent tissues), especially for 5hmC in gene bodies. (C, D) Classification of two independent validation batches using 5hmC classifier derived from plasma cfDNA from the discovery batch. (E) Classification of an independent set of colon cancer tumor tissues using 5hmC biomarkers detected from the discovery batch of tissue samples (tumors vs adjacent tissues). (F) The predicted cancer probability (i.e., score) based on 5hmC classifier from plasma cfDNA shows a significant trend associated with clinical stage. Patients after surgery show predicted scores undistinguishable from healthy individuals. (G) The 5hmC cfDNA classifier for colorectal cancer is disease- and potentially cancer type-specific, showing decreasing predicted probability in cfDNA from stomach, liver, pancreatic and thyroid cancer patients. AUC, area under curve; CAC, cancer patients; HEA, healthy controls; NOR, patients with benign tumor.

An additional validation batch (32 patients vs 37 controls) was independently collected and tested, achieving 88% sensitivity and 89% specificity (AUC = 0.94, Figure 3D). The classifier was further tested on a set of US samples (Supplementary information, Table S2), all of European descent, collected at the University of Chicago Medical Center. The classifier detected 4 out of 5 patients as cancer positive (80% sensitivity) and called 1 out of 6 healthy controls (83% specificity) in this small cohort. Therefore, although the classifier was trained on Chinese patients, it could capture a general signal of 5hmC changes in plasma cfDNA in colorectal cancer.

For comparison, we applied similar approaches to evaluate the performance of 5hmC biomarkers derived from colorectal tumor tissues. We selected 219 differential loci at gene bodies called at 5% FDR and 1.2-fold change between tumor and adjacent tissues (Supplementary information, Table S6) from 30 patients of the discovery batch. The 5hmC tissue biomarkers showed a sensitivity of 86% and a specificity of 100% (AUC = 0.96, Figure 3E) in 14 patients from the tissue validation batch, suggesting that the 5hmC biomarkers from plasma cfDNA exhibit performance comparable to that from tissue gDNA. The colorectal cancer classifiers derived from plasma cfDNA and tissue gDNA profiles are detailed in Supplementary information, Table S7.

Disease sensitivity and specificity of plasma cfDNA-derived 5hmC markers

We next assessed the ability of the 5hmC biomarkers derived from plasma cfDNA to classify cancer stages in a subset of patients with available records. The 5hmC classifier assigned incremental numbers of cancer individuals (predicted cancer probability > 0.5) for patients having undergone surgery for treatment (0/2), patients at cancer stage I (4/6) and patients at stage II and III (38/40) (Cochran-Armitage test for trend P = 0.00038, Figure 3F). The classifier had good but reduced power to call Stage IV patients (12/18), who in general suffered from metastasis to various tissues and are expected to have more complex tumor DNA profiles in circulation due to metastasis.

We further assessed disease specificity of the classifier in patients with colon-related benign diseases (n = 49) and patients with colorectal (n = 71), gastric (n = 61), liver (n = 25), pancreatic (n = 34) and thyroid (n = 46) cancer. Compared to the 86% call rate in colorectal cancer patients, only 8% of patients with benign colon diseases were predicted as having cancer (Figure 3G). The classifier also demonstrated certain tissue specificity, with a decreasing cancer call rate in gastric (85%), liver (44%), pancreatic (29%) and thyroid (28%) cancer patients (Figure 3G). The lower sensitivity in calling the other cancers is not due to any intrinsic difficulty in classifying those cancers, as we achieved a much greater sensitivity in liver and pancreatic cancer using 5hmC markers derived from plasma cfDNA from patients with these respective cancers (data not shown). These results indicate that distantly related cancers can be readily distinguished by the corresponding classifiers through joint testing, while classification of closely related cancers such as colorectal and gastric cancer may be facilitated by additional diagnostic criteria.

A subset of cancer patients had records of classical biomarkers and epidemiological risk factors (Supplementary information, Table S1), with which we compared plasma cfDNA-derived 5hmC biomarkers for cancer detection sensitivity. The detection sensitivity of carcinoembryonic antigen (CEA, 32%), alpha-fetoprotein (AFP, 0%), carbohydrate antigen 125 (CA125, 13%), CA15-3 (0%), CA19-9 (19%), CA72-4 (17%), cytokeratin 19 (49%), neuron-specific enolase (NSE, 21%), and relationship to overweight (body mass index ≥ 25 kg/m², 34%), smoking (9%), alcohol consumption (7%) and previous history of cancer (0%) were all less than 50%. By calling cancer if any conventional biomarker or risk factor is positive, the upper bound detection sensitivity of the combined classical biomarkers and epidemiological factors only reached 54%, a sensitivity rate much lower than the 86% that we could achieve using 5hmC markers. In addition, compared with the methylated SEPT9 (septin 9) gene, a blood-based epigenetic biomarker for colon cancer, our cfDNA 5hmC biomarkers registered a significantly further improved overall sensitivity (0.86 vs 0.48 based on public data)³².

5hmC markers derived from plasma cfDNA in gastric cancer

Next, we analyzed gastric cancer using plasma cfDNA samples. In the discovery batch, 5hmC loci in 7 gastric cancer patients were compared to 18 healthy controls across genomic feature categories (Supplementary information, Table S4). Using the detected differential 5hmC loci, 25 gastric cancer patients could be generally separated from 35 healthy individuals in the validation batch (Figure 4A). Again, 5hmC changes in gene bodies showed relatively higher correlations between the discovery and validation batches compared to other genomic features (Figure 4B). A model-based classifier was generated using the 1 431 differential loci in gene bodies identified at 5% FDR and 1.2-fold change in the discovery batch (Supplementary information, Table S8). This was applied to the validation batch, achieving 92% sensitivity and 91% specificity (AUC = 0.93, Figure 4C). Further assessment of the gastric cancer classifier in an additional, independently collected validation batch (29 patients vs 37 controls) achieved 90% sensitivity and 97% specificity (AUC = 0.97, Figure 4D). The classification performance of the 5hmC biomarkers derived from cancer cfDNA was comparable to that from tumor gDNA samples: a classifier based on 161 differential 5hmC loci in gene bodies detected in 19 pairs of tumors and adjacent tissues in the discovery batch (Supplementary information, Table S9) was applied on 33 pairs of tissues in the validation batch, achieving 82% sensitivity and 94% specificity (AUC = 0.93, Figure 4E). The gastric cancer classifiers derived from plasma cfDNA and tissue gDNA profiles are detailed in Supplementary information, Table S10.

Open in a separate window

Figure 4

Performance of 5hmC biomarkers for gastric cancer. (A) The heat map shows clustering of cfDNA samples from both the discovery and validation batches, using the 1 431 differential gene bodies detected in plasma cfDNA from the discovery batch. (B) Correlation of 5hmC changes in cancer between the discovery and validation batches is higher in plasma cfDNA (cancer patients vs healthy individuals) than in tumor gDNA (tumors vs adjacent tissues), especially for 5hmC in gene bodies. (C, D) Classifying two independent validation batches using 5hmC classifier derived from plasma cfDNA from the discovery batch. (E) Classifying an independent set of gastric cancer tumor tissues using 5hmC biomarkers detected from the discovery batch of tissue samples (tumors vs adjacent tissues). (F) The predicted cancer probability (i.e., score) based on the 5hmC classifier from plasma cfDNA shows a trend associated with clinical stage. The one patient after surgery shows a predicted probability undistinguishable from healthy individuals. (G) The 5hmC cfDNA classifier for gastric cancer is disease- and potentially cancer type-specific, showing decreasing predicted probability in cfDNA from colorectal, liver, pancreatic and thyroid cancer patients.

The 5hmC gastric cancer classifier derived from plasma cfDNA showed a trend of increasing cancer call rate with cancer clinical stage (P = 0.11, Figure 4F). The classifier also demonstrated disease and tissue specificity, with 0% cancer call rate for benign gastric diseases, and with decreasing cancer call rate in patients with colorectal (61%), liver (28%), pancreatic (6%) and thyroid (0%) cancer patients (Figure 4G).

The detection sensitivity of classical biomarkers and epidemiological factors for gastric cancer was 13% (CEA), 6% (AFP), 6% (CA125), 3% (CA15-3), 14% (CA19-9), 29% (CA72-4), 36% (cytokeratin 19), 13% (NSE), 18% (overweight), 25% (smoking), 10% (alcohol) and 7% (previous history of cancer) (Supplementary information, Table S1). The upper bound sensitivity combining these markers and factors was 70%, which again is lower than the 80% call rate observed in our studies.

Tissue origin of the cancer-associated 5hmC changes observed in plasma cfDNA

To demonstrate the tumor relevance of the plasma cfDNA from cancer patients, we sought to examine its source in patient-derived xenograft (PDX) mouse models. PDX mouse models were derived from tumors of three colorectal patients and three gastric patients, each with three independent xenograft animals. Plasma cfDNA of PDX mice was collected at 12-15 weeks of age, from which the 5hmC-containing fragments were enriched and sequenced using the same protocol as with human plasma cfDNA. The proportion of cfDNA derived from the tumor, estimated as the proportion of sequencing reads uniquely mapped to the human genome, was significantly increased in mice grafted with gastric tumors (P = 0.0020) and showed an increased trend in mice grafted with colorectal tumors with fewer passages (P = 0.16, Figure 5A).

Open in a separate window

Figure 5

The origin of cancer-associated 5hmC changes observed in plasma cfDNA. (A) The proportions of human reads in plasma cfDNA captured with 5hmC-Seal are shown for PDX mice grafted with tumor from three gastric cancer patients (stomach_1-3), three colorectal cancer patients (colon_1-3) and for PDX mice without graft (control_1-3). Vertical bars represent s.d. estimated from three replicate PDX mice for each patient. The PDX mice grafted with gastric tumor had greater number of passages (6-10) than those grafted with colorectal tumor (2-5). (B) The correlation of the 5hmC profile between tumor-derived, PDX plasma cfDNA and donor tumor gDNA depends on the number of passages of the PDX mouse. The size of the points is proportional to the size of grafted tumor, and the density of color denotes the growth rate of the grafted tumor. (C) Using the correlation distance of the top five genes that had the greatest 5hmC level in PDX plasma cfDNA, donor tumor gDNA and PDX plasma cfDNA from the same individual patient were clustered together. (D) Genes with elevated 5hmC level in cancer patient plasma cfDNA are enriched in genes with high 5hmC levels in PDX plasma cfDNA. To estimate fold enrichment, the fifth percentile genes in descending order of the log₂ fold change in patient cfDNA were compared against the corresponding fifth percentile genes in descending order of the 5hmC level in PDX mice (stomach pdx high, colon pdx high). Similarly, genes with decreased 5hmC level in cancer patient plasma cfDNA are enriched in genes with low 5hmC levels in PDX plasma cfDNA (stomach pdx low, colon pdx low). Dashed line denotes no enrichment.

Only the sequencing reads mapped to the human genome were further analyzed. Pearson's correlation of 5hmC profile between plasma cfDNA of PDX mice and gDNA of donor tumors significantly depends on the number of passages (P = 0.037, Figure 5B). This suggests a quantitative relationship between tumor growth and the experimental capture of tumor 5hmC in plasma cfDNA, as the size (P = 0.0096) and growth rate (P = 0.0080) of tumors grafted in PDX mice increase with passage numbers (Figure 5B). Using the top five genes with the greatest 5hmC levels in PDX plasma cfDNA, donor tumor and the derived PDX from the same individual patient were clustered together (Figure 5C), indicating donor tumor tissue as the human origin of the PDX cfDNA.

PDX allowed us to study tumor-derived cfDNA without confounding background cfDNA. Genes with greater 5hmC levels in tumor-sourced PDX plasma cfDNA were more likely to be the genes with elevated 5hmC level in plasma cfDNA in cancer patients. Indeed, we found that genes with increased 5hmC levels in patient plasma cfDNA were enriched in those genes with greater 5hmC levels in PDX plasma cfDNA, whereas genes with decreased 5hmC levels in patient plasma cfDNA were enriched in genes with lower 5hmC levels in PDX plasma cfDNA (Fisher's exact tests P < 1×10⁻⁹, Figure 5D). In contrast, genes with 5hmC level changed between tumor and adjacent tissues showed no such enrichment pattern (Figure 5D).

Preparation of cfDNA samples

cfDNA samples were prepared from peripheral blood collected from patients and healthy controls. Briefly, 4 ml of peripheral blood was collected from each subject using EDTA anticoagulant tubes, and the plasma sample was prepared within 6 h by centrifuging twice at 1 350× g for 12 min, and then centrifuging at 13 500× g for 12 min. The prepared plasma samples (about 2 ml/subject) were immediately stored at −80 °C. The plasma cfDNA was isolated using the QIAamp Circulating Nucleic Acid Kit (Qiagen) according to the manufacturer's protocol. Within each experimental batch, samples were randomized with respect to disease status in the following library preparation and sequencing profiling.

Isolation of gDNA from tissues

Tissue samples, including tumor and adjacent tissue samples, from patients were stored at −80 °C after surgical removal. About 10-25 mg tissue was collected using a scalpel after sample thawing. gDNA from tissues was isolated using the ZR Genomic DNA Tissue Kits (Zymo Research) according to the manufacturer's protocol.

5hmC-Seal-seq library preparation and sequencing

Seal-seq libraries for 5hmC profiling were prepared following our previously patented technology²⁹

In this method, the T4 bacteriophage β-glucosyltransferase is used to transfer an engineered glucose moiety containing an azide group onto the hydroxyl group of 5hmC across the human genome. The azide group is then chemically modified with biotin for affinity enrichment of 5hmC-containing DNA fragments. First, the gDNA is fragmented using an enzymatic reaction. Next, the fragmented gDNA or the cfDNA were repaired and installed with the Illumina compatible adaptors. The glucosylation reactions were performed in a 25 μl reaction containing 50 mM HEPES buffer (pH 8.0), 25 mM MgCl₂, purified DNA, 100 μM N₃-UDP-Glc and 1 μM βGT at 37 °C for 1 h. The reaction mix was subject to purification using a Micro Bio-Spin 30 Column (Bio-Rad) into ddH₂O. Subsequently, 1 μl DBCO-PEG4-DBCO (Click Chemistry Tools, 4.5 mM stock in DMSO) was added to the reaction mixture. The reactions were incubated at 37 °C for 2 h. Next, the DNA was purified ona Micro Bio-Spin 30 Column (Bio-Rad). The purified DNA was incubated with 5 μl C1 Streptavidin beads (Life Technologies) in 2× buffer (1× buffer: 5 mM Tris pH 7.5, 0.5 mM EDTA, 1 M NaCl) for 15 min according to the manufacturer's instruction. The beads were subsequently washed eight times for 5 min with 1× buffer. All binding and washing was done at room temperature with gentle rotation. The captured DNA fragments were amplified with 14-16 cycles of PCR amplification. The PCR products were purified using AMPure XP beads according to the manufacturer's instructions. The DNA concentration of each library was measured with a Qubit fluorometer (Life Technologies) and sequencing was performed on the Illumina Hi-Seq or NextSeq 500 platform.

RNA-seq library preparation and sequencing

Tumor and adjacent tissue samples including two colon samples and one stomach sample were collected and RNA was isolated using the ZR-Duet DNA/RNA Miniprep Kit (Zymo Research). Total isolated RNA was used to construct the library with the NEBNext Ultra RNA Library Prep Kit for Illumina following the manufacture's protocol. Sequencing reactions were executed on the NextSeq 500 platform using paired-end mode, yielding at least 32 M reads per sample.

Establishment of patient-derived tumor xenografts

The animal protocol for this study was reviewed and approved by the Ethical Committee of Medical Research, Shanghai Huashan Hospital of Fudan University. It utilized 6-8 week old BALB/c nu/nu mice weighing 16-20 g on receipt (Shanghai Laboratory Animal Center). The fresh pathological tissue fragments were placed in sterile tissue culture medium on ice and brought immediately to the animal facility. Tumor graft samples were cut into multiple 1 × 1 × 1 mm fragments in complete media. Tumor material was implanted into female BALB/c nu/nu mice under isoflurane anesthesia; all possible efforts were made to minimize suffering. A skin incision (0.3 cm) was subsequently made on the right mid-back. One tumor piece (1-3 mm) was inserted into each pocket and the skin was closed. Mice were regularly checked. When thetumor diameter reached 1.5 cm, mice were euthanized and tumors were excised, cut into 1 × 1 × 1 mm fragments again, and passaged to successive generations of three mice. The remaining tumor was snap frozen in liquid nitrogen and stored at −80 °C; plasma was also separated from blood sampled via the mouse eyeball. In this study, the gastric cancer and colorectal cancer patient-derived tumor xenografts (PDX) were randomly selected in our existing PDX model library, while the control group was BALB/c nu/nu mice of 12-14 weeks old.

5hmC enrichment analysis

We designed two similar spike-in probes with unique sequences, named 5hmC spike-in and no5hmC spike-in.

5hmC spike-in:

5′-CTGTCATGGTGACAAAGGCATCC*GGCAGAAATGCCCACACAGCCTCTTTAACCAGCACGCC
AACCGCCTCTGCTTCGGCCCTGGTCACGCAGCTGACAAGGTCTTCATAATAGAGAAATCCTG-3′, C* — 5hmC modifications.

no5hmC spike-in:

5′-CTGTCATGGTGACAAAGGCATCGCAGCGAAATGCCCACACAGCCTCTTTAACCAGCACGCC
AACCGCCTCTGCTTCGGCCCTGGTCACGCAGCTGACAAGGTCTTCATAATAGAGAAATCCTG-3′

These sequences cannot map to the human reference genome. Six cfDNA sequencing libraries were constructed from the same cfDNA (10 ng) sample, which were divided into control and experimental groups, each having three duplicates. About 100 million copies of 5hmC and no5hmC spike-ins were then mixed with the experiment sample before library preparation. The control group did not include the 5hmC pull-down step, whereas the experimental group included the 5hmC pull-down procedure. After sequencing, we extracted spike-in reads and calculated the enrichment ratios. The average ratio of 5hmC spike-in to no5hmC spike-in in the control group was 0.72, while the ratio in the experimental group was 40.36.

Technical stability analysis for 5hmC Seal-seq library preparation

The designed spike-in probes were utilized to improve the robustness and sensitivity of 5hmC-Seal. In total 20 000 copies of 5hmC and no5hmC spike-ins were pre-mixed and then added into the same cfDNA samples before constructing libraries. Different spike-in samples were designed as follows: 20 ng cfDNA with 2 repeats, 10 ng cfDNA with 10 repeats, 5 ng cfDNA with 2 repeats, 2 ng cfDNA with 2 repeats, 1 ng cfDNA with 2 repeats.

Total 5hmC quantification in cfDNA and gDNA

The enzymatic digestion protocol for each gDNA and cfDNA sample was the same. gDNA or cfDNA (all in 8 μl H₂O) was first denatured by heating at 95 °C for 5 min and then transferred into ice water for 2 min. Next, 1 L of 10× S1 nuclease buffer (30 nM CH₃COONa, pH 4.6, 260 mM NaCl, 1 mM ZnSO₄) and 180 units (1 μl) of S1 nuclease were added into the DNA solution. The mixture (10 μl) was then incubated at 37 °C for 4 h. Then 34.5 μl of H₂O, 5 μl of 10× alkaline phosphatase buffer (50 mM Tris-HCl, 10 mM MgCl₂, pH 9.0), 0.5 μl of alkaline phosphatase were added into the DNA digestion solution. The incubation was continued at 37 °C for an additional 4 h.

The CE-ESI-MS experiments were carried out with the CESI-8000 capillary electrophoresis (CE) system from Beckman Coulter coupled with a Sciex Tripel Quad 5500 Mass Spectrometer (Sciex) through a modified Nanosprayed II interface. Bare fused silica capillaries etched with porous tips were made available by Beckman Coulter; these could be inserted into the sheathless nanospray interface. The separation capillary was 100 cm long with an internal diameter of 30 μm and an outside diameter of 150 μm. The capillary was flushed with methanol for 10 min at 100 psi, followed by water, 0.1 M sodium hydroxide, 0.1 M hydrochloric acid and water for 10 min each at 100 psi, and finally by the background electrolyte (BGE) of 10 % acetic acid (pH 2.2) for 10 min at 100 psi before first use. The BGE was also used as conductive liquid in the conductive liquid capillary. Before each run, the conductive liquid capillary was rinsed with BGE for 5 min at 100 psi. Samples for detection were stored at 5 °C in the CE system. Hydrodynamic injections were used in this study, and about 100 nl sample was injected into the separation system for each analysis. A voltage of +25 kV was applied during the separation and the current was between 3.0 and 3.2 μA. The electrospray voltage was optimized to get the best nanospray stability and efficiency, and +1.7 kV was good enough for this study. The quantification calibration curves of 5′-dC, 5′-mdC and 5′-hmdC were constructed using mixture solution of their standards in different concentration. The resulting solutions of our DNA samples were directly measured by CE-ESI-MS. The concentrations of these three nucleosides in each sample were calculated based on the calibration curves. And the 5′-mdC/dC and 5′-hmdC/(dC + 5′-mdC) results of each sample were then calculated.

Sequencing data processing and detection of differential loci

Read-through sequences within raw sequencing reads were trimmed using Trimmomatic version 0.35³⁵. Low-quality bases at the 5′ (Phred quality score < 5) and 3′ (5 bp-sliding window Phred score < 15) were also trimmed. Reads with a minimum length of 50 bp were aligned to the human genome assembly GRCh37 using Bowtie2 version 2.2.6³⁶ with the end-to-end alignment mode. For paired-end sequencing data, read pairs were concordantly aligned with fragment length < 500 bp and with up to 1 ambiguous base, and four mismatched bases per 100 bp length. Alignments with Mapping Quality Score (MAPQ) ≥ 10 were counted for overlap with genomic features using featureCounts of subread version 1.5.0-p1³⁷, without strand information. Autosomal feature counts with > 10 mean counts across samples were then normalized and compared between-group using DESeq2 version 1.12.3³⁸. Since gender is not a significant covariate for both autosomal gene expression³⁹ and DNA methylation⁴⁰, while aging has been linked to DNA methylation⁴¹, age at sample collection/surgery was included as a categorical variable (< 20, 20-55, > 55 yr) in the negative binomial generalized linear model implemented in DESeq2. As for the experimental batches (discovery and validation sets), samples were processed in less than 1 week by 1-3 technicians, the identity of the technician was included in the model to adjust for potential technical correlation. When comparing tumor and adjacent tissues, patient identity was nested under technician identity. A FDR⁴² of 5% was used to identify differential 5hmC loci.

For PDX mouse plasma cfDNA data, sequencing reads were trimmed and aligned to a composite assembly of mixed human and mouse (GRCm38) genome. Unique alignments (MAPQ ≥ 10) were separated into human and mouse reads by chromosome name.

For RNA-seq data, sequencing reads were trimmed and aligned to GRCh37 annotated with GENCODE release 19, using STAR version 2.5.1b⁴³. Unique alignments with ≥ 90% match over reads were summarized by featureCounts. For the correlation analysis presented in Supplementary information, Figure S6A, 5hmC data were also summarized over exon regions as in RNA data. For genes having > 10 mean counts across samples, log₂ fold change between tumor and adjacent tissues were estimated by DESeq2 adjusting for patient identity.

Refining 5hmC biomarkers and evaluating performance

Cancer prediction models were trained using the differential 5hmC loci (e.g., gene bodies) detected in the discovery batch. We applied elastic net regularization on a logistic linear regression model⁴⁴, using the glmnet library in the R Statistical Package⁴⁴:

Where x is a J (1 ... j) by I (1 ... i) matrix of 5hmC level at gene j for sample i. The model is solved by

where P_α is the blend of the ridge (α = 0) and the lasso (α = 1) penalty. The parameter λ controls for the overall strength of penalty, while the parameter α controls for the relative proportion between the ridge and lasso penalty. For a given α, λ is estimated by cross validation and selected as the largest λ at which the mean cross-validated error is within one s.e. of the minimum. As we expect that 5hmC loci with a larger effect size are more robustly detected by the assay and therefore more reproducible across experiments, we included in the model a penalty factor ν_j so that

where f_j is the log₂ fold change of gene j estimated in the training batch.

The parameter α reflects the model assumption (i.e., a large number of small effects or a small number of large effects). In small data sets like our discovery batch, the selection of α based on residual errors may lead to an over-fitted model. Instead, our model assumption was guided by the validation batch, so that α was searched to maximize AUC in the validation batch over a grid of values from 0.05 to 0.95. The model, derived from the training batch using the selected α, was applied in all classifications as described in Figures 3 and ​and44.

The normalization of data in cancer classification adopts the regularized log transformation implemented in DESeq2, which estimates a global mean dispersion trend to shrink the variance at low count genes that are associated with high Poisson noise, so that variance is stabilized across genes in the log transformed data. In Figure 3 and Supplementary information, Figure S6, data from the validation batches were all normalized to the reference distribution derived from the training batch and used directly in cancer classification, i.e., we essentially ignored any remaining batch effect, which could result from library preparation and the sequencing run. This is because under a real clinical setting, the batch effect estimated between testing samples and training samples will generally be biased due to the highly unbalanced case/control proportion in testing samples (low incidence of cancers). Batch effect may introduce some deviation from the 0.5 probability cutoff in cancer calling. External spike-ins may be used to estimate batch effect in future investigations.

ROC curves⁴⁵ were generated to evaluate the performance of a prediction algorithm, using the pROC⁴⁶ library in the R package. Sensitivity and specificity were estimated at the score cutoff that maximizes the sum of sensitivity and specificity using the ROCR⁴⁵ library in the R package.

Statistical analyses

For Figure 5A, the P-value was estimated by a linear mixed effects model: proportion of human reads (square root transformed) ~ xenograft status (none|colorectal|gastric) + γ (tumor donor identity) + ε, random effect γ was introduced to control for correlation among replicate xenografts. For Figure 5B, P-value was estimated by a linear mixed effects model: Pearson's r with tumor donor ~ number of passages + γ (tumor donor identity) + ε.

Annotating 5hmC loci with genomic features and functional analysis

The genomic features analyzed included promoters (3 kb upstream of gene start), gene bodies (gene start to stop sites annotated by GENCODE release 24)⁴⁷, CpG islands (annotated CpG islands by UCSC Table Browser plus +/ 1 kb region), and the ENCODE³¹ features (DNase I hypersensitive sites, TF binding sites and histone modifications)³¹. Each type of ENCODE feature from the profiled ENCODE cell lines was integrated into a single list of features by collapsing overlapped and nearby (< 150 bp) peaks. The ENCODE features analyzed in the 5hmC genomic distribution (Supplementary information, Figures S2 and S3) were as originally annotated without collapsing, with 20 000 features randomly sampled in the interquartile size distribution for each feature category. Pathway enrichment analysis of genes with cancer patient-associated 5hmC loci were explored based on the Kyoto Encyclopedia of Genes and Genomes (KEGG)⁴⁸ using the NIH/DAVID tool⁴⁹.

This work was supported in part, by grants from the National Institutes of Health (U01CA217078 to MBB and WZ, R33CA225338 to CH and WZ, P30CA060553 Career Development Fund to WZ, R01HG006827 to CH), The University of Chicago Ludwig Center for Metastasis Research (to CH), and the Research Special Fund for Public Welfare Industry of Health (201402001 to YZ). We thank Dr H Pickersgill (Life Science Editors) for editing this manuscript.