paper.lyx

#LyX 2.3 created this file. For more info see http://www.lyx.org/
\lyxformat 544
\begin_document
\begin_header
\save_transient_properties true
\origin unavailable
\textclass article
\use_default_options true
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman "default" "default"
\font_sans "default" "default"
\font_typewriter "default" "default"
\font_math "auto" "auto"
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100 100
\font_tt_scale 100 100
\use_microtype false
\use_dash_ligatures true
\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\float_placement tbph
\paperfontsize default
\spacing single
\use_hyperref true
\pdf_bookmarks true
\pdf_bookmarksnumbered false
\pdf_bookmarksopen false
\pdf_bookmarksopenlevel 1
\pdf_breaklinks false
\pdf_pdfborder false
\pdf_colorlinks false
\pdf_backref false
\pdf_pdfusetitle true
\papersize default
\use_geometry true
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine biblatex
\cite_engine_type numerical
\biblio_style plain
\biblatex_bibstyle numeric
\biblatex_citestyle numeric
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\use_minted 0
\index Index
\shortcut idx
\color #008000
\end_index
\leftmargin 2.5cm
\rightmargin 2.5cm
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\is_math_indent 0
\math_numbering_side default
\quotes_style english
\dynamic_quotes 0
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header

\begin_body

\begin_layout Title
Correcting Methylation Calls in Clinically Relevant Low-Mappability Regions
\end_layout

\begin_layout Author
Caiden M.
 Kumar
\begin_inset script superscript

\begin_layout Plain Layout
\begin_inset CommandInset ref
LatexCommand ref
reference "enu:New England Biolabs"
plural "false"
caps "false"
noprefix "false"

\end_inset


\end_layout

\end_inset

, Devon P.
 Ryan
\begin_inset script superscript

\begin_layout Plain Layout
\begin_inset CommandInset ref
LatexCommand ref
reference "enu:Genedata-AG"
plural "false"
caps "false"
noprefix "false"

\end_inset


\end_layout

\end_inset

, Bradley W.
 Langhorst
\begin_inset script superscript

\begin_layout Plain Layout
\begin_inset CommandInset ref
LatexCommand ref
reference "enu:New England Biolabs"
plural "false"
caps "false"
noprefix "false"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Enumerate
\begin_inset CommandInset label
LatexCommand label
name "enu:New England Biolabs"

\end_inset

New England Biolabs
\end_layout

\begin_layout Enumerate
\begin_inset CommandInset label
LatexCommand label
name "enu:Genedata-AG"

\end_inset

GeneData AG
\end_layout

\begin_layout Abstract
DNA methylation is an important component in vital biological functions
 such as embryonic development, carcinogenesis, and heritable regulation.
 Accurate methods to assess genomic methylation status are crucial to its
 effective use in many scenarios, especially in the detection and diagnosis
 of disease.
 Methylation aligners, such as Bismark and bwa-meth, frequently assign significa
ntly higher MapQ values than can be supported by the uniqueness of the region
 reads are mapped to.
 These incorrectly high MapQs result in inappropriate methylation calling
 in repetitive regions.
 We observe reads that should map to separate locations (possibly having
 different methylation states) actually end up mapping to the same locus,
 causing apparent mixed methylation at such loci.
 Methylation calling can be improved by using Bismap mappability data to
 filter out insufficiently unique reads.
 However, simply filtering out Cs in insufficiently unique regions is not
 adequate as it is prone to over-filtering Cs in small mappability dips.
 These Cs can in fact often be called using reads anchored in a nearby mappable
 region.
 We have created a new feature for the MethylDackel methylation caller to
 perform read-based filtering.
 This new methylation calling method resolves some of the apparent mixed
 methylation to either 0% or 100% methylation and removes many unsupportable
 methylation calls.
 We examined methylation calls with and without read-based filtering in
 or near the 7830 genes containing ClinVar variants in a methylation sequencing
 data set from the NA12878 cell line.
 Use of this improved method corrected 41,143 mixed methylation Cs to 0%
 methylation, and 22,345 to 100% methylation throughout the genome.
\end_layout

\begin_layout Section*
Introduction
\end_layout

\begin_layout Standard
As DNA methylation status can have a significant biological function 
\begin_inset CommandInset citation
LatexCommand cite
key "Schubeler2015"
literal "false"

\end_inset

, it is important that there be an accurate way of calling methylation on
 a genome.
 Although there are multiple varieties of DNA methylation, perhaps the most
 significant type in eukaryotes is methylation of cytosine to 5-methylcytosine
 
\begin_inset CommandInset citation
LatexCommand cite
key "breiling_epigenetic_2015"
literal "false"

\end_inset

.
 Data on DNA cytosine methylation state can be gathered using a methylation
 sequencing technique (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:seq"
plural "false"
caps "false"
noprefix "false"

\end_inset

), like bisulfite sequencing 
\begin_inset CommandInset citation
LatexCommand cite
key "Chen2017"
literal "false"

\end_inset

.
 In bisulfite sequencing, unmethylated cytosines are deaminated to uracil
 by the addition of sodium bisulfite.
 5-methylcytosines are not affected.
 Since uracil sequences as thymine and 5-methylcytosine sequences as cytosine,
 positions of unmethylated Cs in a reference sequence can be identified
 by C->T transitions
\begin_inset CommandInset citation
LatexCommand cite
key "Chen2017"
literal "false"

\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename fig.svg
	lyxscale 50
	width 100col%

\end_inset


\begin_inset Caption Standard

\begin_layout Plain Layout
Overview of Methylation Sequencing and Data Analysis Methods.
 
\begin_inset CommandInset label
LatexCommand label
name "fig:seq"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
It is also possible to use an enzymatic method, EM-seq, employing TET2 to
 oxidize 5-methylcytosine and an APOBEC enzyme to deaminate unmodified cytosines
 to uracil.
 While sodium bisulfite treatment produces double strand breaks and other
 DNA damage, this enzymatic method deaminates with more precision 
\begin_inset CommandInset citation
LatexCommand cite
key "vaisvila_enzymatic_2021"
literal "false"

\end_inset

 and longer, better mappable fragments.
\end_layout

\begin_layout Standard
Whichever method is used to deaminate cytosines, methylation sequence data
 is typically aligned to a reference genome using a methylation-aware aligner
 
\begin_inset CommandInset citation
LatexCommand cite
key "Garrett-Bakelman2015"
literal "false"

\end_inset

, which is specifically designed to handle the C->T transitions in methylation
 sequencing data when aligning the reads to a reference.
 Once aligned, a methylation caller is needed to analyze alignments and
 report the methylation status of particular cytosines (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:seq"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 The resulting data shows the methylation status of each callable cytosine
 in the genome and can therefore be used to find and study biologically
 significant DNA methylation sites.
\end_layout

\begin_layout Standard
There are two commonly used methylation callers, MethylDackel 
\begin_inset CommandInset citation
LatexCommand cite
key "ryan_methyldackel_nodate"
literal "false"

\end_inset

 and bismark_methylation_extractor 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btr167"
literal "false"

\end_inset

.
 The Bismark aligner creates XM, XR, and XG tags that are used by the bismark_me
thylation_extractor to call methylation.
 Bismark_methylation_extractor requires these larger bam files and is not
 compatible with bam files produced by other aligners.
 MethylDackel extract is not tightly coupled to a specific aligner or tag
 scheme, performing the filtering, quality assessment, and reference comparisons
 needed for methylation calling internally.
 MethylDackel efficiently loads a complete or targeted reference genome,
 examining alignments in chunks that are processed in parallel with all
 the context needed to differentiate methylation from C->T mutation.
 Flexible and optional thresholds are available for base-quality, positional
 trimming, adjacent base context (e.g.
 CpG or CHH), alignment flags, and MapQ.
\end_layout

\begin_layout Standard
A read must be unambiguously placed if it is to provide information about
 a specific region of a genome.
 Reads that equally match more than one area of the reference genome should
 not be used to assess methylation of any given C.
 To avoid calling Cs using reads derived from multiple genomic loci, methylation
 aligners (and read aligners in general) assign a MapQ value to each read
 alignment (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:seq"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 According to the SAM specification 
\begin_inset CommandInset citation
LatexCommand cite
key "the_sam/bam_format_specification_working_group_sequence_2018"
literal "false"

\end_inset

, MapQ is defined as: 
\begin_inset Quotes eld
\end_inset


\begin_inset Formula $-10\log_{10}\Pr\{\mathrm{mapping\ position\ is\ wrong}\}$
\end_inset

, rounded to the nearest integer
\begin_inset Quotes erd
\end_inset

.
 MapQ indicates how uniquely placed a read alignment is, that is, in how
 many other places could the read align to the reference.
 A low MapQ means that the read may align in many places throughout the
 genome (for instance, a read of centromeric satellite DNA would likely
 have a very low MapQ).
 A high MapQ indicates that the read likely aligns where it is placed and
 nowhere else in the genome.
\end_layout

\begin_layout Standard
A methylation caller can use accurate MapQ values to filter out reads with
 multiple placements in the genome, allowing the resulting methylation calls
 to accurately reflect their specific loci.
 Unfortunately, current aligners frequently produce implausible MapQ scores.
\end_layout

\begin_layout Standard
In this paper, we describe the problem that low-mappability creates for
 methylation calling, the extent of miscalled methylation in real datasets,
 and our implementation of an improved method that measurably reduces incorrect
 mixed-methylation calls improving differential methylation assessment,
 the use of methylation as a biomarker for disease progression, and understandin
g of development processes.
\end_layout

\begin_layout Section*
Results
\end_layout

\begin_layout Standard
While evaluating the methylation aligner bwa-meth 
\begin_inset CommandInset citation
LatexCommand cite
key "2014arXiv1401.1129P"
literal "false"

\end_inset

, we observed a significant number of reads with unexpectedly high MapQ
 values in repetitive regions (e.g.
 centromeres, Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:chrom-overview"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 Another commonly used aligner Bismark
\begin_inset CommandInset citation
LatexCommand cite
key "krueger_bismark:_2011"
literal "false"

\end_inset

 produced even more calls in repetitive regions (data not shown).
 After observing high MapQ reads in the centromere and larger repetitive
 regions, we investigated to see if smaller regions might also be too repetitive
 to support the high aligner MapQ estimates observed.
 We identified repetitive regions using data from Bismap 
\begin_inset CommandInset citation
LatexCommand cite
key "karimzadeh2018"
literal "false"

\end_inset

, a tool that counts the number of occurrences of every single K-mer of
 a particular length (in this case, k=100) in the genome to create a mappability
 score (ranging from 0 to 1 in increments of 0.01) for every base in the
 GRCh38 reference.
 Each individual K-mer is considered either mappable or not, and the overall
 mappability score is derived (indirectly) from how many different K-mers
 are mappable and how many are not.
 Bismap takes the effect of C->T conversion into account and therefore produces
 data which is applicable in the context of methylation sequencing 
\begin_inset CommandInset citation
LatexCommand cite
key "karimzadeh2018"
literal "false"

\end_inset

.
 Reads entirely contained within a region of low mappability should not
 have high MapQ values due to their repetitiveness, however we observed
 many high-MapQ reads in regions with very low or zero Bismap mappability
 (e.g.
 Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:example"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename chroms_ideogram.png
	lyxscale 30
	width 90col%

\end_inset


\begin_inset Caption Standard

\begin_layout Plain Layout
Cs in low mappability regions (blue), as well as those filtered out by per-read
 filtering that still passed MapQ filtering (red), are represented across
 all chromosomes.
 Per-read filtering removes unsupported calls while retaining millions of
 Cs that a naïve per-C filtering approach would needlessly eliminate.
 In this figure, the horizontal banded rectangles represent the chromosomes,
 the blue track above them shows all calls in low-mappability regions, and
 the red track corresponds to miscalls accurately removed by per-read filtering.
\begin_inset CommandInset label
LatexCommand label
name "fig:chrom-overview"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename fig1_svg.svg
	lyxscale 50
	width 90col%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
MapQ filtering does not remove all poorly mapped reads.
 This is an example of a coverage spike in a low-mappability region: clearly
 not sufficiently mappable, but still called with MapQ filtering.
 Mappability is only high in a small portion of this view (top track), but
 using MapQ alone alone enable calls of many Cs (MapQ Only track).
 the Per-read Kept track shows that only 2 sites have enough mappability
 to call and successfully eliminates the remaining miscalled Cs (Per-read
 ∆ track).
 Per C filtering eliminates all Cs in this region including the 2 that should
 be callable (Per-C Kept track).
 
\begin_inset CommandInset label
LatexCommand label
name "fig:example"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset

 
\end_layout

\begin_layout Standard
While low MapQ can indicate repetitiveness, even stringent MapQ thresholds
 cannot reliably select reads for safe methylation calling in regions containing
 repetitive DNA (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:overview-1"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 We considered excluding methylation calls on Cs found in low-mappability
 regions (e.g.
 using bedtools), but rejected this approach because it is prone to both
 over- and under-filtering.
 Cs in short, unique regions would be kept (under-filtered) even if the
 surrounding DNA is repetitive (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:discarded"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 However, this situation was not commonly observed.
 Over-filtering of Cs in short repetitive regions is a larger problem though.
 In this scenario, a C located in a small dip in mappability would be eliminated
 (over-filtered) despite coverage from read pairs anchored in nearby unique
 regions (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:kept"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 When calling using EM-seq reads aligned with bwa-meth to GRCh38, under-filterin
g is rare (between 
\begin_inset Formula $9.4*10^{-5}\%$
\end_inset

 and 
\begin_inset Formula $1.2*10^{-4}\%$
\end_inset

 of called Cs) but over-filtering is much more common (between 
\begin_inset Formula $1.2\%$
\end_inset

 and 
\begin_inset Formula $1.4\%$
\end_inset

 of called Cs) (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:overview-1"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 Although ~1% of called Cs might seem like a small fraction, it accounts
 for ~14.0-16.2 million Cs, while ~
\begin_inset Formula $10^{-4}\%$
\end_inset

 of Cs corresponds to ~1079-1388 Cs.
 In cases where reads from multiple low-mappability regions are placed onto
 one region (creating a coverage spike), filtering by coverage could possibly
 also remove the problematic region, but in cases where reads are simply
 mixed among multiple low-mappability regions (not creating a coverage spike),
 coverage filtering would not be effective at removing the affected regions.
 In addition, if some of the reads used to call a C are uniquely mappable,
 but others are not, per-read filtering can remove only those reads that
 are not uniquely mappable, rather than being forced to discard the entire
 locus as is required for per-C filtering.
\end_layout

\begin_layout Standard
To reliably filter out only problematic read pairs (those where both mates
 are placed in low mappability regions) we modified the MethylDackel methylation
 caller 
\begin_inset CommandInset ref
LatexCommand ref
reference "methyldackel-patch"
plural "false"
caps "false"
noprefix "false"

\end_inset

to accept a bigWig file of low-mappability regions to exclude from analysis.
 This per-read filtering approach precisely eliminates only those reads
 in repetitive regions and has a negligible effect on execution speed (see
 computational methods for details).
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename fig4_svg.svg
	lyxscale 70
	width 100col%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
An overview of the impact of the various filtering scenarios in the q arm
 of GRCh38 chromosome 1.
 The six tracks show the effect of per-C and per-read mappability filtering
 on the methylation calls on this chromosome arm and are representative
 of rest of the genome.
 
\begin_inset Quotes eld
\end_inset

Mappability
\begin_inset Quotes erd
\end_inset

 is the Bismap mappability for the chromosome.
 
\begin_inset Quotes eld
\end_inset

MapQ filtering only
\begin_inset Quotes erd
\end_inset

 is the calls after commonly used MapQ filtering only (no per-C or per-read
 mappability filtering).
 
\begin_inset Quotes eld
\end_inset

Kept by Per-read Filtering
\begin_inset Quotes erd
\end_inset

 shows the methylation calls kept after per-read mappability filtering.
 
\begin_inset Quotes eld
\end_inset

Removed by Per-read Filtering
\begin_inset Quotes erd
\end_inset

 shows those calls removed entirely by per-read filtering for being insufficient
ly mappable.
 
\begin_inset Quotes eld
\end_inset

Under-filtered Cs
\begin_inset Quotes erd
\end_inset

 are not sufficiently mappable according to per-read filtering, but are
 retained by per-C filtering.
 
\begin_inset Quotes eld
\end_inset

Over-filtered Cs
\begin_inset Quotes erd
\end_inset

 are sufficiently mappable according to per-read filtering, but are removed
 by per-C filtering.
\begin_inset CommandInset label
LatexCommand label
name "fig:overview-1"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename fig3_svg.svg
	lyxscale 50
	width 90col%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
An example of a C incorrectly discarded by per-C filtering (red box).
 The C is in a small mappability dip, causing it to be discarded by bedtools
 filtering, while per-read filtering correctly retains the C, as it can
 be called using reads extending into the high-mappability regions on either
 side.
\begin_inset CommandInset label
LatexCommand label
name "fig:discarded"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename fig2_svg.svg
	lyxscale 50
	width 90col%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
An example of a C incorrectly kept by per-C filtering (red box).
 As the C is itself just inside the edge of a mappable region, bedtools
 retains it, even though it is only supported by two insufficiently mappable
 reads.
 Our per-read filtering correctly removes this sort of C, as there is not
 enough data here to call it confidently.
\begin_inset CommandInset label
LatexCommand label
name "fig:kept"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
To focus on the incorrect alignments and methylation calls that have biological
 and medical significance, we examined methylation calls in and just upstream
 of GENCODE genes 
\begin_inset CommandInset citation
LatexCommand cite
key "Harrow2012"
literal "false"

\end_inset

 that contain variants listed in the ClinVar 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/nar/gkx1153"
literal "false"

\end_inset

 database of disease-associated variants.
 Using our alignment filtering approach we successfully avoided miscalling
 264 thousand Cs in these important regions in a 200 ng EM-seq sample which
 have insufficient mappability to support methylation calling.
 Briefly, to identify the number of methylation calls correctly removed
 (miscalls), we took the list of unfiltered Cs for the regions we were intereste
d in and removed from it all Cs we kept after per-read filtering.
 The result was a list of every C that was correctly filtered out (Figure
 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:Miscall-counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
\end_layout

\begin_layout Subsection*
Filtering Effect
\end_layout

\begin_layout Standard
We compared the effect of per-read mappability filtering on number of C's
 called, as well as number of C's called in low-mappability regions (as
 defined by per-C filtering), on EM-seq libraries of 3 sample masses aligned
 with the bwa-meth aligner.
 We saw a significant positive filtering effect on all 3 sample input amounts,
 with a mean of 7.67 million Cs filtered out by read-based filtering in each
 sample (see Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:Miscall-counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 The miscalled Cs were scattered across the genome, as shown in Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:chrom-overview"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
\end_layout

\begin_layout Standard
We also examined how many mixed methylation (0% < %Methyl < 100%) Cs were
 resolved to either 0% or 100% methylation, and found that approximately
 60,000-90,000 methylation calls in the EM-seq data that were previously
 mixed were resolved by per-read filtering to either fully unmethylated
 or fully methylated (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:genome-resolved-calls"
plural "false"
caps "false"
noprefix "false"

\end_inset

), showing that unlike per-C filtering, per-read filtering can correct methylati
on calls without always having to entirely delete them.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename Miscalls Overview.png
	lyxscale 45
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:Miscall-counts"

\end_inset

Left: Percents of miscalled Cs in the 10, 50, and 200 ng EM-seq and WGBS
 samples, aligned with bwa-meth and Bismark, where a 
\begin_inset Quotes eld
\end_inset

miscall
\begin_inset Quotes erd
\end_inset

 is a methylation call filtered out entirely by per-read filtering.
 The vertical axis is the percent of miscalled Cs in the sample.
 Right: Percents of miscalls within ClinVar genes in the same samples.
 The vertical axis is the percent of ClinVar miscalled Cs in the sample.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename Resolved Calls.png
	lyxscale 40
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:genome-resolved-calls"

\end_inset

Counts of genome-wide resolved calls in the 10, 50, and 200 ng EM-seq and
 WGBS samples, aligned with bwa-meth and Bismark, where a 
\begin_inset Quotes eld
\end_inset

resolved call
\begin_inset Quotes erd
\end_inset

 is a methylation call that is mixed in the unfiltered data, but either
 0% or 100% in the filtered data.
 The vertical axis is the number of resolved calls in the sample.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsection*
Effect on Medically Relevant Genes
\end_layout

\begin_layout Standard
To estimate biological significance of these miscalls, we used genes in
 GENCODE containing variants in ClinVar (
\begin_inset Quotes eld
\end_inset

ClinVar genes
\begin_inset Quotes erd
\end_inset

).
 Focusing on genes that may contain miscalls, we filtered the list to only
 include those genes which had at least one region, either in the gene or
 in a 2 kb region upstream of the annotated start site (to capture regulatory
 regions), where the mappability is too low to allow reliable methylation
 calling.
 Examination of this data in IGV
\begin_inset CommandInset citation
LatexCommand cite
key "robinson_integrative_2011"
literal "false"

\end_inset

 and numerically revealed that there were miscalls in or upstream of thousands
 of unique genes in the EM-seq/bwa-meth data, with a median number of 64
 miscalls per gene, indicating that indeed miscalls may affect biologically
 significant methylation.
 Many of these miscalls can be resolved to 0 or 100% methylation by the
 per-read filtering method described in this paper (Figures 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:genome-resolved-calls"
plural "false"
caps "false"
noprefix "false"

\end_inset

, 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:resolved_gene_counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename Miscalled Genes.png
	lyxscale 40
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:resolved_gene_counts"

\end_inset

Counts of ClinVar genes containing resolved (from intermediate to 0 or 100%
 methylation) miscalls in the 10, 50, and 200 ng EM-seq and WGBS samples,
 aligned with bwa-meth and Bismark.
 The vertical axis is the number genes containing miscalls that were resolved
 by mappability-informed methylation calling.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
It is evident that the sites of most methylation miscalls fall into two
 broad types.
 First, are repetitive regions that simply have a mappability of 0 for a
 long stretch of sequence.
 All Cs in such regions are classified as miscalls due to the lack of mappable
 regions to anchor reads with.
 Second, are shorter regions at the border between high mappability and
 low mappability regions.
 These regions are where per-C filtering produces over- and under-filtered
 Cs.
 Over-filtered Cs occur when a C is in the low mappability region, but there
 is a mappable read anchored in the high mappability region that allows
 it to be called.
 Under-filtered C's occur when a C is just inside a mappable region, but
 there is no sufficiently mappable read to call it.
 Examples of both from a 200 ng EM-seq sample aligned to GRCh38 using bwa-meth
 are provided in Figures 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:igv-1"
plural "false"
caps "false"
noprefix "false"

\end_inset

 and 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:igv-2"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename chr2_110091420-110092421.svg
	width 90col%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:igv-1"

\end_inset

A view of a cluster of miscalls in a large low-mappability region containing
 the MALL gene (ENSG00000144063.3) 
\begin_inset CommandInset citation
LatexCommand cite
key "perez_structural_1997"
literal "false"

\end_inset

, variants of which are associated with human disease
\begin_inset CommandInset citation
LatexCommand cite
key "fragoza_extensive_2019"
literal "false"

\end_inset

,
\begin_inset CommandInset citation
LatexCommand cite
key "kang_atypical_2015"
literal "false"

\end_inset

.
 The reads are shown paired, with a line connecting reads 1 and 2 in a pair,
 with CpG methylation sites highlighted.
 Miscalls are shown in the track above the coverage graph.
 Note that every CpG call (as well as calls in other contexts) are miscalled
 in this region, as there are no high-mappability regions to anchor reads
 with.
 Two of these Cs show mixed evidence of methylation state, which may be
 due to DNA originating from multiple loci all mapping in this location
 with falsely high MapQ.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename chr1_148819799-148820800.svg
	width 90col%

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:igv-2"

\end_inset

A view of a cluster of miscalls near the boundary between high and low mappabili
ty regions in the gene PDE4DIP (ENSG00000178104.19), regulation of which
 is reported to relate to dementia
\begin_inset CommandInset citation
LatexCommand cite
key "potashkin_bioinformatic_2020"
literal "false"

\end_inset

, schizophrenia
\begin_inset CommandInset citation
LatexCommand cite
key "camargo_disrupted_2007"
literal "false"

\end_inset

, pineoblastoma
\begin_inset CommandInset citation
LatexCommand cite
key "snuderl_recurrent_2018"
literal "false"

\end_inset

 and other human disease
\begin_inset CommandInset citation
LatexCommand cite
key "overlack_direct_2011,auer_rare_2015"
literal "false"

\end_inset

.
 The read alignments are shown paired, with a line connecting reads 1 and
 2, and highlighted for CpG methylation.
 Note that some methylation calls in the low-mappability region are not
 miscalls, as they can be called using read pairs extending into a high-mappabil
ity region.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsection*
EM-seq vs.
 WGBS
\end_layout

\begin_layout Standard
Due to the reduced DNA damage and longer possible read lengths that are
 found in in EM-seq libraries, we compared with WGBS sequence libraries
 produced using the same 3 sample masses, to see if there is a difference
 between the two methods with regard to the number of miscalls.
 WGBS libraries have fewer overall Cs called
\begin_inset CommandInset citation
LatexCommand cite
key "vaisvila_enzymatic_2021"
literal "false"

\end_inset

 and therefor fewer miscalled Cs at all input masses than EM-seq libraries.
 In ClinVar gene regions, where both methods provide more similar overall
 coverage, the percent of miscalls is also similar between EM-seq and to
 WGBS libraries ( Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:Miscall-counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 We expected that longer reads could allow for more low-mappability Cs to
 be called using reads and mates extending into higher-mappability regions.
 Importantly, we observed that EM-seq libraries are able to resolve many
 more miscalled Cs throughout the genome (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:genome-resolved-calls"
plural "false"
caps "false"
noprefix "false"

\end_inset

) and in ClinVar regions than WGBS libraries 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:resolved_gene_counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
\end_layout

\begin_layout Subsection*
Bwa-meth vs.
 Bismark
\end_layout

\begin_layout Standard
For comparison, we examined the Bismark methylation aligner 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btr167"
literal "false"

\end_inset

 in the same manner as the bwa-meth analysis.
 We observed that compared to bwa-meth, Bismark produced more inappropriately-hi
gh MapQ reads in repetitive regions (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:chrom-overview"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 In addition, Bismark and bwa-meth use different systems to define values
 for MapQ.
 Bismark (as a result of using Bowtie2 for alignment) reports MapQ based
 on number of reference mismatches, producing values between 0 and 42 
\begin_inset CommandInset citation
LatexCommand cite
key "noauthor_qc_nodate,johnurbangenome_biofinysics_2014"
literal "false"

\end_inset

.
 Bwa-meth (using BWA-MEM for alignments) follows the SAM specification in
 estimating probabilities in the MapQ field.
 Overall, we observed a higher percent of miscalls with Bismark than we
 did with bwa-meth, as shown in Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:Miscall-counts"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
 This difference was noticeable in both overall miscalls and miscalls in
 ClinVar genes.
\end_layout

\begin_layout Subsection*
GRCh38 vs.
 T2T
\end_layout

\begin_layout Standard
We also examined the effect of using the T2T reference 
\begin_inset CommandInset citation
LatexCommand cite
key "nurk_complete_2021"
literal "false"

\end_inset

 in place of GRCh38, since the T2T reference is much more complete than
 GRCh38, including repetitive regions such as centromeres, telomeres, and
 rDNA, which could play a role in where methylation is called and what areas
 of the genome have unique enough sequence to support calling.
 We found that, as we expected, there were far fewer miscalls in the T2T
 data as compared to the GRCh38 data, with less than half the number of
 miscalls in T2T as in GRCh38.
 Considering both strands, total calls were roughly similar between the
 two reference sequences with more Cs available in the T2T reference.
 ( Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:t2t"
plural "false"
caps "false"
noprefix "false"

\end_inset

, left).
 Use of per-read filtering has a smaller effect with the more complete reference
, but thousands of Cs are still affected (Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:t2t"
plural "false"
caps "false"
noprefix "false"

\end_inset

, right).
 This indicates that using the T2T reference can allow for better calling
 of methylation with fewer miscalls, however without a high quality annotation
 available, this may not yet be a practical option.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename T2T Overview.png
	lyxscale 40
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:t2t"

\end_inset

Left: Total call counts before filtering in the 10, 50, and 200 ng EM-seq
 samples, aligned with bwa-meth to GRCh38 and T2T.
 The vertical axis is the number of calls in the sample.
 Right: Percents of miscalled Cs in the same samples, where a 
\begin_inset Quotes eld
\end_inset

miscall
\begin_inset Quotes erd
\end_inset

 is a methylation call filtered out entirely by per-read filtering.
 The vertical axis is the percent of miscalled Cs in the sample.
 In both figures, the colors define which sequencing method and aligner
 were used, as detailed in the legend.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Section*
Materials and Methods
\end_layout

\begin_layout Subsection*
DNA Methylation Sequencing
\end_layout

\begin_layout Subsection*
Materials
\end_layout

\begin_layout Standard
Previously published
\begin_inset CommandInset citation
LatexCommand cite
key "vaisvila_enzymatic_2021"
literal "false"

\end_inset

 libraries from 10, 50, and 200 ng of genomic DNA from the NA12878 cell
 line (Coriel) were analyzed.
 This input DNA was supplemented with a small amount of lambda phage DNA
 (NEB #N3011), and a pUC19 plasmid (NEB #N3041) treated with M.SssI CpG Methyltra
nsferase (NEB #M0226).
\end_layout

\begin_layout Subsection*
Whole Genome Bisulfite Libraries
\end_layout

\begin_layout Standard
The whole genome bisulfite (WGBS) libraries were prepared using the Ultra
 II DNA library prep kit (NEB E7645) before being Bisulfite converted using
 the Zymo EZ DNA Methylation-Gold bisulfite conversion kit according to
 the protocol described in Vaisvila et al.
 
\begin_inset CommandInset citation
LatexCommand cite
key "vaisvila_enzymatic_2021"
literal "false"

\end_inset

.
\end_layout

\begin_layout Subsection*
EM-seq Libraries
\end_layout

\begin_layout Standard
The EM-seq libraries were prepared using the NEBNext Ultra II DNA library
 prep kit (NEB E7645) (using the NEBNext EM-seq adapter), then EM-seq converted
 using TET2 and ABOBEC3A, as described in Vaisvila et al.
 
\begin_inset CommandInset citation
LatexCommand cite
key "vaisvila_enzymatic_2021"
literal "false"

\end_inset

.
\end_layout

\begin_layout Subsection*
Sequencing
\end_layout

\begin_layout Standard
Methylation libraries were sequenced with more diverse libraries (~10%)
 on 2 flow cells of an Illumina NovaSeq 6000 
\begin_inset CommandInset citation
LatexCommand cite
key "novaseq-6000"
literal "false"

\end_inset

 using the S4 chemistry.
 We analyzed about 1.3 billion 99 bp paired-end reads for the for each library.
\end_layout

\begin_layout Standard
\align center
\begin_inset Float table
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Tabular
<lyxtabular version="3" rows="7" columns="3">
<features tabularvalignment="middle">
<column alignment="center" valignment="top">
<column alignment="center" valignment="top">
<column alignment="center" valignment="top">
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Sequencing Method
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Sample Mass
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Read Count
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
WGBS
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
10 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,268,442,060
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
WGBS
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
50 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,269,601,144
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
WGBS
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
200 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,269,737,220
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
EM-seq
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
10 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,298,351,192
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
EM-seq
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
50 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,296,831,526
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
EM-seq
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
200 ng
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
1,296,793,470
\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
Sequencing read counts after downsampling and adapter/quality trimming broken
 down by sequencing method and sample mass.
 Read counts are reported as counts of individual reads, not of pairs.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsection*
Computational Methods
\end_layout

\begin_layout Standard
The following data and tools were used for the analysis:
\end_layout

\begin_layout Standard
The GRCh38.p11 analysis set (hereafter referred to as 
\begin_inset Quotes eld
\end_inset

GRCh38
\begin_inset Quotes erd
\end_inset

) supplemented with phage T4, phage lambda, phage Xp12, and pUC19 contigs
 was used throughout 
\begin_inset CommandInset citation
LatexCommand cite
key "church_modernizing_2011"
literal "false"

\end_inset

.
 A VCF file of disease-associated variant sites was downloaded from ClinVar
 (see supplemental data 
\begin_inset CommandInset ref
LatexCommand ref
reference "The-ClinVar-VCF"
plural "false"
caps "false"
noprefix "false"

\end_inset

 for file date) and a GFF file of the GENCODE v31 gene annotations were
 used.
 The mappability data was the 100bp multi-read bigWig file downloaded from
 Bismap (see supplemental materials 
\begin_inset CommandInset ref
LatexCommand ref
reference "The-Bismap-file"
plural "false"
caps "false"
noprefix "false"

\end_inset

 for link).
 
\end_layout

\begin_layout Standard
\begin_inset CommandInset label
LatexCommand label
name "nextflow script"

\end_inset

The Nextflow scripts used to analyze this data can be found at 
\begin_inset Flex URL
status open

\begin_layout Plain Layout

https://github.com/nebiolabs/low_bismap_methyl_calls
\end_layout

\end_inset

.
 A detailed diagram of the analysis pipeline is found in Supplemental figures
 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:dag"
plural "false"
caps "false"
noprefix "false"

\end_inset

 and 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:dag-1"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
 Tools used include bedtools 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btq033"
literal "false"

\end_inset

, samtools 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btp352"
literal "false"

\end_inset

, GNU awk, MethylDackel, bigWigToBedGraph 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btq351"
literal "false"

\end_inset

, BEDOPS 
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/bts277"
literal "false"

\end_inset

, GNU sort, GNU head, bwa-meth, and Bismark.
 GNU sort and head are required since some of the functionality needed (specific
ally, the ability to parallelize sorting and the ability to use a negative
 value with the -n option for head to count lines from the end of the file)
 is not present in BSD sort and head.
\end_layout

\begin_layout Section*
Discussion
\end_layout

\begin_layout Standard
Reads placed with falsely high confidence elicit cascading detrimental effects
 on methylation calling, differential methylation assessment, and assessment
 of phenotypes associated with methylation status.
 To remove the spurious, unsupported methylation calls that result from
 such reads, we suggest that reads with both mates in low mappability regions
 (as determined by Bismap) should be excluded from methylation calling.
 Additionally, because of the more accurate MapQ values, decreased run time,
 and more flexibility to separate methylation calling from alignment, we
 also recommend the use of bwa-meth for alignment and MethylDackel with
 MapQ > 10 and the new mappability feature enabled for methylation calling.
\end_layout

\begin_layout Section*
\begin_inset CommandInset label
LatexCommand label
name "sec:Supplemental-Materials"

\end_inset

Supplemental Materials
\end_layout

\begin_layout Standard
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
setcounter{figure}{0} 
\backslash
renewcommand{
\backslash
thefigure}{S.
\backslash
arabic{figure}}
\end_layout

\end_inset


\end_layout

\begin_layout Subsection*
Detailed Computational Methods
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename analysis_clinvar.svg
	lyxscale 90
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
Detailed description of the procedure used to annotate the T2Tv1.1 reference
 using the GRCh38 ClinVar and GENCODE data.
 Hexagons are input files.
 Ovals are intermediate files.
 Gray rhombuses are processing steps with multiple arguments.
 Parallelograms are output files.
 Labels on edges are processing steps with one input and one output.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename procfig_1.svg
	lyxscale 65
	width 60page%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:dag"

\end_inset

Detailed description of preprocessing steps used to identify biologically
 relevant low-mappability regions.
 Ovals are intermediate files.
 Gray rhombuses are processing steps with multiple arguments.
 Dotted lines indicate inputs that go to the methylation calling steps in
 Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:dag-1"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename procfig_2.svg
	lyxscale 65
	width 70page%
	scaleBeforeRotation

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "fig:dag-1"

\end_inset

Detailed description of steps for calling methylation and intersecting the
 resulting calls with the low mappability regions determined in the preprocessin
g steps.
 Ovals are intermediate files.
 Gray rhombuses are processing steps with multiple arguments.
 Parallelograms are output files.
 Dotted lines indicate inputs that come from the preprocessing steps in
 Figure 
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:dag"
plural "false"
caps "false"
noprefix "false"

\end_inset

.
 Parallelograms are output files.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
The sequencing reads were aligned using the bwa-meth aligner with default
 options.
 Methylation calling was performed on the resulting BAMs using MethylDackel
 v0.6.0.
 MethylDackel uses a default value of 10 as a minimum MapQ which should
 include mostly single locus reads, however inaccurate MapQs can lead to
 reads being incorrectly included in methylation calling.
 In order to eliminate reads in low-mappability regions, a patch was created
 for MethylDackel which allows it to take as an input a bigWig file which
 is then used to filter out read pairs (this patch currently only supports
 paired-end reads) where neither mate intersects a high-mappability region.
 The patch allows for user configuration of the low mappability threshold
 and the number of bases which must be equal to or above that threshold
 in order for the read pair to be kept (the defaults are a low mappability
 threshold of 0.01 and to require 15 bases that are greater than or equal
 to that threshold in a single read).
 The filtering algorithm has been optimized both by loading the mappability
 data into memory before calling, and through the use of a custom run-length
 compressed binary file format that we here term BBM (Binary BisMap
\begin_inset CommandInset ref
LatexCommand ref
reference "The-BBM-format"
plural "false"
caps "false"
noprefix "false"

\end_inset

).
 This format can store the mappability data for the hg38/GRCh38 human genome
 in 143 MB, compared to 1.11 GB when stored as a bigWig, giving a compression
 ratio for this dataset of 7.78:1.
 MethylDackel's trimming options were set to trim off the first 1 bp of
 read 1 and the last 2 bp of read 2.
\end_layout

\begin_layout Standard
The ClinVar VCF and GENCODE GFF were combined using bedtools and BEDOPS
 into a single BED file listing all GENCODE annotation regions that overlap
 one or more ClinVar variants (these regions will be referred to as 
\begin_inset Quotes eld
\end_inset

ClinVar regions
\begin_inset Quotes erd
\end_inset

), which will be called the 
\begin_inset Quotes eld
\end_inset

ClinVar regions BED
\begin_inset Quotes erd
\end_inset

.
 Each region in the ClinVar regions BED was expanded 2 kb upstream using
 bedtools slop to capture methylation upstream of the gene.
 The Bismap bedGraph (which did not contain zeroes) and an FAI index of
 the GRCh38 reference genome FASTA file were processed with bedtools, awk,
 and GNU sort to obtain a bedGraph containing all mappability data, including
 zeroes.
 This file will be referred to as the 
\begin_inset Quotes eld
\end_inset

Bismap complete bedGraph
\begin_inset Quotes erd
\end_inset

.
 It was then combined with the ClinVar regions BED using bedtools map to
 create a file containing the minimum and mean mappability for each gene
 with a ClinVar variant.
\end_layout

\begin_layout Standard
The Bismap complete bedGraph was then filtered using awk to produce a file
 containing only low mappability regions (mappability < 0.01).
 The file with minimum and mean mappability for every ClinVar region was
 filtered likewise on minimum mappability.
 The two resulting files (one of low-mappability regions, one of ClinVar
 regions with low minimum mappability) were combined to produce a file of
 all low mappability regions that are in ClinVar regions.
\end_layout

\begin_layout Standard
Alignments of 2x101-bp paired-end EM-seq and whole genome bisulfite reads
 were processed using MethylDackel (with and without the custom patch) using
 a minimum MapQ cutoff of 10 and the default settings mentioned above and
 combined with the file of all low mappability regions in ClinVar regions
 to produce a list of all methylation calls in low mappability ClinVar regions
 (this will be referred to as the 
\begin_inset Quotes eld
\end_inset

low mappability ClinVar calls file
\begin_inset Quotes erd
\end_inset

).
 A low-mappability methylation call, as used here, is defined as a methylation
 call that is in a region with Bismap mappability less than 0.01, as defined
 by the position of the C alone (positions of anchoring reads are not considered
 when determining this).
\end_layout

\begin_layout Standard
The low mappability ClinVar calls file was intersected with the file of
 ClinVar regions with low minimum mappability to produce a file of ClinVar
 regions with low mappability calls.
 The -wa option for bedtools intersect was used here, which writes a copy
 of the ClinVar region to the output file for each low mappability call
 in the region, in order that this file would contain multiple copies of
 each region, one per low mappability call in the region.
 These duplicates were then used to count low mappability calls by feeding
 the data to a custom Python script (found in the Nextflow script in the
 supplemental materials 
\begin_inset CommandInset ref
LatexCommand ref
reference "nextflow script"
plural "false"
caps "false"
noprefix "false"

\end_inset

) which counted and combined the duplicates, producing a list of all ClinVar
 regions with low mappability calls and how many low mappability calls are
 in each region.
 A similar intersection was performed for all ClinVar regions with low minimum
 mappability (i.e.
 ClinVar regions that could potentially be affected) to produce a file of
 all calls (regardless of mappability of the C) in such regions, and a list
 of such regions with a count of calls per region.
\end_layout

\begin_layout Standard
Since this analysis was run for three different input masses and two sequencing
 protocols, the low mappability calls files were also processed through
 a custom Python script (see supplemental materials) to add a field specifying
 which input mass, sequencing protocol, and MethylDackel filtering setting
 were used.
 As the input mass would have been difficult to parse out of the BAM name
 due to inconsistent file name formatting, a CSV was created mapping BAM
 file names to input masses and given as input to this step.
 The same field was added to the lists of all ClinVar regions with low mappabili
ty calls and counts described previously.
 We also produced similarly-annotated files for all methylation calls (regardles
s of mappability), all methylation calls in GENCODE genes (with low minimum
 mappability) which contained ClinVar variants, all methylation calls in
 low-mappability regions that were also in GENCODE genes which contained
 ClinVar variants, counts of calls in in GENCODE genes which contained ClinVar
 variants regardless of mappability, methylation calls that are either 0%
 or 100% (regardless of mappability), and methylation calls which are neither
 0% nor 100% (also regardless of mappability).
\end_layout

\begin_layout Standard
To assess miscalls (i.e.
 which Cs are filtered out by per-read filtering), as well as under/over-filteri
ng with respect to per-C filtering (discussed above), and resolution of
 mixed methylation calls to 0% or 100%, pairs of these files were compared
 using bedtools to produce files showing the difference between filtered
 and unfiltered data.
 All files were examined and compared in Tableau
\begin_inset script superscript

\begin_layout Plain Layout
®
\end_layout

\end_inset

.
\end_layout

\begin_layout Standard
To compare the behavior of Bismark with bwa-meth, this analysis was re-run
 with the Bismark aligner using default settings and deduplicated using
 deduplicate_Bismark according to the documented protocol 
\begin_inset CommandInset citation
LatexCommand cite
key "krueger_felixkruegerbismark_2021"
literal "false"

\end_inset

.
 The name of the aligner was added to the field specifying the input mass,
 protocol, and MethylDackel filtering settings present in all output files
 containing methylation calls or counts of such calls in ClinVar regions
 (using the same custom Python script).
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename realign.svg
	width 70page%

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption Standard

\begin_layout Plain Layout
Detailed description of BAM realignment steps.
 Hexagons are input files.
 Ovals are intermediate files.
 Gray rhombuses are processing steps with multiple arguments.
 Parallelograms are output files.
 Labels on edges are processing steps with one input and one output.
\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
To compare the need for and the effect of mappability filtering on the more
 complete T2T v1.1 reference genome 
\begin_inset CommandInset citation
LatexCommand cite
key "nurk_complete_2021"
literal "false"

\end_inset

 (hereafter referred to as 
\begin_inset Quotes eld
\end_inset

T2T
\begin_inset Quotes erd
\end_inset

), the analysis was re-run using the T2T reference in place of GRCh38 (supplemen
ted with phage T4, phage lambda, phage Xp12, and pUC19 just as with GRCh38).
 To use the ClinVar and GENCODE datasets, which were respectively originally
 obtained as VCF and GFF files containing genomic coordinates in GRCh38,
 with the T2T reference, the positions were transferred over to the T2T
 reference using CrossMap 
\begin_inset CommandInset citation
LatexCommand cite
key "zhao_crossmap_2014"
literal "false"

\end_inset

.
 This lift-over was done in two steps with separate chain files.
 First, the GRCh38 regions were lifted over using CrossMap to T2Tv1.0 using
 a chain file published on GitHub by Nico Alavi 
\begin_inset CommandInset citation
LatexCommand cite
key "alavi_burgshrimpsliftover_t2t_2021"
literal "false"

\end_inset

.
 Then, the lifted-over T2Tv1.0 annotations were lifted over again to T2Tv1.1
 using a chain file from the T2T Consortium 
\begin_inset CommandInset citation
LatexCommand cite
key "noauthor_marblchm13_2021"
literal "false"

\end_inset

.
 Specifically, the ClinVar regions BED, with the ClinVar and GENCODE data
 already combined, was lifted over.
 This converted file was used in place of the ClinVar regions BED for analysis
 with the T2T reference.
 As there was no Bismap mappability file for the T2T reference, we downloaded
 the Bismap tool and generated a 
\begin_inset Formula $k=100$
\end_inset

 multi-read Bismap mappability dataset for this reference genome.In total,
 this analysis was run on all combinations of sequencing method (WGBS or
 EM-seq), aligner (Bismark or bwa-meth), input mass (10, 50, or 200 ng),
 filtering (read-based filtering or no read-based filtering), and reference
 genome (GRCh38 or T2T), using a Nextflow
\begin_inset CommandInset citation
LatexCommand cite
key "di_tommaso_nextflow_2017"
literal "false"

\end_inset

 pipeline (see supplemental materials 
\begin_inset CommandInset ref
LatexCommand ref
reference "nextflow script"
plural "false"
caps "false"
noprefix "false"

\end_inset

) which managed the execution of all the analysis tools needed.
\end_layout

\begin_layout Subsection*
External Analysis Resources 
\end_layout

\begin_layout Standard
Individual tools used in the analysis include sambamba
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btv098"
literal "false"

\end_inset

, bigWigToBedGraph
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btq351"
literal "false"

\end_inset

, MethylDackel, awk, python, GNU sort, GNU head, bedtools
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/btq033"
literal "false"

\end_inset

, and BEDOPS
\begin_inset CommandInset citation
LatexCommand cite
key "doi:10.1093/bioinformatics/bts277"
literal "false"

\end_inset

.
 Versions of all tools used are specified in the nextflow 21.04.0 scripts
 using conda
\begin_inset CommandInset citation
LatexCommand cite
key "noauthor_miniconda_2021"
literal "false"

\end_inset

 dependency resolution.
\end_layout

\begin_layout Standard
The custom Python script that adds the input mass, sequencing protocol,
 MethylDackel filtering setting, and aligner name can be found at 
\begin_inset Flex URL
status open

\begin_layout Plain Layout

https://github.com/nebiolabs/low_bismap_methyl_calls
\end_layout

\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset CommandInset label
LatexCommand label
name "methyldackel-patch"

\end_inset

The pull request for the patch adding mappability support to MethylDackel
 can be found at 
\begin_inset Flex URL
status collapsed

\begin_layout Plain Layout

https://github.com/dpryan79/MethylDackel/pull/80
\end_layout

\end_inset

.
 It was merged into MethylDackel in version 0.5.0.
\end_layout

\begin_layout Standard
\begin_inset CommandInset label
LatexCommand label
name "The-BBM-format"

\end_inset

 The BBM format definition is available at 
\begin_inset Flex URL
status open

\begin_layout Plain Layout

https://github.com/dpryan79/MethylDackel/blob/master/BBM_Specification.md
\end_layout

\end_inset

.
 
\end_layout

\begin_layout Standard
\begin_inset CommandInset label
LatexCommand label
name "The-Bismap-file"

\end_inset

The Bismap file used for this analysis was downloaded from 
\begin_inset Flex URL
status open

\begin_layout Plain Layout

https://www.pmgenomics.ca/hoffmanlab/proj/bismap/trackhub/hg38/k100.Bismap.MultiTrac
kMappability.bw
\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset CommandInset label
LatexCommand label
name "The-ClinVar-VCF"

\end_inset

The ClinVar VCF was the July 22, 2019 version of ClinVar's variants VCF,
 with a file name of ClinVar_20190722.vcf.gz
\end_layout

\begin_layout Standard
The chromosome ideogram figure was generated using karyoploteR
\begin_inset CommandInset citation
LatexCommand cite
key "gel_karyoploter_2017"
literal "false"

\end_inset

.
\end_layout

\begin_layout Subsection*
Acknowledgments
\end_layout

\begin_layout Standard
We are grateful to the many people who contributed to the EM-seq project
 and made these libraries available.
 This work would not have been possible with out access to an extensible,
 open-source methylation caller.
 Many thanks to Chaithanya Ponnaluri for critical reading of this manuscript,
 to Matt Campbell and Ariel Erijman for helpful discussion, and to the organizer
s and supporters of the summer internship program at New England Biolabs
 that enabled this project.
\end_layout

\begin_layout Standard
\begin_inset CommandInset bibtex
LatexCommand bibtex
btprint "btPrintCited"
bibfiles "4641904citation"
options "plain"

\end_inset


\end_layout

\end_body
\end_document