%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

\documentclass[doublespacing,linenumbers]{bmcart}
%\documentclass[doublespacing]{bmcart}

%%% Load packages
\usepackage{amsthm}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage[numbers,sort&compress]{natbib}
\usepackage[colorlinks=false]{hyperref}
\usepackage[utf8]{inputenc} %unicode support
\usepackage{graphicx}
\usepackage{caption}
\usepackage{booktabs}
%\usepackage{xcolor}

% NB: These line *must* be included when submitting to BMC.
\def\includegraphic{}
\def\includegraphics{}

%%% Put your definitions there:
\startlocaldefs
%    \newcommand\red[1]{{\color{red}#1}}
\endlocaldefs

\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
    \dochead{Research}

    \title{Evaluation of Whole-Genome DNA Methylation Sequencing Library Preparation Protocols}

    % Author list
    \author[ % Jacob Morrison
        addressref={aff1},
        email={jacob.morrison@vai.org}
    ]{\inits{J.}\fnm{Jacob} \snm{Morrison}}
    \author[ % Julie M Koeman
        addressref={aff2},
        email={julie.koeman@vai.org}
    ]{\inits{J.M.}\fnm{Julie M.} \snm{Koeman}}
    \author[ % Benjamin K. Johnson
        addressref={aff1},
        email={ben.johnson@vai.org}
    ]{\inits{B.K.}\fnm{Benjamin K.} \snm{Johnson}}
    \author[ % Kelly K. Foy
        addressref={aff1},
        email={kelly.foy@vai.org}
    ]{\inits{K.K.}\fnm{Kelly K.} \snm{Foy}}
    \author[ % Wanding Zhou
        addressref={aff3,aff4},
        email={wanding.zhou@pennmedicine.upenn.edu}
    ]{\inits{W.}\fnm{Wanding} \snm{Zhou}}
    \author[ % David W. Chesla
        addressref={aff5},
        email={david.chesla@spectrumhealth.org}
    ]{\inits{D.W.}\fnm{David W.} \snm{Chesla}}
    \author[ % Larissa L. Rossell
        addressref={aff5},
        email={larissa.rossell@spectrumhealth.org}
    ]{\inits{L.L.}\fnm{Larissa L.} \snm{Rossell}}
    \author[ % Emily J. Siegwald
        addressref={aff5},
        email={emily.siegwald@spectrumhealth.org}
    ]{\inits{E.J.}\fnm{Emily J.} \snm{Seigwald}}
    \author[ % Marie Adams
        addressref={aff2},
        corref={aff2},
        noteref={n1},
        email={marie.adams@vai.org}
    ]{\inits{M.}\fnm{Marie} \snm{Adams}}
    \author[ % Hui Shen
        addressref={aff1},
        corref={aff1},
        noteref={n1},
        email={hui.shen@vai.org}
    ]{\inits{H.}\fnm{Hui} \snm{Shen}}

    % Affiliation list
    \address[id=aff1]{ % Dept of Epigenetics, VARI
        \orgdiv{Department of Epigenetics},
        \orgname{Van Andel Research Institute},
        \street{333 Bostwick Avenue NE},
        \city{Grand Rapids},
        \state{MI},
        \postcode{49503},
        \cny{USA}
    }
    \address[id=aff2]{ % Genomics Core, VARI
        \orgdiv{Genomics Core},
        \orgname{Van Andel Research Institute},
        \street{333 Bostwick Avenue NE},
        \city{Grand Rapids},
        \state{MI},
        \postcode{49503},
        \cny{USA}
    }
    \address[id=aff3]{ % CHOP
        \orgdiv{Center for Computational and Genomic Medicine},
        \orgname{The Children's Hospital of Philadelphia},
        \street{3501 Civic Center Boulevard},
        \city{Philadelphia},
        \state{PA},
        \postcode{19104},
        \cny{USA}
    }
    \address[id=aff4]{ % Penn
        \orgdiv{Department of Pathology and Laboratory Medicine},
        \orgname{University of Pennsylvania},
        \city{Philadelphia},
        \state{PA},
        \postcode{19104},
        \cny{USA}
    }
    \address[id=aff5]{ % Spectrum
        \orgdiv{Spectrum Health Office of Research and Education},
        \orgname{Spectrum Health System},
        \street{15 Michigan Street NE},
        \city{Grand Rapids},
        \state{MI},
        \postcode{49503},
        \cny{USA}
    }

    % Any extra notes to add
    \begin{artnotes}
        \note[id=n1]{These authors are co-senior authors.}
    \end{artnotes}

\end{fmbox} % comment this for two column layout

\begin{abstractbox}
    \begin{abstract} % No more than 350 words
        \parttitle{Background} % context and purpose of study
        With rapidly dropping sequencing cost, the popularity of whole-genome DNA methylation sequencing has been on the rise. Multiple library preparation protocols exist, but a systematic evaluation and benchmarking of their performance against each other is currently lacking. We have performed 22 whole-genome DNA methylation sequencing experiments on fresh frozen human samples, and extensively benchmarked common library preparation protocols for whole-genome DNA methylation sequencing, including three traditional bisulfite-based protocols and a new enzyme-based protocol. Additionally, different input DNA quantities were compared for two kits compatible with a reduced starting quantity. In addition, we also present bioinformatic analysis pipelines for sequencing data from each of these library types.

        \parttitle{Results} % main findings
        An assortment of metrics were collected for each kit, including raw read statistics, library quality and uniformity metrics, cytosine retention, and CpG beta value consistency between technical replicates. Overall, the NEBNext Enzymatic Methyl-seq kit performed quantitatively better than the other three protocols at two different DNA input amounts. Additionally, the results for the different input amounts were generally consistent across all metrics.

        \parttitle{Conclusions} % brief summary and potential implications
        Based on these results, we recommend use of the NEBNext Enzymatic Methyl-seq kit for whole-genome DNA methylation sequencing. Further, a general bioinformatic pipeline is applicable across the four protocols, with the exception of extra trimming needed for the Swift Biosciences's Accel-NGS Methyl-Seq protocol to remove the Adaptase sequence.
    \end{abstract}

    \begin{keyword} % keywords
        \kwd{DNA Methylation}
        \kwd{Epigenetics}
        \kwd{Whole genome bisulfite sequencing}
        \kwd{Enzymatic methylation sequencing}
        \kwd{Fallopian tube}
    \end{keyword}
\end{abstractbox}

%\end{fmbox}% uncomment this for two column layout

\end{frontmatter}

\section*{Background}
The dynamic interplay between various epigenetic modifications influence cellular differentiation, lineage specification, tissue development, and can also promote oncogenic states through changes in histone modifications and DNA methylation. In mammals, DNA methylation typically occurs at the 5' position of CpG dinucleotides (denoted 5-mC) and remains the single best studied epigenetic mark, partly due to its robustness through most storage conditions and histological preparations, such as formalin-fixed, paraffin-embedded (FFPE) samples \cite{laird_2003}.

DNA methylation can be assayed with various principles \cite{laird_meth_analysis}, such as bisulfite conversion \cite{clark_1994}, restriction enzyme digestion, or differential affinity for methylated DNA binding proteins. Array-based and sequencing-based methods built upon these principles have been developed and benchmarked \cite{bock_2010,harris_2010,bock_2016}. The current gold standard approach to examine genome-wide DNA methylation composition and differences is through chemical modification of unmethylated cytosines using sodium bisulfite \cite{clark_2006}. Bisulfite deaminates unmethylated cytosines (Cs) to uracils that are converted to thymines (Ts) during PCR amplification. Methylated cytosines (mCs) remain unaltered through this process. The end result yields stable genetic differences between methylated (C) and unmethylated cytosines (T), reflecting the underlying DNA methylation landscape, effectively turning the epigenetic difference into a genetic difference, which can then be studied using conventional genome-scale methods, such as microarrays or sequencing. Various generations of bisulfite-based microarrays have been used to profile hundreds of thousands of human samples due to the low cost and easy, standardized data processing and analysis. With dropping sequencing cost in recent years, the popularity of sequencing-based methods has been on the rise \cite{zhou_nature_genetics}. 

Whole-genome bisulfite sequencing (WGBS) provides the most comprehensive single base resolution DNA methylation maps. It was successfully applied to \textit{Arabidopsis thaliana} in 2008 \cite{lister_2008, cokus_2008} and then to humans in 2009 \cite{lister_2009}. In these early methods, adapter-ligated library material undergoes bisulfite conversion, leading to sheared and degraded library fragments and overall lower quantities and diversity of sequenceable material. A post-bisulfite adapter tagging (PBAT) method \cite{pbat_protocol_nar, pbat_method} was developed to overcome this hurdle, effectively decreasing the input range to nanogram level. Notably, this method has been used for single-cell WGBS profiling \cite{smallwood_scwgbs}. However, PBAT libraries typically do not maintain strand specificity, thus leading to ambiguity in quantifying strand specific methylation states, such as imprinted regions, which are critical for allele-specific gene expression patterns \cite{tucci_2019}. More recently, Swift Biosciences has produced a kit that is an alternative approach to the post-bisulfite library preparation. The alternative approach maintains the low DNA input capabilities of PBAT, while also adding strand specificity by including a low-complexity sequence on the 3’ end of the ssDNA during library preparation that serves as a scaffold for sequencing adapter attachment (Accel-NGS Methyl-Seq protocol, Swift Biosciences).

The conditions needed for bisulfite conversion are known to be harsh on the DNA and cause degradation. In recent years, it has become clear that this conversion can also be achieved with an enzymatic approach. 5-mCs can be converted to 5-hydroxymethylcytosine (5-hmC), then to 5-formylcytosine (5-fC), and eventually to 5-carboxylcytosine (5-caC) by the ten-eleven translocation (TET) family dioxygenases \cite{em_seq}. Further, the apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like 3A (APOBEC3A) deaminates methylated and unmethylated Cs into thymines and uracils, respectively \cite{em_seq}. While both 5-mCs and Cs are affected by this process, the TET-oxidized methylcytosines, 5-hmC, 5-fC, and 5-caC are minimally impacted \cite{em_seq}. Based on this principle, an enzymatic methyl-seq (EM-seq) method was recently developed by New England Biolabs. Their method uses TET2 to oxidize methylated cytosines and subsequent APOBEC3A treatment to convert unmethylated cytosines to uracils, while maintaining strand specificity \cite{em_seq}. WGBS and EM-seq, collectively referred to as whole-genome methylation sequencing (WGMS), both convert a 5-mC/C difference to a C/T difference; therefore, analysis tools developed for WGBS are also applicable to EM-seq.

In this study, we extensively benchmark the performance of three most commonly used protocols for bisulfite-based whole-genome DNA methylation profiling including the KAPA Hyper Prep kit (Kapa), the Miura and Ito post-bisulfite adapter tagging (PBAT) method, and the Swift Biosciences Accel-NGS Methyl-Seq DNA library kit (Swift), as well as the new EM-seq protocol from New England Biolabs, the NEBNext Enzymatic Methyl-seq kit (NEB), on fresh frozen human tissue samples. For each technique, we evaluate input quantities, read mapping statistics, library complexity, insert size, cytosine retention, as well as reproducibility between replicates. We also present bioinformatic analysis pipelines for each of these library types.

\section*{Results}
Benchmarking studies often use cell lines due to the largely isogenic background and reduced biological variance to probe the reproducibility of the method. However, it is often of interest to apply these approaches to complex tissue sources for primary research. Thus, we have devised an approach that uses frozen normal primary human solid tissue to benchmark these kits in a more common research scenario. The results presented are based on two primary patient fallopian tube samples (denoted A and B). DNA from each sample was prepared using one of the four library preparation protocols, as described in Table \ref{tab:data_summary}. With the exception of the PBAT protocol, the DNA from each protocol was processed twice to produce technical replicates (denoted 1 and 2). In addition to generating samples using the suggested amount of DNA input, a smaller DNA input (10 ng each) was used in the NEB and Swift protocols to test the effectiveness of low DNA inputs on these protocols. (New England Biolabs states they can go as low as 10 ng \cite{neb_technote}, while Swift Biosciences states they can go as low as 100 pg \cite{swift_protocol}.) In addition to summarizing the samples used in this analysis, Table \ref{tab:data_summary} shows which protocols generate strand-specific libraries.

The first set of metrics compared between the four preparations is related to the quality of raw reads received from the sequencer, including base quality and adapter contamination, and the effects of trimming on those reads. In general, the raw base qualities, percentage of reads with adapter contamination, and percentage of bases trimmed are comparable between the Kapa, NEB, and Swift protocols (Figure \ref{fig:raw_trim} and Additional file 1: Figure S1). However, the PBAT protocol suffers from a higher percentage of low-quality bases along the length of the read, leading to a higher percentage of trimmed bases during the trimming stage. Because the PBAT protocol does not contain an amplification step during library preparation, the higher percentage of trimmed bases relative to the other preparations has a larger effect on the amount of usable data available from this protocol.

Following adapter and quality trimming of the raw reads, the reads were mapped to the human genome and the fraction of optimally aligned (defined as MAPQ $\geq40$), sub-optimally aligned (MAPQ $<40$), and not aligned read fragments was calculated. Figure \ref{fig:lib_quality}A shows, for Sample A, the NEB and Swift protocols had about the same fraction of read fragments that were optimally aligned ($\sim85$\%), regardless of the amount of input DNA. The Kapa protocol was slightly behind ($\sim80$\%), while the PBAT protocol was closer to $\sim75$\%. The lower percentage of optimally aligned reads for the PBAT protocol, when coupled with the substantially lower number of read fragments, means there is much less data to use when performing analyses using this protocol.

Another metric that can speak to the quality of a library preparation is the insert size, which relates to the size of the sequenced DNA fragments. For each Kapa, Swift, and NEB sample, DNA was sheared as a single sample. The shearing process should leave each sample with the same starting length of DNA, but it can be seen in Figure \ref{fig:lib_quality}B that the bisulfite preparations (Kapa and Swift) retained shorter fragments. On the other hand, the enzymatic preparation (NEB) retained a wider range of fragment lengths, which are generally longer than the bisulfite fragments. The retention of a smaller range of shorter fragments by the bisulfite preparations suggest degradation of the DNA during library generation, which is consistent with the known tendency of bisulfite conversion to degrade samples. Whereas the other samples were sheared prior to conversion, the PBAT samples were only sheared by the bisulfite process itself. Figure \ref{fig:lib_quality}B shows this process is quite destructive, leaving shorter fragments than the other bisulfite conversion protocols.

With regards to the fraction of reads with MAPQ $\geq40$ marked as duplicates (Figure \ref{fig:lib_quality}C), the PBAT protocol and standard input NEB and Swift protocols each had about 10\% of duplicate reads. The low input NEB protocol was slightly higher at $\sim15$\%. The Kapa and low input Swift protocols had the highest percentage of reads, sitting closer to 25\%.

The library complexity, seen in Figure \ref{fig:lib_quality}D for Sample A, is a metric that can be used to determine if a library has reached a saturation point, where sequencing deeper will only gain a marginal amount of unique (i.e. not duplicate) reads. While none of the samples were sequenced to saturation, the Kapa and low input Swift samples showed a lower level of complexity compared with the NEB and standard input Swift samples. Due to the PBAT samples having much fewer reads, it is difficult to ascertain where the PBAT library complexity ranks compared to data derived from the NEB or Kapa protocols at higher read depths. At the PBAT sequencing depth in this study, it appears to have a trend similar to the Kapa protocol complexities, implying an overall lower level of library complexity relative to Swift and NEB data.

In order to determine the uniformity of reads distributed across the genome, the ratio of the observed read coverage to the expected coverage was calculated across several regions (Figures \ref{fig:cpg_coverage}A--\ref{fig:cpg_coverage}B and Additional file 1: Figure S6). In general, there is consistent coverage across the genic, intergenic, and repeat-masked regions, with each sample having less than a 5\% departure from expected uniformity (closer to 1.0 is better) in these regions, as shown in Figure S6a--S6c. In contrast, exonic regions, all CpGs, and CpG islands show greater heterogeneity across kits and larger departures from uniform coverage (Figures \ref{fig:cpg_coverage}A-- \ref{fig:cpg_coverage}B and Figure S6d). Here, the PBAT protocol does well, having a much higher observed rate of coverage compared to expected. The other protocols all have lower coverage than would be expected, with the NEB samples having the closest to uniform coverage and the Kapa samples having the lowest coverage of these regions. The undercoverage of CpGs by the Kapa protocol is consistent with Nair and Luu \textit{et al.} \cite{nair_2018}. One example of the differences in coverage uniformity across the protocols can be found in the EPCAM promoter region, as seen in Figure \ref{fig:epcam_coverage}.

When performing WGMS, a library preparation's coverage of cytosines, particularly in a CpG context, is important to assess the DNA methylation landscape of a given sample. When comparing the protocols used in this analysis, the Kapa protocol has the lowest percentage of CpGs covered in a number of different regions on these samples (Figures \ref{fig:cpg_coverage}C--\ref{fig:cpg_coverage}D and Additional file 1: Figure S7), which is consistent with previous work \cite{nair_2018}. Generally, the other preparations have coverage over 85\%, with the exception of the Swift samples' coverage of CpG islands, which is in the 75-80\% range. It should be noted that the low DNA input runs of the NEB protocol produced coverage levels that are consistent with the standard DNA input runs across samples and technical replicates.

The read-averaged cytosine retention for all protocols reflect what would be expected for a WGMS run, namely CpG retention above 70\% and the other cytosine contexts around 1\% (Figure \ref{fig:read_avg_cret}). Moreover, the technical replicates showed consistent retention, while the two biological replicates showed a bigger difference. Across all protocols, Sample A consistently had higher CpA retention, but not CpC or CpT retention. It is known that biological CpH retention tends to occur in the CpA dinucleotide context \cite{schultz_2015}. The higher CpA retention, in contrast to CpC and CpT retention, shows Sample A likely has true CpH methylation, a process previously thought to be largely restricted to embryonic stem cells and neurons \cite{he_2015}. Our results also confirm that, unlike CpA methylation, CpC and CpT retention likely do not reflect true biological methylation, but incomplete conversion, at least in mammalian samples. The Swift and NEB protocols generally showed the lowest amount of such technical artifact (both below 0.5\%), except for one replicate of the low input Swift preparation. Similar results can be seen for the base-averaged cytosine retention, which can be found in Additional file 1 Figure S8.

The M-bias plots \cite{m_bias} (Additional file 1: Figure S9) show a consistently lower CpG retention across the entire read length for both reads 1 and 2 of the PBAT protocol. For the Kapa protocol, there is a consistently higher CpG retention rate for read 1 than the other preparations. However, read 2 tends to be more in line with the retention rate seen in the others. For CpH retention, each protocol has approximately the same level of retention, with slight deviations in the first 5 bp. Due to the adapter trimming performed (see Methods section for details), retention rates can behave erratically at the end of reads depending on the base content of the adapter that is trimmed. It should also be noted that, by default, the aligner used in this analysis (see Methods section) does not include cytosines in the first and last 3 bp of a read when determining methylation beta values, so erratic behavior on the ends of the reads seen in Figure S9 is not included in methylation-related metrics.

In order to compare the consistency of CpG beta values, Spearman correlation coefficients were calculated between technical replicates of each protocol (Figure \ref{fig:tech_corr}), as well as between preparations of the same sample (Additional file 1: Table S1). The correlations between the NEB technical replicates, both the standard and low input samples, had the highest correlation at 0.9 or higher. The Kapa protocol had the lowest correlation (i.e. less consistency in beta values between replicates) of just under 0.75. The standard input Swift sample did better than the low input sample, with correlations of 0.873 and 0.814, respectively. A correlation for the PBAT protocol could not be calculated because PBAT did not have a technical replicate. When comparing preparations to one another, the Kapa protocol had low correlations, with the maximum correlation of 0.62 to the Swift protocol (Table S1). Only two of the other ten correlations were 0.8 or below, both of which were between the PBAT and Swift protocols. Of the NEB and Swift protocols, the NEB protocol had the highest cross-protocol correlations, with overall greater correlations compared to the Swift protocols.

Because of the different approach to cytosine conversion used by the NEB protocol relative to others (enzymatic versus chemical), the difference in beta values between the standard input NEB and Swift protocols were compared to look for bias in the methylation level of the NEB protocol. After calculating the beta value differences for Sample A, only four CpGs with $|diff| > 0.5$ for both Technical Replicate 1 and 2 were found (see Additional file 1: Figure S10). However, in Sample B, only one (which is not one of the four CpGs in Sample A) was seen, so this difference was taken to be sample dependent and not due to an inherent bias in the enzymatic conversion process.

\section*{Discussion}
After comparing the four different WGMS library preparations, we noted that the NEBNext Enzymatic Methyl-seq kit (EM-seq) performed quantitatively better in almost every metric, both with a standard amount (200 ng) and a low amount (10 ng) of input DNA. The NEB samples had higher quality libraries, with larger insert sizes and higher library complexity than the other kits, while having at least a comparable, if not more, uniform distribution of reads across the genome. The larger insert sizes lead to better mappability during alignment (i.e., a higher fraction of optimally aligned reads). Moreover, it also facilitates read-level analysis, such as epi-polymorphism \cite{landan_polymorphism} and allele-specific DNA methylation analyses. The higher library complexity implies more unique information could be gleaned from deeper sequencing of these libraries. Furthermore, the low input NEB samples showed little reduction in library complexity relative to the standard input NEB samples. Because of this, it may be possible to push the lower threshold of DNA input below New England Biolabs's recommended limit of 10 ng. When comparing methylation-related metrics, the NEB samples were among the highest in percentage of CpGs covered, with consistent cytosine retention across all reads and a sub-one percent CpH retention level. Among traditional bisulfite-based methods, the Swift Biosciences Accel-NGS Methyl-Seq DNA library kit performed best, which is consistent with a previous study that compared Swift with two other bisulfite-based protocols, TruSeq and QIAseq \cite{zhou_2019}. Its results were comparable enough to the NEB protocol, particularly at the standard input level (100 ng), such that conclusions drawn from both protocols should be fairly consistent. However, its performance trailed as the input amount dropped. 

All kitted protocols had comparable operability in the lab. Kapa and Swift took roughly 7 hours to process, which can fit in a single workday. The NEB enzymatic protocol was a 9-hour protocol but had several convenient stopping points noted to easily break processing into two days. Actual hands-on time for all three protocols was about 4.5 hours. The PBAT protocol was the only technically difficult and time consuming one, mostly because of the requirement to make and quality-control solutions. This protocol took between 14 and 16 hours to complete, about 7 of which were hands on.

One notable difference separating the four DNA methylation library generation methods is library uniformity, particularly in CpG island regions (CGIs). Kapa library coverage was underrepresented and PBAT overrepresented in these regions, while Swift and NEB generated relatively uniform libraries (Figures \ref{fig:cpg_coverage} and \ref{fig:epcam_coverage}). We propose that mechanistic differences in library generation account for this phenomenon, with both the relative timing of bisulfite conversion and the choice of random priming versus adapter ligation affecting the CGI coverage in the downstream dataset.  In the Kapa protocol, sheared and adapter-ligated libraries undergo bisulfite conversion, causing preferential strand breakage in CGI regions and inhibiting post-bisulfite conversion PCR steps, as both priming regions are no longer present on the library fragment.  In the PBAT protocol, libraries are not sheared. Instead, the bisulfite conversion induces random DNA breaks \cite{munson_2007,tanaka_2007,grunau_2001,mill_2009}, again more prevalently in the CGI regions, making these small fragments more available for continued library generation than larger fragments.  Additionally, CGIs are G:C rich; during the random priming step of the PBAT library creation process, these higher G:C regions generate more stable primer annealing conditions and increase the likelihood of extension and incorporation into the final library due to a slightly higher regional melting temperature and exhibition of G:C clamping features \cite{dieffenbach_1993}. The more even CGI coverages exhibited by the Swift and NEB libraries are also protocol related: in the Swift protocol, bisulfite conversion precedes manual shearing and adapter ligation, avoiding the strand-breakage issue and generating relatively diverse libraries, and NEB’s process completely avoids the deleterious bisulfite conversion step through enzymatic reactions.

The clinical utility of DNA methylation has been long recognized, explored, and established, particularly for cell-free DNA (cfDNA)-based early detection of neoplasm. The current gold-standard technologies for DNA methylation profiling are bisulfite-based. Harsh bisulfite treatment is known to cause heavy degradation of the template, which is often scarce to begin with for clinical cfDNA samples. cfMeDIP-seq \cite{shen_2018} has been developed to overcome this limitation, pushing the lower input limit to 1-10 ng. However, it is associated with similar limitations of common affinity-based methods \cite{laird_meth_analysis,shen_2018}. As the enzymatic conversion by the NEB protocol still demonstrates a library complexity at 10ng comparable to that of regular bulk DNA-based methods, the NEB protocol may prove to be a good alternative approach for clinical early detection work. In addition, with its template-preserving nature and excellent performance at low-input levels, the NEB protocol could hold more promise for single-cell DNA methylation profiling.

\section*{Conclusion}
In this study, we compared four commonly used WGMS library preparation protocols, including three bisulfite-based protocols and one enzyme-based protocol. We found the NEBNext Enzymatic Methyl-seq kit performed quantitatively better than the other three protocols at two amounts of DNA input (200 ng and 10 ng). Furthermore, the results from the two amounts of DNA input were comparable across biological and technical replicates. Based on these results, we recommend use of the NEBNext Enzymatic Methyl-seq kit for whole-genome DNA methylation sequencing.

\section*{Methods}
\subsection*{Fallopian Tube Sample Preparation}

Fallopian tubes from two primary patients were delivered in saline from a local hospital. Upon arrival, samples were washed using sterile 1x Phosphate Buffered Saline and snap frozen in liquid nitrogen. At a later date, the fallopian tubes were thawed and minced. DNA was extracted from the minced fallopian tubes using the tissue protocol from Qiagen's DNeasy Blood \& Tissue Kit (69504). Following extraction, dsDNA was quantified using Invitrogen’s Qubit 3.0 Fluorometer.

\subsection*{Whole-Genome Methylation Sequencing Libraries}

DNA controls pUC19 and unmethylated lambda DNA (0.0005\% and 0.01\%, respectively) were added to each high molecular weight genomic DNA sample. The DNA was then sheared to approximately 350 bp in average size for all prepared libraries, with the exception of the post-bisulfite adapter tagging (PBAT) libraries, where the DNA was not initially sheared.

Libraries were prepared from the KAPA Hyper Prep kit (v6.17) (Roche Sequencing, Cat. \#KK8504) with an input of 300 ng of sheared DNA following the manufacturer's protocol with the following modifications. Illumina TruSeq Nano adapters at a concentration of 10 $\mu$M were used. The post ligtation cleanup elution was reduced to 20 $\mu$L and the entire DNA elution went into the EZ DNA Methylation-Gold kit (Zymo Research, Cat. \#D5005). The bisulfite converted DNA was eluted in 20 $\mu$L and 10 cycles of library amplification were performed using the KAPA HiFi HotStart Uracil+ ReadyMix (Roche Sequencing, Cat. \#KK2800).

Libraries were prepared from the Accel-NGS Methyl-Seq DNA library kit (v3.0) (Swift Biosciences, Cat. \#30024) with an input of either 10 ng or 100 ng of sheared DNA following manufacturer's protocol with the following modifications. The DNA was bisulfite converted using the EZ DNA Methylation-Gold kit (Zymo Research, Cat. \#D5005) with an elution volume of 15 $\mu$L. Following adapter ligation, either 8 cycles (10 ng DNA input) or 4 cycles (100 ng DNA input) of library amplification were performed.

Libraries were prepared from the NEBNext Enzymatic Methyl-seq kit (New England Biolabs, Cat. \#E7120L) using an input of either 10 ng or 200 ng of sheared DNA and libraries were made according to the manufacturer's protocol. The denaturation method used was 0.1 N sodium hydroxide, according to the protocol, and either 8 cycles (10 ng DNA input) or 4 cycles (200 ng DNA input) of PCR amplification were performed.

Libraries were prepared from the PBAT method described by Miura and Ito \cite{pbat_method} using an input of 100 ng of genomic DNA that went directly into the EZ DNA Methylation-Gold kit (Zymo Research, Cat. \#D5005) with an elution volume of 21 $\mu$L. Then, 10 $\mu$L of the bisulfite converted DNA was used for making the PBAT libraries as previously described in \cite{pbat_method} with the modification of using KAPA HiFi HotStart Uracil+ ReadyMix (Roche Sequencing, Cat. \#KK2800) for the DNA template extension step.

KAPA pure beads (Roche Sequencing, Cat. \#KK8001) were used for cleanup steps for all prepared libraries.

Quality and quantity of the finished libraries were assessed using a combination of the Agilent High Sensitivity DNA chip (Agilent Technologies, Inc., Cat. \#5067-4626), QuantiFluor\textsuperscript{\textregistered} dsDNA System (Promega Corp., Cat. \#E2670), and Kapa Illumina Library Quantification qPCR assay (Roche Sequencing, Cat. \#KK4824). 100 bp paired-end sequencing was performed on an Illumina NovaSeq6000 sequencer using an S4, 300 bp sequencing kit (Illumina Inc., San Diego, CA, USA), with 10\% PhiX. Base calling was done by Illumina RTA3 and output of NCS was demultiplexed and converted to FASTQ format with Illumina Bcl2fastq (v1.9.0).

\subsection*{Alignment and Methylation Extraction}

Upon receipt of the FASTQ files, the files were trimmed using Trim Galore \cite{trim_galore} version 0.6.4\_dev (using Cutadapt version 2.10). Default inputs were used, other than: \texttt{--illumina --trim-n --paired --cores 4 --fastqc --fastqc\_args "--noextract"}. Additionally, the FASTQ files for the Swift samples included \texttt{--clip\_R2 14}, due to the Adaptase\textsuperscript{\texttrademark} method used by Swift Biosciences \cite{swift_trimming}.

The trimmed FASTQ files were aligned to GRCh38 \cite{gencode} using BISCUIT \cite{biscuit} version 0.3.16. An index for the reference genome was created using \texttt{biscuit index GRCh38.p13.genome.fa}, followed by aligning each sample to the indexed reference. The alignment step used the default options for \texttt{biscuit align}, with these exceptions: \texttt{-M -t 20 -R sample\_specific\_read\_group}. Each sample received its own read group (\texttt{-R} tag). The aligned reads were duplicate marked using Samblaster \cite{samblaster} version 0.1.25, with the \texttt{-M} flag and defaults. The reads were then coordinate sorted and indexed using Samtools \cite{samtools} version 1.10 (with htslib 1.10.2). Default options were used, with the following exceptions \texttt{-@ 20 -m 5G -o sample\_name.sorted.markdup.bam -O BAM} (sort) and \texttt{-@ 20} (index).

The extraction of cytosine methylation information proceeded as follows. Pileup VCF files were generated using \texttt{biscuit pileup} with default parameters. \texttt{bgzip} and \texttt{tabix} (included with htslib version 1.10.2) were used to compress and index the VCF files. Default parameters were used for \texttt{bgzip} and \texttt{tabix}, with the exception of \texttt{-p vcf} in the call for \texttt{tabix}. The VCF files were then processed through \texttt{biscuit vcf2bed}, \texttt{bedtools sort} (bedtools \cite{bedtools} version 2.29.2), and \texttt{biscuit mergecg} to create coordinate sorted BED files containing CpG methylation beta value information. For each command, the default parameters were used. After creating the BED files, they were compressed and indexed using \texttt{bgzip} and \texttt{tabix}, with default parameters being used in both cases, with \texttt{tabix} also including \texttt{-p bed}.

Any code not explicitly stated in this and subsequent sections can be found on GitHub \cite{source_code}.

\subsection*{FASTQ and Alignment Quality Control}

Quality control data were collected from a number of sources and viewed using MultiQC \cite{multiqc} version 1.9. Quality control for the FASTQ files was generated by FastQC \cite{fastqc} version 0.11.9 and Cutadapt during the trimming process. Command arguments used can be found in the previous section. Statistics on the percentage of duplicate marked reads were produced by Samblaster. Library complexity was calculated by Preseq \cite{preseq} version 2.0.3 using these options, \texttt{c\_curve -B -P -v -o sample.complex.ccurve.txt}, where ``sample'' is the name of each sample. Quality controls from Samtools were generated with \texttt{samtools stats} and \texttt{samtools flagstat}, with default parameters and \texttt{-@ 20}. BISCUIT includes a custom BASH script to generate quality control statistics related to data aligned by BISCUIT. This script was run with this command, \texttt{QC.sh -v samp.pileup.vcf.gz -o samp\_QC hg38\_assets GRCh38.p13.genome.fa samp samp.sorted.markdup.bam}. In each case, ``samp'' corresponds to the name of the processed sample.

The hg38\_assets mentioned in the BISCUIT quality control script command can be found in a zip file on the BISCUIT GitHub release page \cite{qc_assets}.

\subsection*{Library Protocol Comparison Analysis}

In order to collect statistics related to the raw reads stored in the FASTQ files, a custom Python (version 3.7.6) script was written. It uses the \texttt{gzip}, \texttt{glob}, and \texttt{time} Python base modules and these additional Python packages: \texttt{argparse} (version 1.1), \texttt{numpy} (version 1.18.1), and \texttt{biopython} (version 1.76). Statistics regarding the trimmed reads were collected from log files generated by Cutadapt, as described previously.

For a number of the analyses, the aligned BAM files were subsampled before calculating the corresponding metric. The BAMs were subsampled using \texttt{samtools view -hbu -F 0x4 -q 40 sample.bam | samtools view -hbu -s FRAC -}. ``FRAC'' was calculated as
\begin{equation}
\frac{150,000,000}{Number~Mapped~Fragments~with~MAPQ\geq40},
\end{equation}
for each sample. The subsampled BAMs were sorted and indexed using Samtools. Pileup VCF files and merged CpG BED files were generated in a similar manner to the original BAMs, as described previously.

Using the subsampled BAMs, the average coverage and percentage of covered CpGs within different genomic regions, including CpG islands, exons, genes, and repeat-masked regions, were calculated using custom scripts. The CpGs that fell in each region were determined by intersecting a BED file containing CpG coordinates and the coverage at those locations with a BED file containing the region's coordinates. The coverage was determined using \texttt{bedtools genomecov}, while the intersection was done using \texttt{bedtools intersect}. The average coverage was calculated by taking a weighted average of the coverage for each CpG. The percentage of covered CpGs was calculated by
\begin{equation}
\frac{Number~of~CpGs~in~Region~with~Coverage>1}{Total~Number~of~CpGs~in~Region}.
\end{equation}
The scripts used to calculate these values made use of GNU parallel \cite{gnu_parallel}.

The observed coverage to expected coverage ratio was calculated as $ratio = Obs / Exp$, where:
\begin{align}
Obs &= \frac{Number~of~Bases~Mapped~to~Region}{Total~Number~of~Mapped~Bases} \\
Exp &= \frac{Sum~of~Mappability~Scores~in~Region}{Total~Sum~of~Mappability~Scores~in~Genome}.
\end{align}
This formulation assumes all bases are not equally accessible when sequencing. The expected coverage takes into account this difference in accessibility by including mappability scores based on the Bismap k100 multi-read mappability scores \cite{bismap}. The observed coverage does not include these scores, as mappability is assumed to be inherently included when performing DNA sequencing. Because Bismap did not include a mappability score for every base in the genome, the expected and observed coverage calculations were restricted to those bases that included a mappability score. Only mapped reads (FLAG field in BAM does not include \texttt{0x4} flag) with MAPQ score $\geq 40$ were included in this calculation. The values for each observed to expected ratio were calculated using custom BASH scripts that used \texttt{bedtools} and \texttt{GNU awk} (version 4.0.2).

In order to calculate correlations between samples, the genome was binned into 100 kb bins and the average beta value calcuated for CpGs in each window. The bins were determined via \texttt{bedtools makewindows -w 100000 -g GRCh38.p13.genome.fa.fai | sort -k1,1 -k2,2n}. The ``.fai'' file was generated via \texttt{samtools faidx GRCh38.p13.genome.fa}. The average beta value for each bin was calculated via \texttt{bedtools map -a bins.bed -b sample.bed -c 4 -o mean | gzip}. \texttt{bins.bed} contains the 100 kb bins, while \texttt{sample.bed} is the merged CpG BED file generated from the subsampled BAMs. Correlations between the 100 kb averaged beta values were the calculated using the Spearman correlation coefficient, as implemented in Python's \texttt{scipy} package (version 1.4.1).

The analysis for methylation bias in the NEBNext Enzymatic Methyl-seq kit was performed by extracting CpGs that had more than 20 reads covering them for both the NEB and Swift protocols for Sample A replicate 1 and 2 or the NEB and Swift protocols for Sample B replicate 1 and 2. Methylation bias was determined by the requiring both Equations \ref{eqn:bias1} and \ref{eqn:bias2} to be true.
\begin{align} \label{eqn:bias1}
\lvert\beta_{NEB} - \beta_{Swift}\rvert &> 0.5~(Replicate~1) \\
\lvert\beta_{NEB} - \beta_{Swift}\rvert &> 0.5~(Replicate~2) \label{eqn:bias2}
\end{align}
This was done separately for the standard DNA input NEB and Swift protocols of Samples A and B.

The CpG and other genomic region BED files mentioned in this section can be found on the GitHub release page for the analysis code.

\begin{backmatter}

\section*{Additional Files}
	\subsection*{Additional file 1}
    \textbf{Figure S1:} Raw read statistics for each protocol. (a)--(e) Statistics for read 1. (f)--(j) Statistics for read 2.
    \textbf{Figure S2:} Library quality metrics for each protocol for Sample B. (a) The percentage of optimally, sub-optimally, and not aligned read fragments for each protocol. Note, read fragments treat reads 1 and 2 as separate entities, as it is possible that one read in the pair is mapped, while the other is not. (b) Insert size distribution. (c) The library complexity, which is a function of the duplicate rate.
    \textbf{Figure S3:} Bioanalyzer traces are shown for each sample.
    \textbf{Figure S4:} The distribution of MAPQ scores for Sample A and Sample B.
    \textbf{Figure S5:} The cumulative coverage across CpGs and all base for reads with MAPQ $\geq$ 40. (a)--(b) Cumulative coverage for Sample A. (c)--(d) Cumulative coverage for Sample B. Note, all libraries were downsampled to be comparable to PBAT (150 million reads, see Methods for details); therefore, any differences are not confounded by sequencing depth.
    \textbf{Figure S6:} The ratio of observed coverage to expected coverage for various genomic element categories. (a) Ratio for genic regions. (b) Ratio for intergenic regions. (c) Ratio for repeat-masked regions. (d) Ratio for exonic regions.
    \textbf{Figure S7:} Percentage of CpGs covered by at least one unique read with MAPQ $\geq40$ for various genomic element categories. (a) Percentage of covered CpGs in exonic regions. (b) Percentage of covered CpGs in genic regions. (c) Percentage of covered CpGs in repeat-masked regions. Note, all libraries were downsampled to be comparable to PBAT; therefore, any differences are not likely confounded by sequencing depth.
    \textbf{Figure S8:} Base-averaged cytosine retention by dinucleotide context, namely (A) CpA, (B) CpC, (C) CpG, and (D) CpT. In each panel two technical replicates are shown for each biological replicate. The x-axis denotes percent retention, with a scale of 0-5\% for CpH panels and 0-100\% for the CpG panel.
    \textbf{Figure S9:} CpG and CpH retention by read position for each biological sample. (a)--(b) The CpG retention by read position in Sample A. (c)--(d) The CpH retention by read position in Sample A. (e)--(f) The CpG retention by read position in Sample B. (g)--(h) The CpH retention by read position in Sample B. The length of read 2 in the Swift samples is shorter due to additional trimming required for these libraries (see Methods section for details). Note, all libraries were downsampled to be comparable to PBAT; therefore, any differences are not likely confounded by sequencing depth.
    \textbf{Table S1:} Spearman correlation coefficients between protocols for Sample A Technical Replicate 1. Note, all libraries were downsampled to be comparable to PBAT; therefore, any differences are not likely confounded by sequencing depth. Overall low correlation values are due to low coverage from downsampling.
    \textbf{Figure S10:} Difference in beta values between the NEB and Swift protocols. Red x's denote CpGs where the difference for both technical replicates exceeds 0.5.

\section*{Acknowledgements}
Computation for the work described in this paper was supported by the High Performance Cluster and Cloud Computing (HPC3) Resource at the Van Andel Research Institute.

\section*{Funding}
This study was funded by the National Institutes of Health/National Cancer Institute [R37 CA230748] and the Ovarian Cancer Research Fund (now Ovarian Cancer Research Alliance) [373933].

\section*{Abbreviations}
\textbf{DNA:} Deoxyribonucleic acid;
    \textbf{PCR:} Polymerase chain reaction;
    \textbf{FFPE:} Formalin-fixed, paraffin-embedded;
    \textbf{WGBS:} Whole-genome bisulfite sequencing;
    \textbf{EM-seq:} Enzymatic methyl-seq;
    \textbf{WGMS:} Whole-genome (DNA) methylation sequencing;
    \textbf{PBAT:} Post-bisulfite adapter tagging;
    \textbf{C:} Cytosine;
    \textbf{T:} Thymine;
    \textbf{mC:} Methylated cytosine;
    \textbf{5-mC:} 5-methylcytosine;
    \textbf{5-hmC:} 5-hydroxymethylcytosine;
    \textbf{5-fC:} 5-formylcytosine;
    \textbf{5-caC:} 5-carboxylcytosine;
    \textbf{TET:} Ten-eleven translocation;
    \textbf{APOBEC3A:} Apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like 3A

\section*{Availability of data and materials}
Sequencing data is from human tissue and will be available through dbGaP. Software used in processing the data is available on GitHub at \url{https://github.com/jamorrison/wgms_kit_comparison/releases/tag/v1.0.3}.

\section*{Ethics approval and consent to participate}
This study received ethics approval from the Van Andel Research Institute Institutional Review Board (IRB \#19017).

\section*{Competing interests}
The authors declare that they have no competing interests.

\section*{Consent for publication}
Not applicable.

\section*{Authors' contributions}
MA and HS conceptualized the study. JM led the study and performed the analyses. JM, MA and HS wrote the manuscript. DWC, LLR, and EJS provided the samples. KKF, JMK, and MA prepared samples and performed sequencing. BKJ and WZ provided intellectual input during analysis. All authors read and approved the manuscript.

\section*{Authors' information}
Not applicable.

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{vancouver}
\bibliography{wgms_kit_comparison_v1.3}

% Figures
\section*{Figures}
\begin{figure}[!htb] % Figure 1
	\caption{\textbf{Raw read statistics for each protocol.} \textbf{(A)} Percentage of bases with base quality $<$ 20 for read 1. \textbf{(B)} Percentage of bases with base quality $<$ 20 for read 2.} \label{fig:raw_trim}
\end{figure}

\begin{figure}[!htb] % Figure 2
	\caption{\textbf{Library quality metrics for each protocol for Sample A.} \textbf{(A)} The percentage of optimally, sub-optimally, and not aligned read fragments for each protocol. Note, read fragments treat reads 1 and 2 as separate entities, as it is possible that one read in the pair is mapped, while the other is not. \textbf{(B)} Insert size distribution. \textbf{(C)} Duplicate rate for reads with MAPQ $\geq40$. \textbf{(D)} The library complexity, which is a function of the duplicate rate. Metrics for Sample B are shown in Additional file 1: Figure S2.} \label{fig:lib_quality}
\end{figure}

\begin{figure}[!htb] % Figure 3
	\caption{\textbf{Library uniformity as measured by coverage of various genomic element categories.} Ratio of observed coverage to expected coverage for \textbf{(A)} all CpGs and \textbf{(B)} CpG islands. \textbf{(C)} Percentage of all CpGs covered by at least one unique read with MAPQ $\geq40$. \textbf{(D)} Percentage of CpGs in CpG islands covered by at least one unique read with MAPQ $\geq40$. Note, for (C) and (D), all libraries were downsampled to be comparable to PBAT (150 million reads, or $\sim$4.8X coverage, per sample, see Methods for details); therefore, any differences are not likely confounded by sequencing depth. As expected, this coverage will be substantially higher at increased depth.} \label{fig:cpg_coverage}
\end{figure}

\begin{figure} % Figure 4
	\caption{\textbf{The EPCAM promoter region as a representative example for data generated with the protocols.} The aligned reads tracks are taken from the Integrated Genomics Viewer (IGV) \cite{igv_2011} in the bisulfite mode mode where red represents an unconverted cytosine and blue represents a converted cytosine. Each panel represents one sample, with A and B denoting the biological replicates and 1 and 2 the technical replicates for each library construction protocol. The shown region is 1500 bp upstream and downstream of exon 1. The location of a CpG island is indicated with a green box on the bottom. Note, the strands for the PBAT samples have been flipped in silico before being displayed to account for the strand definition in the Miura and Ito protocol. The strands in the PBAT protocol are opposite from what is expected by IGV, as well as the definition used by the other protocols.} \label{fig:epcam_coverage}
\end{figure}

\begin{figure}[!htb] % Figure 5
	\caption{\textbf{Read-averaged cytosine retention by dinucleotide context, namely (A) CpA, (B) CpC, (C) CpG, and (D) CpT.} In each panel two technical replicates are shown for each biological replicate. The x-axis denotes percent retention, with a scale of 0-5\% for CpH panels and 0-100\% for the CpG panel.} \label{fig:read_avg_cret}
\end{figure}

\begin{figure}[!htb] % Figure 6
	\caption{\textbf{The NEB protocol has the highest correlation of beta values between Sample A technical replicates.} The Spearman correlation coefficient, $r_s$, between the two replicates is listed in each figure, along with the number of 100kb bins used in calculating the coefficient. Note, all libraries were downsampled to be comparable to PBAT; therefore, any differences are not likely confounded by sequencing depth. Overall low correlation values are due to low coverage from downsampling.} \label{fig:tech_corr}
\end{figure}

% Tables
\section*{Tables}
\begin{table}[!htb] % Table 1
	\caption{\textbf{Library preparation protocols used in analysis.} Throughout the text, the ``Short Name'' entry is used to describe which protocol is being discussed, rather than the full name. Sample and technical replicates names are included in parentheses.} \label{tab:data_summary}
	\begin{center}
		\resizebox{\textwidth}{!}{\begin{tabular}{l | c c c c c }
		\toprule
		Short Name & Protocol & Strand-specific? & DNA Input [ng] & \# Samples & \# Tech. Reps. \\
		\midrule
		Kapa & KAPA Hyper Prep & Yes & 300 & 2 (A/B) & 2 (1/2) \\
		NEB & NEBNext Enzymatic Methyl-seq & Yes & 200 & 2 (A/B) & 2 (1/2) \\
		PBAT & Miura F., Ito T. (2018) & No & 100 & 2 (A/B) & 1 (1) \\
		Swift & Accel-NGS Methyl-Seq DNA & Yes & 100 & 2 (A/B) & 2 (1/2) \\
		\midrule
		Low NEB & NEBNext Enzymatic Methyl-seq & Yes & 10 & 2 (A/B) & 2 (1/2) \\
		Low Swift & Accel-NGS Methyl-Seq DNA & Yes & 10 & 2 (A/B) & 2 (1/2) \\
		\bottomrule
		\end{tabular}}
	\end{center}
\end{table}

\end{backmatter}
\end{document}
