pmc Nucleic Acids Res Nucleic Acids Res 4 nar 0411011 nar Nucleic Acids Research 0305-1048 1362-4962 Oxford University Press PMC12135180 PMC12135180.1 12135180 12135180 40464685 10.1093/nar/gkaf482 gkaf482 1 AcademicSubjects/SCI00010 Narese/39 Methods Methyl-CODEC enables simultaneous methylation and duplex sequencing https://orcid.org/0000-0001-9059-6335 Liu Ruolin Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States https://orcid.org/0000-0002-4957-9004 Darbeheshti Farzaneh Dana–Farber Cancer Institute , Boston , MA 02142 , United States Walsh Laurel Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States Li Rachel Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States https://orcid.org/0000-0002-4059-5198 Bae Jin H Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States Zeggar Hayet Radia Dana–Farber Cancer Institute , Boston , MA 02142 , United States Narayan Azeet Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States https://orcid.org/0000-0002-0639-5803 Xiong Kan Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States https://orcid.org/0000-0001-7598-2406 Makrigiorgos G Mike Dana–Farber Cancer Institute , Boston , MA 02142 , United States Adalsteinsson Viktor A Broad Institute of MIT and Harvard , Cambridge , MA 02215 , United States To whom correspondence should be addressed. Email: mike_makrigiorgos@dfci.harvard.edu Correspondence may also be addressed to Viktor A. Adalsteinsson. Email: viktor@broadinstitute.org Ruolin Liu and Farzaneh Darbeheshti should be regarded as Joint First Authors. 10 6 2025 04 6 2025 53 10 489042 gkaf482 20 5 2025 04 4 2025 01 2 2025 04 06 2025 04 06 2025 11 06 2025 © The Author(s) 2025. Published by Oxford University Press on behalf of Nucleic Acids Research. 2025 https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License ( https://creativecommons.org/licenses/by-nc/4.0/ ), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com . Abstract DNA mutations and methylation often contribute to disease development in a synergistic manner. While duplex sequencing is the most accurate method for detecting DNA mutations, it typically lacks the ability to simultaneously assess methylation or requires many reads. Here, we developed Methyl-CODEC to enable simultaneous methylation sequencing and duplex sequencing using single read pairs. To achieve this, Methyl-CODEC links an enzymatically deaminated sense strand to the reverse complement of the antisense strand, which is protected from conversion by using conversion-resistant dCTPs in the strand linking step. Methyl-CODEC shows high concordance with standard enzymatic or bisulfite-based whole genome methylation sequencing, while also uniquely preserving the original DNA sequence. We show that hydroxy-methyl-dCTP is superior in this regard relative to other conversion-resistant dCTPs. Methyl-CODEC improves genetic sequencing accuracy, enables better read alignment for next-generation sequencing, and distinguishes C > T mutations from unmethylated Cs. It also identifies rare mutations including those producing methylated Cs, which are enriched in CpG contexts. Methyl-CODEC opens new horizons for enhanced detection of biomarkers in cancer and molecular medicine. Graphical Abstract Graphical Abstract National Cancer Institute 10.13039/100000054 2R01CA221874-04 pmc-status-qastatus 0 pmc-status-live yes pmc-status-embargo no pmc-status-released yes pmc-prop-open-access yes pmc-prop-olf no pmc-prop-manuscript no pmc-prop-legally-suppressed no pmc-prop-has-pdf yes pmc-prop-has-supplement yes pmc-prop-pdf-only no pmc-prop-suppress-copyright no pmc-prop-is-real-version no pmc-prop-is-scanned-article no pmc-prop-preprint no pmc-prop-in-epmc yes pmc-license-ref CC BY-NC Introduction DNA molecules carry not only genetic information within their canonical nucleobases but also crucial epigenetic modifications. Together, these genetic and epigenetic alterations influence cellular behavior in a synergistic manner, and a comprehensive assessment of both is essential for effective disease detection and treatment. While duplex sequencing [ 1 ] is considered the gold standard for highly accurate genetic sequencing, high numbers of reads are often required when both strands are sequenced separately. For epigenetic analyses, methods such as whole genome bisulfite sequencing (WGBS) [ 2 ] and enzymatic methyl sequencing (EM-seq) [ 3 ] are widely utilized. Both methods convert unmethylated cytosine (C), which is the predominant state of cytosine in the genome (>95%), into uracil (U) and subsequently into thymine (T). This process results in requiring a “three-letter” genome for read alignment, which can hinder alignment accuracy and efficiency. Paired-end sequencing with longer read lengths provides an advantage in mappability [ 4 ] and can partially mitigate the alignment challenges associated with converted DNAs [ 5 ], such as those in WGBS or EM-seq. Furthermore, TET-assisted pyridine borane sequencing [ 6 ] and direct methylation sequencing (DM-seq) [ 7 ] have been introduced, offering improved mapping due to selectively converting 5-methylcytosine (5mC) instead of unmethylated C. However, a significant limitation of all these methods is their difficulty in distinguishing C to T conversion via enzymatic or chemical deamination from C > T mutation, which is the most prevalent mutation in nature [ 8 , 9 ]. Other technologies capable of simultaneous epigenetic and genetic sequencing have been developed. Third-generation sequencing technologies, such as Oxford Nanopore Technologies [ 10 , 11 ] and Pacific Biosciences [ 12 ], enable direct readouts of both epigenetic and genetic sequences. However, their genetic sequencing accuracy is typically lower than that of short-read sequencing and falls far below the accuracy of duplex sequencing, except when consensus reads are utilized, but this compromises throughput. Additionally, they often require high input amounts of unamplified DNA, which constrains their applicability to samples with limited DNA, such as liquid biopsies. For short-read sequencing, methods for simultaneous epigenetic and genetic sequencing have been developed that make a copy of single-stranded DNA and use the original strand for methylation detection and the copied strand for mutation detection. Methods like five-letter sequencing [ 13 ] and Methyl-SNP-seq [ 14 ] allow simultaneous sequencing of epigenetic and genetic data in a single workflow, but their detection of rare mutations is limited by single-strand accuracy (error rate ∼10 −4 ). MethylSaferSeqS [ 15 ] copies both strands of a DNA duplex using a biotinylated primer to permit the separation of the original and copied strands. This method produces two separate libraries and enables both duplex sequencing and methylation sequencing but requires sequencing each strand on its own, which limits its application for whole-genome duplex sequencing. For assessing regional aggregated methylation levels alongside genome sequencing in a single workflow, RIMS-seq2 [ 16 ] employs limited deamination to achieve 1% C to T transitions at 5mC sites. Yet, it does not provide single-base methylation resolution and has significant limitations in detecting rare somatic C > T mutations. At present, methods for simultaneous epigenetic and genetic sequencing have limitations in accuracy, efficiency, or breadth. Here, we introduce Methyl-CODEC for simultaneous methylation and duplex sequencing using single read pairs. The process follows the strategy we originally described for CODEC ( C oncatenating O riginal D uplex for E rror C orrection), where the sense strand is linked with the antisense strand complement. CODEC produces DNA molecules where the sequence of both strands can be compiled in single sequence reads, thus enabling genome-wide, duplex sequencing accuracy [ 17 ]. Methyl-CODEC, on the other hand, enables the linkage of an enzymatically methyl-converted sense strand with the copied antisense strand of the native sequence, and vice versa. This retains the duplex sequencing accuracy of CODEC while providing uncompromised alignment and single-base methylation information, allowing for simultaneous, genome-wide detection of methylation and mutations in single DNA duplexes using single read pairs. Materials and methods Samples and oligonucleotides Universal methylated and non-methylated standard DNA of HCT 116 DKO samples were purchased from ZYMO Research (catalog no. D5014). Genomic DNA of the NA12878 cell line was obtained from Coriell Institute of Medical Research (catalog no. GM12878). DNA fragmentation was performed by a Covaris ultrasonicator to achieve a mean fragment size of 150 bp. The fragmented DNA was purified via AMPure (Beckman Coulter), followed by elution in low TE buffer. Twenty nanograms of fragmented DNA was used as input for library preparation. To prepare a 95% non-methylated library, we mixed fragmented non-methylated standard DNA with fragmented methylated standard DNA in a 19:1 ratio. Methyl-CODEC adapter oligonucleotides containing methylated cytosine, instead of unmodified cytosine, were ordered from Integrated DNA Technologies with polyacrylamide gel electrophoresis purification. Three different sets of adapter oligonucleotides with distinct indices were purchased ( Supplementary Table S1 ). The Methyl-CODEC adapter quadruplex was formed by adding 5 μM of each oligonucleotide to low TE buffer with 100 mM NaCl in a total volume of 40 μl. The mixture was then heated at 85°C for 3 min, gradually cooled to 20°C, and incubated at room temperature for 12 h. Methyl-CODEC library preparation End-repair/dA-tailing and ligation NEBNext Ultra II DNA Library Prep Kit for Illumina (catalog no. E7645) was used according to the manufacturer’s protocol ( Supplementary methods ). Phi29 strand displacing extension For the strand displacement step, Phi29 DNA polymerase (NEB, catalog no. M0269L) was applied by adding the following to 41.5 μl of ligated DNA from the previous step, followed by incubation at 30°C for 20 min: 10 μl of 10× buffer, 1 μl of polymerase, 47 μl of nuclease-free water, and 2 μl of 10 mM customized dNTPs, including dATP, dTTP, dGTP, and conversion-resistant dCTP analogues, evenly. We used four different conversion-resistant dCTPs, including 5-methyl-dCTP (NEB, catalog no. N0356S), 5-carboxy-dCTP (TriLink, catalog no. N-2063), 5-hydroxymethyl-dCTP (Jena Bioscience, catalog no. NU-932S), and 5-propargylamino-dCTP (Jena Bioscience, catalog no. NU-809-1). The Phi29 extended product was purified using 0.75× AMPure XP beads, followed by final elution in 28 μl of nuclease-free water to prepare for the enzymatic methyl-conversion step. Enzymatic methyl conversion NEBNext Methyl-seq Conversion Module (NEB, catalog no. E7125S) was used for DNA conversion based on the manufacturer’s protocol ( Supplementary methods ). Library amplification Methyl-CODEC library amplification was performed by adding 25 μl Q5U Hot Start High-Fidelity DNA Polymerase (NEB, catalog no. M0515) and 5 μl KAPA Library Amplification Primer Mix (Roche, catalog no. KK2623) to 20 μl Methyl-CODEC library molecules. Polymerase chain reaction (PCR) was run via Mastercycler X50 Eppendorf ( Supplementary methods ). The libraries were quantified and qualified by Qubit dsDNA High Sensitivity Kit (Invitrogen, catalog no. Q33230) and Agilent High Sensitivity DNA Kit (catalog no. 5067-4626). Sequencing The four Methyl-CODEC NA12878 samples were sequenced using NovaSeq S4 with no index cycle and the read length of 169, whereas the three Methyl-CODEC HCT116 DKO samples were sequenced using NovaSeq X with eight index cycles and the read length of 161. The four NA12878 samples had 930M (hydroxy-C), 166M (5mC), 999M (carboxy-C), and 419M (propargylamino-C) read pairs. The three HCT116 samples yielded 834M (methylated control), 820M (unmethylated control), and 772M read pairs (mixed). The sequencing metrics are available in Supplementary Table S2 . Methyl-CODEC analysis workflow Methyl-CODEC uses the same demultiplexing and adapter trimming steps as CODEC [ 17 ]. The rest of the workflow proceeds as follow (see also Supplementary Fig. S2 ): Confirming read structure : After adapter trimming, Methyl-CODEC includes an additional step to determine whether a read pair has the correct structure (i.e. one read is from the converted strand and the other is from the protected strand) and to identify its specific structure (i.e. product 1 versus product 2). This is achieved using a two-step algorithm. In the first step, if one read has a fraction of G <5% and the other has a fraction of G >5%, the former is marked as the converted strand/read and the latter as the protected strand/read. If both reads have a fraction of G either all above 5% or all below 5%, the second step is applied. In the second step, overlap alignment is performed, aligning the suffix of the first read to the prefix of the reverse complement of the second read using custom alignment parameters (e.g. no mismatch penalty for conversion induced C/T and G/A mismatches and high penalties for gaps). The C/T and G/A mismatches in the overlap region are counted. If the number of C/T mismatches is greater, the read pair is identified as product 1 (with the first read protected and the second converted). Conversely, if the number of G/A mismatches is larger, the read pair is identified as product 2 (with the first read converted and the second read protected). If no overlap alignment is found, the read pair is marked as a by-product. The first step accelerates the process and accommodates cases where the two reads have little or no overlap. Alignment of the protected reads : Methyl-CODEC utilizes standard alignment algorithms without requiring methylation specific aligners. After confirming the read structure (i.e. product 1 versus product 2), Methyl-CODEC uses standard BWA-MEM [ 18 ] to align the protected read to the entire reference genome (HG38 in this study). If the alignment is unsuccessful, the read pair is marked as a by-product. For successfully aligned reads, the reference sequences of the alignments are extracted, including 500 base pairs as padding. Alignment of converted reads : The converted reads are aligned to the local reference sequences obtained from the protected strand alignments using BWA-SW [ 16 ]. Custom alignment parameters are applied here as well, avoiding penalties for conversion-induced mismatches (i.e. C/T mismatches for product 1 and G/A mismatches for product 2). Reads with an alignment score below 30 are classified as by-products. An entire read pair will be classified as a by-product if either the protected strand or the converted strand has an alignment score below 30 ( Supplementary Fig. S2 ). Deduplication : For all Methyl-CODEC product reads, Picard MarkDuplicates is used to mark PCR duplicates using unique molecular identifiers (UMIs), start and stop positions of the aligned paired-end reads. It is important to note that the UMIs incorporated into the adapters include only non-C bases, while the copied portion incorporates deamination-resistant Cs, ensuring the UMIs remain intact and unaffected by deamination ( Supplementary Fig. S1A and Supplementary Table S1 ). Among all PCR-duplicated read pairs of a single duplex, only the one with the best sum of quality scores was used for downstream analysis in order to confirm that a single read pair is sufficient to achieve duplex accuracy. However, generating a consensus from duplicates may further improve accuracy. Determining methylation status and cytosine contexts : After obtaining deduplicated paired-end alignments, the methylation statuses of cytosines in the overlapped regions are determined. For product 1, the Cs in the protected strands are first scanned. A methylated C or a non-methylated C is identified when the corresponding base on the converted strand is C or T, respectively, providing that the two complementary bases from each strand meet the quality score threshold (Q30). The same process applies to product 2 but focuses on G in the protected strand and G or A in the converted strand. In both situations, the 12 base pairs from the fragment ends are ignored. These methyl-resolved cytosines are then annotated based on their contexts (CpG, CHG, or CHH, H is an IUPAC code that stands for A, C, or T) in a reference-independent manner. The context is called only if the adjacent bases from the two strands agree (accounting for conversion induced C/T, G/A mismatches) and either both meet the quality score threshold (Q30) or match the reference base. Cytosines in the CpG context are annotated with a Bismark-like symbol Z (methylated) or z (unmethylated). CHG and CHH contexts are determined based on the further adjacent bases. Finally, all methyl-CODEC products with methylation status and context information are stored as a BAM file. All steps (i)–(iv) are encapsulated in a single executable file. Alignments are performed in memory using direct calls to the BWA API in C++ code to minimize IO compute time. Methylation extraction After generating BAM alignments with Bismark-like tags, fragment-level and locus-level methylation information was extracted using bismark_methylation_extractor, which also calculated sample-level methylation fractions for CpG, CHG, and CHH contexts. Methylation fractions for CpG islands were calculated by intersecting the locus-level methylation with CpG island coordinates. Processing EM-seq WGMS and WGBS A 20-ng input EM-seq library and a 50-ng input WGBS library were downloaded from the EM-seq paper [ 3 ]. Reads were then aligned using Bismark (default mode), and methylation information was extracted with bismark_methylation_extractor, ignoring the first 12 bases from the read 2 (parameters: --ignore 2 --ignore_3prime 2 --ignore_r2 12 --ignore_3prime_r2 2) to adjust for M-bias [ 19 ]. Simulating paired-end EM-seq WGMS from Methyl-CODEC for comparison Paired-end EM-seq whole genome methylation sequencing (WGMS) data were simulated from Methyl-CODEC products by in silico conversion of the protected strand in a product specific manner (converting C to T for product 1 or G to A for product 2). The resulting simulated paired-end reads had both ends converted while maintaining the corrected pair orientation. They were then aligned using Bismark in two modes (default and local). The local mode allows soft-clipping at the ends, while default mode enforces end-to-end alignment. One NA12878 and three HCT116 samples were simulated. In NA12878, the Methyl-CODEC alignments were compared with Bismark alignments. In three simulated HCT116 samples, after alignment by Bismark in two modes, methylation results were extracted from the original converted strands, excluding the first and last 12 bp in the reads, and then compared with Methyl-CODEC results. Extracting HCT116 methylation levels in CpG sites with and without overlapping germline mutations Germline mutations in the HCT116 parental cell line were identified using GATK HaplotypeCaller [ 20 ]. Germline mutations (SNPs and indels) overlapping CpG sites in the HG38 genome were further filtered to include only C > T and G > A single nucleotide polymorphisms (SNPs). Methylation levels were calculated for CpG sites overlapping C > T and G > A SNPs and for CpGs that were unaffected by germline mutations, i.e. CpGs that were at least 3 bp away from any germline mutations. Somatic mutation detection Methyl-CODEC uses the CODECsuite to call somatic mutations from single read pairs, with modifications to account for methylation conversion: (i) Mutations are called if a base C from a protected strand is paired with a T from the converted strand (product 1), or G from the protected strand is paired with A from the converted strand (product 2). (ii) The conversion-induced mismatches (C/T and G/A) are excluded from all filters. (iii) Soft-clipping below a certain length is allowed at the 5′ end (default is 3 bp, adjustable via command-line options). (iv) Alignment score filter is turned off for the converted strands, though mapping quality filter is still applied. To calculate residual single nucleotide variant (SNV) frequencies by mutation contexts in NA12878, mutations detected in product 2 were reverse complemented (e.g. C > T became G > A) and combined with product 1 mutations. For the denominators (number of A,C,G,T bases) in the residual SNV frequencies, we complemented the bases detected in product 2 (e.g. C became G and G became C) and combined with the bases detected in product 1. Results Methyl-CODEC workflow Methyl-CODEC enables simultaneous epigenetic and genetic sequencing with the accuracy of duplex sequencing by linking a methyl-converted sense strand with the reverse complement of the antisense strand of the native sequence, and vice versa. This is achieved by using conversion-resistant dCTP analogs in place of standard dCTPs during the strand displacement step of CODEC [ 17 ], followed by EM-seq conversion [ 3 ] (Fig. 1A and Supplementary Fig. S1A ). Specifically, after end repair and ligation to the CODEC adapter quadruplex—in which standard dCTPs are replaced by conversion-resistant dCTPs—strand displacing extension occurs using a mix of conversion-resistant dCTP and regular dATP, dGTP, and dTTP. These conversion-resistant dCTPs are protected from being converted to Us (and subsequently to Ts during PCR) by APOBEC ( Supplementary Fig. S1B ). This process allows the copied side of each original strand to remain intact after methylation conversion (referred to as the protected strand, in contrast to the original/converted strand). From each input DNA duplex, Methyl-CODEC creates two products: one with the sense strand converted and the other with the antisense strand converted (Fig. 1A ). This retains complete epigenetic and genetic information from both strands in the library and allows detection of hemi-methylation when both Methyl-CODEC products are sequenced ( Supplementary Fig. S1E ). In practice, if the library is sequenced at low depth, we often uncover only one of the two Methyl-CODEC products from each DNA duplex. This enables simultaneous methylation sequencing of one strand and duplex sequencing of both strands of each DNA duplex using a single read pair. Figure 1. Methyl-CODEC library preparation and analysis workflows. ( A ) Methyl-CODEC captures both genetic and epigenetic information in a single workflow. A modified quadruplex adapter, which is resistant to conversion, is used for ligation. The strand-displacing extension is performed using Phi29 polymerase with conversion-resistant dCTPs, instead of regular dCTPs, followed by enzymatic methyl conversion. This results in two products, each comprising the converted sequence of one strand (dashed line) concatenated with the reverse complement of the other strand with conversion protection (solid line with green background). UMIs, which include only non-C bases, are incorporated into the adapters to facilitate the re-association of the two products as required, although a single Methyl-CODEC product is sufficient to achieve duplex accuracy. Arrows indicate the sequencing directions. ( B ) Methyl-CODEC analysis workflow. (i) The strand information (protected versus converted) is determined for a given read based on G% or C/T or G/A mismatches. (ii) For a given duplex read pair, the protected strand is first aligned to the standard four-letter reference genome (HG38) by BWA-MEM. The local reference sequence is extracted from the alignment of the protected strand. (iii) The converted strand is aligned to the extracted local reference sequence by BWA-SW using a custom setting that does not penalize C/T or A/G mismatches. (iv) PCR duplicates are marked using UMIs and start and stop positions of the aligned paired-end reads, and only the fragment with the best sum of quality scores was used for downstream analysis. (v) A,G,T, unmethylated C, and methylated C are determined using the duplex informed deconvolution table ( Supplementary Fig. S1C ) from comparing the sequences of the two strands. Mutations are called by comparing to the reference. Owing to the use of complementary strands, Methyl-CODEC achieves duplex accuracy for mutation detection. See the “Materials and methods” section and Supplementary Fig. S2 for further details. A correct Methyl-CODEC product has one read of a pair derived from the protected strand and the other read from the converted strand, originating from opposing strands of the same original DNA duplex (i.e. duplex reads) (Fig. 1A ), whereas Methyl-CODEC by-products are generated in similar ways as CODEC by-products. These by-products primarily result from mechanisms such as double ligation, blank ligation, or intramolecular ligation [ 17 ]. Moving forward, we will refer to the correct Methyl-CODEC product simply as “Methyl-CODEC product” to distinguish it from Methyl-CODEC by-products. To determine Methyl-CODEC product reads, the protected strands are first aligned to the entire reference genome using BWA-MEM [ 18 ]. Next, the converted strands are aligned strictly to regions where protected strands were aligned, with an additional 500-bp padding, using custom BWA-SW with zero mismatch penalty for conversion induced C/T or G/A mismatches (Fig. 1B and Supplementary Fig. S2 ). Only reads with a minimum alignment score of 30 for both strands are considered Methyl-CODEC products; other fragments are classified as by-products but retained for further analysis. From deduplicated Methyl-CODEC products, A, C, G, T, and 5mC calls are made based on the deconvolution table ( Supplementary Fig. S1C ) in overlapping regions of the two strands and subject to fragment level and base level filtering, similar to CODEC (see the “Materials and methods” section). Methylation calls are stored as Bismark-like tags [ 21 ] for compatibility of downstream analysis. We first confirmed that the Methyl-CODEC workflow can generate Methyl-CODEC products using the human cell line NA12878. Using four different conversion-resistant analogs (5m-dCTP, carboxy-dCTP, hydroxy-dCTP, and propargylamino-dCTP), 47.6%, 76.6%, 71.2%, and 76.8% of the total reads, respectively, were confirmed as Methyl-CODEC products, with the majority (∼95%) of Methyl-CODEC by-products being intermolecular by-products ( Supplementary Fig. S1D and Supplementary Table S2 ). It remains unclear whether the observed lower yield for 5m-dCTP is due to the intrinsic variability of the experiments or a potential inefficiency of Phi29 polymerase in incorporating 5m-dCTP into the extended strands. Overall, hydroxy-dCTP was selected as the conversion-resistant dCTP variant for Methyl-CODEC due to its most complete protection against enzymatic methyl-conversion, exhibiting the lowest residual C > T frequency among the four variants ( Supplementary Fig. S3 ). Additionally, hydroxy-dCTP resulted in the lowest fraction of unmethylated cytosines (umC) dropouts caused by false C > T mutations ( Supplementary Fig. S4A ) and the highest Pearson correlation (0.84) of CpG methylation with EM-seq among the four analogs ( Supplementary Fig. S4B ). Methyl-CODEC matches EM-seq WGMS and WGBS for methylation sequencing Methyl-CODEC utilizes the same enzymatic methyl-conversion step as EM-seq [ 3 ], while WGBS remains a widely used approach for methylation analysis. Therefore, we compared methylation results from whole genome Methyl-CODEC sequencing with WGMS data from EM-seq and WGBS on the same cell line, NA12878, obtained from public repositories [ 3 ]. For alignment of EM-seq and WGBS data, Bismark [ 21 ] was used. First, we found that genome-wide average cytosine methylation levels across three different contexts were similar (Fig. 2A ). Additionally, the per-CpG island methylation level from Methyl-CODEC showed high concordance with EM-seq and WGBS, achieving Pearson correlations of 0.99 and 0.97, and strong linear fits with adjusted R 2 of 0.97 and 0.95, respectively (Fig. 2B ). Across individual CpG sites genome-wide, Methyl-CODEC also showed a similar distribution (Fig. 2C ) and achieved high Pearson correlations with EM-seq WGMS (0.84) and WGBS (0.82) (Fig. 2D ). Figure 2. Concordance of 5mC readout between Methyl-CODEC and EM-seq or WGBS in NA12878 cell line. ( A ) Genome-wide methylation level at three different sequence contexts (CpG, CHG, and CHH). H is an IUPAC code that stands for A, C, or T. ( B ) Scatter plots comparing methylation level per CpG island ( n = 25 368) between Methyl-CODEC and EM-seq, and between Methyl-CODEC and WGBS. The blue lines show linear fits made using R ggplot2. Pearson correlations and adjusted R -squares are shown in the top left corners. ( C ) Violin plots comparing the genome-wide distribution of CpG methylation levels. ( D ) Pairwise Pearson correlations (top right corner) and genome-wide scatter plots of 9 532 983 CpG methylation, with a mean depth of 3 in all samples. This plot was generated by methylKit [ 44 ]. EM-seq WGMS is known to result in less GC bias compared to WGBS [ 3 ]. We found that Methyl-CODEC closely resembled the GC profile of EM-seq WGMS, maintaining similar coverage across various GC contents ( Supplementary Fig. S5A ). Interestingly, all three technologies exhibited a similar relationship between fragment-level methylation and GC content: methylation level steadily increased, and then decreased as GC content rose, peaking around 50% GC content ( Supplementary Fig. S5B ). Unlike EM-seq and WGBS, Methyl-CODEC uniquely preserves the original DNA sequence for next-generation sequencing (NGS) alignment. To assess the impact of this on read alignment, we first removed duplicates and then simulated paired-end WGMS from NA12878 Methyl-CODEC and aligned them by Bismark (see the “Materials and methods” section). We compared the concordance between the Methyl-CODEC alignments based on the original four-letter DNA sequence, which is how Methyl-CODEC is configured, and the Bismark alignment based on the methyl-converted sequence (three-letter) as per WGBS or EM-seq. After applying mapping quality thresholds for both Methyl-CODEC and Bismark, we observed a sharp difference when using different Bismark alignment parameters ( Supplementary Fig. S6A ). In particular, Bismark’s local alignment mode significantly reduced the number of unaligned reads compared to the default mode (11.3% versus 28.9% unaligned). However, this improvement came with a slight increase in the fraction of misplaced/discordant reads: 0.38% in local mode compared to 0.06% in default mode ( Supplementary Fig. S6A ). Because Methyl-CODEC does not rely upon a “three-letter” genome for read alignment, our results suggest that Methyl-CODEC may enable more accurate read alignment for methylation sequencing, although its impact was limited to a minor subset of reads with likely incorrect alignments. In addition, Methyl-CODEC can distinguish C > T mutations, including germline and somatic or rare mutations, from C to T conversions via enzymatic or chemical deamination, given its retention of the original DNA sequence. To explore this, we prepared and sequenced three additional Methyl-CODEC libraries using human HCT116 DKO samples: unmethylated control DNA, methylated control DNA, and a 1:19 mix of these two. Unmethylated human DNA is purified from the HCT116 DKO cell line, which carries genetic knockouts of the DNA methyltransferases DNMT1 (−/−) and DNMT3b (−/−). The methylated DNA control is derived from HCT116 DKO DNA that has been enzymatically methylated specifically at CpG sites. From the parent HCT116 cell line [ 22 ], we detected 590 281 germline mutations (thus inherited by all HCT116 DKO cell lines), including indels, at CpG sites, with about 80% of these germline mutations being C > T/G > A homozygous and heterozygous SNPs ( Supplementary Fig. S6B ). We hypothesized that for EM-seq and WGBS, this could potentially result in false calling of hundreds of thousands of unmethylated CpG sites, which could lead to an appreciable underestimation of global methylation. Indeed, using the methylated control genome, we observed 97.3% global CpG methylation with Methyl-CODEC, compared to 95.8% (Bismark-default) and 92.4% (Bismark-local) global CpG methylation with simulated paired-end EM-seq WGMS ( Supplementary Fig. S6C ). We then calculated CpG methylation levels both outside and within the germline C > T/G > A sites in the methylated control. Methyl-CODEC consistently maintained global methylation levels in both cases (with 97.3% outside the germline C > T/G > A sites and 97.4% within these sites, closely aligning with the global average of 97.3%). In contrast, simulated EM-seq WGMS aligned by Bismark-default achieved a similar methylation level (97.0%) to Methyl-CODEC outside the germline C > T/G > A regions ( Supplementary Fig. S6D ) but dropped significantly to 31.6% at germline C > T/G > A sites ( Supplementary Fig. S6E ). Of note, this shortfall of EM-seq WGMS could be overcome by also performing whole genome sequencing on the same sample, but this introduces added complexity. For the unmethylated control genome, which still contains low levels of DNA methylation [ 23 ], all methods showed similar global CpG methylation levels at 13.4%, 13.4%, and 13.1%, respectively. While Methyl-CODEC maintained 13% methylation at CpG sites that overlap germline C > T and G > A mutations, Bismark-default and Bismark-local dropped to 4% methylation at those sites, consistent with the observation that Bismark alignment of standard EM-seq cannot distinguish between the conversion of umC versus real C > T mutations. In all, Methyl-CODEC is superior to EM-seq and WGBS in distinguishing C > T mutations from C to T conversions via enzymatic or chemical deamination. Methyl-CODEC enables accuracy equivalent to duplex sequencing and CODEC Duplex sequencing [ 1 ], as well as CODEC, rely on the principle that sequence information obtained from both DNA strands can distinguish genuine mutations from sequencing noise [ 24 ]. To determine whether duplex sequencing accuracy is retained for Methyl-CODEC, we analyzed the NA12878 cell line for mutation detection (see the “Materials and methods” section). This analysis yielded a raw deduplicated depth of 25× (duplex depth of 5×), with 41% of duplexes covered by just single read pairs of Methyl-CODEC products. Our results confirmed that Methyl-CODEC products have the required accuracy across all base contexts, including C > T/G > A mutations mapped to the “G” reference bases (Fig. 3A ). This is achieved because the G nucleotide on the strand complementary to a C (referred to as the “G-strand” in Bis-SNP [ 25 ]) remains unaffected by enzymatic methyl-conversion. On the other hand, the C > T/G > A mutations mapped to “C” reference are prone to false detection, possibly due to (i) PCR or sequencing errors and/or (ii) incomplete protection (Fig. 3B ). Figure 3. Methyl-CODEC achieves duplex sequencing level accuracy (>Q60). ( A ) A scatter plot comparing residual SNV frequencies of Methyl-CODEC (raw deduplicated depth at 25×, duplex depth at 5×) and standard CODEC in NA12878, both using a commercial end-repair/dA-tailing protocol. Methyl-CODEC data is split into pyrimidines (C or T) as the reference and purines (G or A) as the reference, while CODEC data represent the average across both. Error bars indicate 95% binomial confidence intervals by Wilson method. Mutations that produce methylated cytosines (5mC) and umC are labeled next to the corresponding dots. ( B ) A hypothetical double stranded DNA example demonstrating how Methyl-CODEC strand information (product 1 versus product 2) can be used to achieve duplex accuracy for C > T/G > A contexts and how T > umC can result from sequencing errors. “rc” indicates reverse complement. Additionally, we hypothesized that impurities in commercially available dCTP analogs can lead to incomplete protection from enzymatic conversion. Indeed, we found that hydroxy-dCTP gave the best C > T accuracy (a 2-log reduction in error rate compared to 5m-dCTP) when “C” is the reference ( Supplementary Fig. S3 ). The mutation frequency of C > T was close to the raw Illumina sequencing error rate (1 × 10 −3 ), suggesting that incomplete protection was almost eliminated by using hydroxy-dCTP. On the other hand, propargylamino-dCTP yielded worse accuracy, and carboxy-dCTP yielded accuracy between that of 5m-dCTP and hydroxy-dCTP ( Supplementary Fig. S3 ). In addition, T > C mutations were also prone to false detection by PCR or sequencing errors due to unmethylated C being indistinguishable from a single T > C sequencing error on the protected strand (Fig. 3B ), and we did observe the T > C (specifically, unmethylated C) mutation rate was close to the raw Illumina sequencing error rate (Fig. 3A ). For all single base contexts, except for C > T and T > C, the residual SNV frequencies exhibited minimal or no strand bias, indicating all Methyl-CODEC products can be used for their detection with duplex sequencing accuracy. Meanwhile, for C > T and T > C, duplex sequencing accuracy is still attainable when the complementary mutation (e.g. G > A or A > G) is detected, which will be present on one of two Methyl-CODEC products from each DNA duplex. Methyl-CODEC expands genome-wide mutation calling The conventional single base substitution alphabet, e.g. G > C, considers only mutations among the standard bases A, C, G, and T for both reference and alternative alleles. Methyl-CODEC extends this by distinguishing C or cytosine into non-methylated cytosine (umC) and 5mC as distinct alternative alleles, denoted as A > umC, G > umC, T > umC and A > 5mC, G > 5mC, T > 5mC, respectively (Fig. 3A ). In NA12878 Methyl-CODEC data, we found that five out of the six expanded mutation types have residual mutation frequencies close to or lower than those of the corresponding conventional mutation types in CODEC (Fig. 3A ), with T > umC being an exception, which we reason was dominated by sequencing errors on the protected strand (Fig. 3B ), but could be called with duplex accuracy if the other Methyl-CODEC product was sequenced. On the other hand, T > umC cannot be detected by standard enzymatic or bisulfite based methylation sequencing because unmethylated C is converted back to T by these methods. Interestingly, observed frequencies for A > 5mC, G > 5mC, and T > 5mC were significantly lower (between 10 −9 and 10 −7 ), aligning with our expectations, as 5mC predominantly occurs at CpG dinucleotides, which comprise only 1% of the genome. We then investigated the proportion of rare duplex mutations in NA12878 Methyl-CODEC that resulted in the formation of novel cytosine (7.8%) and CpG dinucleotides (2.3%) (Fig. 4A ), after removing C > T and T > umC due to their sub-duplex-sequencing accuracy. Next, we analyzed the trinucleotide contexts of the 148 rare duplex mutations creating 5mC and the 656 rare duplex mutations creating umC, and observed that the majority (82%) of the novel 5mC mutations occurred at newly formed CpG sites, particularly for T > 5mC and A > 5mC, whereas only 17% of the novel umC mutations formed at CpG sites (Fig. 4B ), consistent with established data and models indicating that methylation predominantly occurs at CpG sites in humans. Interestingly, for G > 5mC mutations, a significant portion (48.6%) did not occur at newly formed CpG sites. In addition, regardless of whether the newly formed CpGs are methylated or not, the overall methylation of existing neighboring CpGs remains consistent ( Supplementary Fig. S7 ). As methylation is assessed from only a single strand of each original DNA duplex using a single Methyl-CODEC product, it is possible that methylation accuracy could be further improved if both Methyl-CODEC products from each duplex were recovered. Nonetheless, our results suggest the feasibility of detecting rare mutations and their corresponding methylation status in single DNA duplexes, using single NGS read pairs, with Methyl-CODEC. Figure 4. Expanded mutation calling of Methyl-CODEC in NA12878. ( A ) A bar plot showing, from left to right, the total number of rare duplex SNVs, the number of rare duplex SNVs that produced cytosines, and the number of rare duplex SNVs that created new CpG sites, excluding C > T and T > umC, which are not at duplex accuracy. ( B ) Plots showing the count of rare duplex SNVs that produced a 5mC (left) versus umC (right) across trinucleotide contexts. Reference trinucleotide contexts that can generate novel CpG are indicated on the X -axis. T > umC is excluded due to lack of duplex accuracy. Discussion Methyl-CODEC enables simultaneous methylation and duplex sequencing using single read pairs. It introduces conversion-resistant adapters and dCTP analogs during CODEC strand linking and adds an enzymatic methyl conversion step along with a novel analysis workflow. While prior methods have been developed for combined genetic and epigenetic analysis [ 12–15 ], Methyl-CODEC is the first to achieve both genetic duplex sequencing accuracy and single-base methylation resolution using single read pairs at the whole genome level. To demonstrate this, we compared Methyl-CODEC’s genetic sequencing accuracy to standard CODEC [ 17 ] and its methylation sequencing performance to EM-seq WGMS and WGBS on the same cell line. Our results showed that Methyl-CODEC can achieve duplex sequencing accuracy for all single nucleotide contexts using single read pairs, including for C > T and T > C when G > A and A > G, respectively, are instead evaluated. For methylation sequencing, Methyl-CODEC demonstrated strong concordance with widely used techniques such as WGBS or EM-seq, while uniquely preserving the original DNA sequence for improved NGS alignment and C > T mutation detection. We found that Methyl-CODEC uniquely improved methylation detection at CpG sites (around 590 000 sites) overlapping germline C > T/G > A mutations. Methyl-CODEC also expanded genome-wide mutation calling by distinguishing C into umC and 5mC as distinct alternative alleles. This enabled us to discover genome-wide, low-frequency mutations and methylation co-occurring within single base pairs of single DNA duplexes. Those novel 5mC were found to be enriched in newly formed CpG contexts compared to non CpG contexts, suggesting that most of these novel 5mC sites were genuine and may reveal new biological processes. We believe this could be useful in many areas such as for cancer and minimal residual disease detection using liquid biopsies. A current limitation of Methyl-CODEC in epigenomic sequencing is its inability to distinguish between 5mC and 5hmC, although this could be added in the future. In addition, the current end-repair/dA-tailing step leads to artificial hypomethylation in the resynthesized regions, although we, among others, have described superior end-repair/dA-tailing methods that could limit resynthesis [ 26 , 27 ]. Furthermore, only about half of the coverage is usable for detecting C > T and T > C at duplex accuracy given the need to detect these by way of G > A and A > G, respectively. Although Methyl-CODEC generates two Methyl-CODEC products per original DNA duplex with each of the respective base changes, recovering both in a WGS library is rare. Methyl-CODEC, like CODEC, is conceptually compatible with targeted sequencing workflows. This has several conceptual advantages. For instance, by preserving a native strand sequence, Methyl-CODEC hybrid capture sequencing stands to streamline and mitigate potential biases typically associated with hybrid capture following bisulfite or enzymatic methyl-conversion. Lastly, Methyl-CODEC, like CODEC, requires short fragments (ideally ∼150 bp), whereas longer insert-size libraries can be sequenced using EM-seq or WGBS, potentially resulting in better alignment than observed in the current simulation. Methyl-CODEC has the potential to enable a wide range of clinical applications. In oncology, for example, it can integrate tumor-specific epigenetic and genetic markers to enhance comprehensive cancer detection and management, thereby advancing precision oncology. Duplex sequencing accuracy is pivotal for sensitive detection of minimal residual disease from blood [ 28 , 29 ] and Methyl-CODEC uniquely preserves this accuracy. By capturing epigenetic and genetic modifications in cell-free DNA (cfDNA), Methyl-CODEC could improve liquid biopsy-based cancer screening [ 30–33 ]. For example, while mutations on cfDNA can indicate the presence of cancer [ 34 ], epigenetic patterns detected on the same molecules can help identify the tissue of origin [ 35 ]. This capability can guide targeted imaging or biopsies to confirm the presence and extent of cancer prior to initiating therapy. Additionally, Methyl-CODEC can aid in the classification and selection of treatment strategies for cancers of unknown primary, a particularly challenging group to treat, by combining mutation with methylation analysis [ 36 , 37 ]. Beyond oncology, Methyl-CODEC could have applications in prenatal testing and transplantation, where it can identify the tissue of origin of DNA molecules derived from the fetus or tissue donor [ 38 ]. Methyl-CODEC could also be applied to other biological areas, such as detecting allele-specific methylation, which is crucial for localizing regulatory sequence polymorphisms [ 39 , 40 ] and assessing biological aging [ 41 , 42 ]. A recent study has shown that epigenetic clocks are influenced by somatic mutations, which can alter methylation patterns not only at the mutation site but also in surrounding regions [ 43 ]. Lastly, as Methyl-CODEC does not depend on a reference genome to derive methylation status, we envision applications in infectious diseases or in non-model organisms. In conclusion, by enabling simultaneous methylation and duplex sequencing using single read pairs, Methyl-CODEC opens promising new avenues of broad potential impact for biomedicine. Supplementary Material gkaf482_Supplemental_File