\documentclass[10pt]{article}
\usepackage[hmargin=1.5cm,top=2cm,bottom=2cm]{geometry}
\usepackage{multicol}
\setlength\columnsep{15pt}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{array}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage[auth-sc]{authblk}
\usepackage{longtable}
\usepackage{multirow}
\usepackage{hyperref}
\usepackage{enumerate}
\usepackage[labelfont=bf]{caption}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{mdframed}
\usepackage{graphics}
\usepackage{multirow}
\usepackage{rotating}
\usepackage{array}
\usepackage{capt-of}
\usepackage{lscape}
\usepackage{caption}
\usepackage{breakurl}
\usepackage{todonotes}
\usepackage{hanging}
\usepackage{pagecolor}
\usepackage[final]{pdfpages}
\usepackage[leftFloats,CaptionAfterwards]{fltpage}
\usepackage[numbers,super,sort&compress]{natbib}
\setlength{\bibsep}{0pt plus 0.3ex}
\usepackage{abstract}
\usepackage{enumitem}
\usepackage{soul}
\usepackage{titlesec}
\titleformat{\section}[block]{\large\bfseries\filcenter}{\thesection.}{0.4em}{}
\titleformat{\subsection}[block]{\normalsize\sc\bfseries\filcenter}{\thesubsection.}{0.4em}{}
\titleformat{\subsubsection}[block]{\normalsize\sc\itshape\filright}{\thesubsection.}{0.4em}{}
\setcounter{secnumdepth}{5}

\usepackage[hang]{footmisc}
\setlength\footnotemargin{0em}

\setlength{\skip\footins}{0.75cm}

\makeatletter
\renewcommand\footnoterule{%
  \kern-3\p@
  \hrule\@width \textwidth height 1.5pt
  \kern2.6\p@}
\makeatother

\makeatletter
\def\@biblabel#1{\@ifnotempty{#1}{#1.}}
\makeatother

\newcommand{\filllastline}[1]{
\setlength\leftskip{0pt}
\setlength\rightskip{0pt}
\setlength\parfillskip{0pt}
#1}

\newenvironment{Figure}
{\par\medskip\noindent\minipage{\linewidth}}
{\endminipage\par\medskip}

\title{\bf Long-range single-molecule mapping of chromatin accessibility in eukaryotes}
\renewcommand\Authfont{\scshape\normalsize}
\author[1,*]{Zohar Shipony}
\author[1,*]{Georgi K. Marinov}
\author[5]{Matthew P. Swaffer}
\author[1]{Nasa A. Sinott-Armstrong}
\author[5]{Jan M. Skotheim}
\author[1,2]{Anshul Kundaje}
\author[1,3,4,$\#$]{William J. Greenleaf}
\renewcommand\Affilfont{\itshape\small}
\affil[1]{Department of Genetics, Stanford University, Stanford, CA 94305, USA}
\affil[2]{Department of Computer Science, Stanford University, Stanford, CA 94305, USA}
% \affil[3]{Center for Personal Dynamic Regulomes, Stanford University, Stanford, CA 94305, USA}
\affil[3]{Department of Applied Physics, Stanford University, Stanford, CA 94305, USA}
\affil[4]{Chan Zuckerberg Biohub, San Francisco, California, USA}
\affil[5]{Department of Biology, Stanford University, Stanford, CA 94305, USA}
\affil[*]{These authors contributed equally to this work}
\affil[$\#$]{Corresponding author}
\date{}

\begin{document}
\maketitle

% \centerline{}
% \centerline{}
\begin{abstract}

% https://www.nature.com/nmeth/about/content

\noindent {\normalsize \textbf{Active regulatory elements in eukaryotes are typically characterized by an open, nucleosome-depleted chromatin structure; mapping areas of open chromatin has accordingly emerged as a widely used tool in the arsenal of modern functional genomics. However, existing approaches for profiling chromatin accessibility are limited by their reliance on DNA fragmentation and short read sequencing, which leaves them unable to provide information about the state of chromatin on larger scales or reveal coordination between the chromatin state of individual distal regulatory elements. To address these limitations, we have developed a method for profiling accessibility of individual chromatin fibers at multi-kilobase length scale (SMAC-seq, or \underline{\textbf{S}}ingle-\underline{\textbf{M}}olecule long-read \underline{\textbf{A}}ccessible \underline{\textbf{C}}hromatin mapping \underline{\textbf{seq}}uencing assay), enabling the simultaneous, high-resolution, single-molecule assessment of the chromatin state of distal genomic elements. Our strategy is based on combining the preferential methylation of open chromatin regions by DNA methyltransferases (CpG and GpC 5-methylcytosine (5mC) and N$^6$-methyladenosine (m$^6$A) enzymes) and the ability of long-read single-molecule nanopore sequencing to directly read out the methylation state of individual DNA bases. Applying SMAC-seq to the budding yeast \textit{Saccharomyces cerevisiae}, we demonstrate that aggregate SMAC-seq signals match bulk-level accessibility measurements, observe single-molecule protection footprints of nucleosomes and transcription factors, and quantify the correlation between the chromatin states of distal genomic elements.}
}
\centerline{}
% \centerline{}
\end{abstract}

\renewcommand{\footnotesize}{\fontsize{10pt}{12pt}\selectfont}

\begin{multicols}{2}

\filllastline{The packaging of DNA by nucleosomes into chromatin is a major organizing principle of genome organization in eukaryotes. The majority of the genome is tightly packaged by nucleosomal particles that wrap around DNA (usually $\sim$147bp), thus making it inaccessible to binding by most regulatory proteins. Conversely, regions of open chromatin tend to be strongly associated with regulatory elements (REs), such as enhancers, promoters, and insulators, and nucleosomes often exhibit characteristic occupancy patterns in their vicinity. These biological properties have proven highly useful for identifying candidate such elements (cREs), and in turn for understanding the functional organization of genomes. Regions of open chromatin have greatly increased sensitivity to cleavage by nucleases such as DNAse I, as already noted nearly four decades ago for promoter and enhancer elements around individual genes\cite{Wu1980,Keen1981,McGhee1981}. Subsequent advances in microarray\cite{Dorschner2004,Sabo2006} and DNA sequencing technologies\cite{Boyle2008,Hesselberth2009} have enabled DNAse hypersensitivity-based mapping of cREs genome-wide. Similarly, digestion of DNA is inhibited by nucleosome occupancy, and MNase digestion of chromatin (MNase-seq) has accordingly become a widely used tool to map nucleosome positioning throughout genomes\cite{Schones2008}. More recently, the Tn5 transposase has also been used as a facile probe of chromatin accessibility\cite{Buenrostro2013}. However, while short read-based assays provide immensely useful information about the identity of cREs and positioned nucleosomes, they give little insight into the long-range physical organization of individual chromatin fibers. Because cleavage-based approaches remove the linkage between distal segments of DNA molecules, classical functional genomic assays return estimates of the relative chromatin states only at short localized regions of the genome. To what extent these states are coordinated across distance and what the distribution of} 

\end{multicols}

\begin{figure*}
\begin{center}
\includegraphics[width=18.5cm]{Fig1V5.png}
% \includegraphics[angle=90,origin=c,width=17.5cm]{Fig1.png}
% \includegraphics[width=23.5cm]{Fig1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf The SMAC-seq assay for profiling chromatin accessibility and nucleosome positioning at the multikilobase scale}. 
(a) Outline of the SMAC-seq assay. Intact chromatin is treated with m$^6$A and CpG and GpC 5mC methyltransferases, which preferentially methylate DNA bases in open chromatin regions. HMW DNA is then isolated, subjected to nanopore sequencing and methylated bases are used to reconstruct the open chromatin state within individual molecules. 
(b-h) SMAC-seq faithfully captures chromatin accessibility around promoters and positioned nucleosomes in \textit{S. cerevisiae}; 
(b) MNAse-seq and dSMF profiles around chemically mapped positioned nucleosome dyads; 
(c) DNAse-seq and dSMF profiles around the top 20\% highly expressed genes in \textit{S. cerevisiae};
(d) DNAse-seq and dSMF profiles around the bottom 20\% expressed genes in \textit{S. cerevisiae};
(e) SMAC-seq profile around chemically mapped positioned nucleosomes dyads (shown is the  ``diamide 0 min rep2'' sample); 
(f) SMAC-seq profile around the top 20\% highly expressed genes in \textit{S. cerevisiae};
(g) SMAC-seq profile around the bottom 20\% expressed genes in \textit{S. cerevisiae};
(h) SMAC-seq correlates closely with both DNAse-seq and nucleosome occupancy profiling at the level of individual loci, and provides a combined readout of accessibility and nucleosome positioning. Shown is aggregate SMAC-seq signal along the genome (aggregated over 50-bp windows sliding every 5 bp; see Methods for details) together with DNAse-seq, nucleosome chemical mapping data, and transcriptional activity (measured by PRO-seq and PRO-cap). Large aggregate SMAC-seq signal enrichments match closely with DNAse accessibility peaks, while smaller aggregate SMAC-seq peaks are inversely correlated with positioned nucleosomes; 
(i) SMAC-seq profiles chromatin accessibility in repetitive regions of the genomes that are ``invisible'' to short reads. Shown is the telomeric region of chrXVI 
% \centerline{}
% \rightline{\textit{(legend continued on next page)}}
} 
\label{Fig1}
\end{figure*}
% \clearpage
% \let\thefootnote\relax\footnotetext{50-bp windows sliding every 5 bp; see Methods for details) together with DNAse-seq, nucleosome chemical mapping data, and transcriptional activity (measured by PRO-seq and PRO-cap). Large aggregate SMAC-seq signal enrichments match closely with DNAse accessibility peaks, while smaller aggregate SMAC-seq peaks are inversely correlated with positioned nucleosomes; 
% (i) SMAC-seq profiles chromatin accessibility in repetitive regions of the genomes that are ``invisible'' to short reads. Shown is the telomeric region of chrXVI.}

\begin{multicols}{2}

\noindent  chromatin states looks like at the kilo-- to multikilobase scale is largely unknown. Well established tools exist for mapping direct interactions between distal genomic elements in the context of three-dimensional genome architecture\cite{Lieberman2009,Fullwood2009,Mumbach2016} but they only capture pairwise interactions within large cell populations, not the chromatin state of the interacting regions, and they are often of limited resolution. Similar limitations apply to emerging methods for studying local chromatin structure\cite{Risca2017}. Super-resolution microscopy using highly multiplexed fluorescent probes is a powerful tool for revealing the folding of the chromatin fiber at the single cell level\cite{Boettiger2016}, but this approach is limited to a small number of loci and does not provide high-resolution base-pair level information about regulatory states. Cryo-electron microscopy-based approaches\cite{Ou2017} allow direct observations of the nucleosomal state along individual chromatin fibers, but at present the underlying DNA sequence cannot be linked to these measurements. More recently, long-read sequencing was used to map MNase cleavage events at multinucleosomal lengths\cite{NanoMNAse2018}; however, such approaches only provide information about two points in genomic space leaving the state of chromatin in the intervening sequence unknown.  

These technological limitations prevent the investigation of the manner in which chromatin  states of adjacent elements are correlated and how such coordination might play a role in gene regulation (for example, by creating self-reinforcing or mutually exclusive epigenetic states).  To address these technological limitations, we have developed SMAC-seq, a versatile, single molecule method that directly assays long-range nucleosome positioning and accessibility states within the chromatin fiber. We use this method to study chromatin architecture and co-accessibility states in the yeast \textit{Saccharomyces cerevisiae}, both under normal growth and conditions of cellular stress. We use SMAC-seq to assess the degree of coordination between the positions of nearby nucleosome particles, enumerate mutually exclusive regulatory states along individual loci, and observe coordinated changes in nucleosome positioning and chromatin accessibility upon transcriptional activation. SMAC-seq also enables high-resolution footprinting of transcription factor (TF) occupancy, and provides strand-specific information about the exposure of DNA bound to nucleosomes and other proteins. We expect future applications of, improvements on, and extensions of the SMAC-seq approach to enable novel insights into the dynamics of chromatin states in the context of a wide variety of experimental systems and biological questions.

\section*{Results}

\subsection*{SMAC-seq maps chromatin accessibility and nucleosome positioning at the multi-kb scale}

SMAC-seq is built on the conceptual foundations of the NOMe-seq assay\cite{Kelly2012,Nabilsi2014}. In its original form, NOMe-seq relies on the preferential modification of bases within accessible DNA with M.CviPI (a GpC-specific 5mC methyltransferase), followed by bisulfite conversion and Illumina-based sequencing readout. A more recent variation, dSMF\cite{Krebs2017} (\textbf{d}ual-enzyme \textbf{S}ingle \textbf{M}olecule \textbf{F}ootprinting), utilized an additional CpG-specific 5mC methyltransferase (M.SssI) to map promoter states in \textit{Drosophila} cells in finer detail. Such an approach is applicable in the \textit{Drosophila} context thanks to the absence of endogenous CpG methylation in flies.

We build on the dSMF approach by adding a third DNA-modifying enzyme, the m$^6$A methyltransferase EcoGII, and adapting the protocol to nanopore sequencing (Figure \ref{Fig1}A). Briefly, isolated nuclei are treated with the three enzymes, which preferentially modify DNA within accessible chromatin. High-molecular weight (HMW) DNA is then isolated and subjected to long-read single-molecule sequencing using the Oxford Nanopore platform. Using the ability of nanopore sequencing to directly read modified DNA bases\cite{Simpson2017,Rand2017}, we then obtain methylation maps for individual DNA molecules on a multikilobase scale (see the Methods section for details), which we then interpret in terms of chromatin accessibility.

The addition of m$^6$A methylation is of key importance to making SMAC-seq a high-resolution broadly applicable assay. Many eukaryote genomes are endogenously methylated at 5mC positions\cite{Feng2010,Zemach2010}. This is usually in a CpG context, in particular in metazoans, but it is not always limited to CpG dinucleotides. For example, in plants methylation in CHG and CHH contexts is also a common occurrence\cite{Lister2008}, thus even results obtained with M.CviPI alone are confounded by endogenous methylation. In addition, CpG and GpC dinucleotides are relatively rare in the genome; the average resolution achieved by the combination of CpG and GpC methyltransferases is $>$10 bp in \textit{Drosophila} and $\sim$15 bp in yeast. The addition of m$^6$A dramatically increases the resolution of SMAC-seq, to $\sim$3 bp in all main model organisms (Supplementary Figures \ref{FigS1}, \ref{FigS26}, and \ref{FigS27}).

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=17.75cm]{Fig2V8-part1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq provides a single-molecule linked-read view of the chromatin landscape}. 
(a) Unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrIII (``aggregate'' signal ``Sample 1''). 
(b) Unfiltered nanopore reads fully spanning a 6.6-kilobase neighborhood encompassing several genes on chrIV (``aggregate'' signal from ``Sample 1''). 
In both cases, accessibility is shown at 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the average SMAC-seq track.
% 500 reads sample, BI
}
\label{Fig2part1}
\end{figure*}

We initially developed and optimized the method in the budding yeast \textit{Saccharomyces cerevisiae} under normal growth conditions, as \textit{Saccharomyces cerevisiae} has no endogenous DNA methylation (allowing the simultaneous use of all three enzymes) and has a small genome ($\sim$12 Mbp), making it possible to achieve very high depth of nanopore sequencing coverage. To verify the specificity and efficiency of the enzymatic treatments, we carried out both dSMF experiments (i.e M.CviPI + M.SssI treatment) on yeast chromatin and M.CviPI + M.SssI + EcoGII reactions on naked yeast DNA (as well as untreated naked DNA controls), and subjected the resulting material to whole-genome bisulfite sequencing (see the Methods section for more details). We observe $\geq$95\% methylation in both the CpG and GpC contexts when starting with naked DNA, $\leq$10\% when working with chromatin, and nearly 0\% on naked DNA  (Supplementary Figure \ref{FigS62}). Comparing  dSMF profiles to DNase-seq and MNase-seq profiles around transcription start sites (TSSs) and positioned nucleosomes (obtained from previous H4S47C chemical mapping studies\cite{Brogaard2012}) revealed the expected nucleosome depletion and nucleosomal patterns, respectively (Figure \ref{Fig1}b-d).

We then carried out a SMAC-seq experiment on unsynchronized \textit{S. cerevisiae} cells using all three enzymes (experimental details described in the Methods section). We isolated HMW DNA and performed nanopore sequencing on the MinION platform. To call methylated bases, we first used Albacore for raw base calling, and then applied the Tombo\cite{Stoiber2017} algorithm in ``de novo'' mode (running on top of the Minimap aligner\cite{Li2016}) for ``resquiggling'' of raw nanopore signal to the genomic sequence, and calling of methylated bases. After mapping, we obtained reads with a median length of $\sim$1.5 kbp from this initial experiment (``Sample 1''; Supplementary Table \ref{TableS1}), which allows the capture of multiple promoter regions per fragment for much of the yeast genome (Supplementary Figure \ref{FigS61}). We also analyzed our initial dataset with Nanopolish\cite{Simpson2017}, an alternative algorithm for calling methylated bases, which is capable of identifying 5mC events in CpG and GpC context.

Unlike bisulfite-based conversion followed by sequencing-by-synthesis, nanopore-based direct measurement of nucleotide modifications does not provide unambiguous binary calls for methylated bases. Instead, methylation probabilities are obtained for each base. We therefore examined multiple strategies for binarizing these methylation calls within each read. We observed that while per-base methylation probabilities are skewed towards the two ends of the $[0,1]$ interval, a substantial number of bases have probabilities in between those extremes (Supplementary Figure \ref{FigS42}). We find that binarization at $p=0.5$ delivers the most optimal results (Supplementary Figures \ref{FigS40} and \ref{FigS41}) in terms of signal-to-noise ratio. Using a simple average of these binarized per-base methylation values, we compared SMAC-seq to dSMF, MNase-seq, and DNase-seq profiles, as well as signal from ChIP-seq for RNA Polymerase (Pol2) and transcription initiation factors around known chromatin features such as active promoters and well positioned nucleosomes (Figure \ref{Fig1}e-g, Supplementary Figures \ref{FigS8} and \ref{FigS9}). SMAC-seq faithfully reproduces nucleosomal positioning throughout the genome and nucleosome depletion around promoters, and, strikingly, has a larger observed dynamic range than dSMF data (possibly due to the higher mapping efficiency of long nanopore reads). Using comparisons with dSMF data we estimate the false positive rate of methylation base calling to be on the order of 20-25\% for Tombo and around 10-15\% for Nanopolish (Supplementary Figure \ref{FigS60}). We also examined potential sequence biases inherent to the combination of methylation enzymes and base calling algorithm. We find modest differences in methylation levels for different $k$-mers in the genome (less than two-fold for $k=6$; Supplementary Figures \ref{FigS3} and \ref{FigS4}).

In practice, the biologically relevant length scale of accessibility measurements is usually larger than an individual base. Furthermore, given noise intrinsic to single molecule sequencing methods, we reasoned that sharing methylation information between adjacent bases should improve the reliability of overall accessibility measurements. We therefore developed a simple Bayesian procedure to aggregate the methylation probabilities at individual bases and derive accessibility calls at the level of single reads over windows of arbitrary size (described in detail in the Methods section and thereafter referred to as ``aggregate'' signal as opposed to ``average'' signal, which refers to simple probability averaging and binarization).

We also note that sometimes we observe a subpopulation of reads that appear to be either entirely fully methylated or fully methylated over large segments of their length (Supplementary Figure \ref{FigS2} and Figure \ref{Fig2part1}a). We interpret these as originating from naked DNA molecules most likely deriving from dead cells. As such reads can confound many analyses, in particular when measuring coaccessibility within single reads, we devised a procedure for filtering them out (the resulting sets of reads are referred to as ``filtered reads''; (Supplementary Figures \ref{FigS38} and \ref{FigS39}). However, as discussed below, in certain situations chromatin is indeed largely nucleosome-free over specific regions in vivo, and in such cases filtering out fully methylated reads removes real biological signal. For these unique special case loci, we do not eliminate reads based on a very high fraction of accessibility (see Methods).% We note this issue prominently as it can potentially be a significant confounding factor in future experiments if not accounted for properly.

We then compared average SMAC-seq profiles against chemical maps of positioned nucleosomes (generated using H4S47C substitutions and copper-induced cleavage\cite{Brogaard2012}), DNase-seq, and maps of transcriptional activity at the level of individual loci in the genome (Figure \ref{Fig1}h). Qualitatively, we observe that large peaks in SMAC-seq signal profiles match very closely with DNase-seq peaks, while smaller ``bumps'' in the SMAC-seq signal profile are inversely correlate with positioned nucleosomes, consistent with labeling of linker DNA. We also observe positive correlation between average SMAC-seq methylation levels and DNase-seq and ATAC-seq coverage over promoter regions (Supplementary Figure \ref{FigS59}). Thus SMAC-seq simultaneously probes both regions of ``open'' chromatin, as well as the position of nucleosomes.

The long reads of nanopore sequencing allow SMAC-seq to provide accessibility maps for the whole yeast genome, not just for the portions of it that are uniquely mappable with short reads (Supplementary Figure \ref{FigS29}). For example, SMAC-seq maps chromatin and nucleosomes in the otherwise not uniquely mappable telomere of chrXVI (Figure \ref{Fig1}i), which is revealed to contain several active promoters and numerous well positioned nucleosomes. We also used SMAC-seq to characterize chromatin accessibility around multiple transposable elements (Supplementary Figures \ref{FigS28}), for several of which we observe open chromatin peaks around their promoters.

\subsection*{SMAC-seq provides single-molecule accessibility profiles on individual chromatin fibers}

\begin{figure*}
\begin{center}
\begin{minipage}[c]{0.80\linewidth}
\includegraphics[width=14.5cm]{Fig2V8-part2.png}
\end{minipage}\hfill
\begin{minipage}[c]{0.20\linewidth}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq's single-molecule readout provides insights into the distribution and relationship between mutually exclusive chromatin yeast rDNA states.} 
(a) SMAC-seq reveals the distribution of alternative chromatin states of rDNA arrays. Shown are all reads covering the \textit{RDN37-1} array in the \textit{RDN1} locus in the ``diamide 30 min rep1'' experiment (unfiltered reads, ``aggregate'' signal). See Supplementary Figures \ref{FigS56}--\ref{FigS80} for additional details. ChIP-seq and ChIP-exo tracks were generated by including and normalizing all multimappers rather than the usual unique-only policy (See the Methods section for more details). 
(b) Normalized mutual information profiles for the \textit{RDN37-1} array show anti-correlation between the accessibility peaks immediately upstream of the \textit{35S} TSS and the nucleosome-free state over the \textit{35S} transcriptional unit.
(c) High-resolution SMAC-seq profiles reveal regulatory protein footprints in the immediate vicinity of the 35S TSS and the Reb1 binding site in the rDNA NTS region (shown are 3000 randomly sampled reads using 10-bp aggregate SMAC-seq signal at 1-bp resolution).
\label{Fig2part2}}
\end{minipage}
\end{center}
\end{figure*}

Unlike short-read methods for probing chromatin, SMAC-seq allows the profiling of nucleosome positioning and open chromatin within individual long molecules. To demonstrate this capacity,  we investigated all reads spanning the 4-kb neighborhood around the centromere of chrIII (Figure \ref{Fig2part1}a). Centromeres in \textit{S. cerevisiae} are specified by a precisely defined sequence element, are occupied by a single strongly positioned nucleosome containing the H3 histone variant Cse4, and are also associated with the CBF3 complex, an essential component of the kinetochore\cite{Cole2011}. Yeast centromeric nucleosomes are thought to be nearly perfectly positioned\cite{Cole2011,Henikoff2014} and thus represent a good case system to study nucleosome positioning at the single-molecule level. We indeed observe strong nucleosomal positioning using SMAC-seq, with nearly all individual reads exhibiting the expected from the presence of a strongly positioned centromere nucleosomal pattern. We also find hints of substructure within the centromeric nucleosome in the form of accessibility traces inside the protected centromeric region and potential protection footprints in its immediate open chromatin vicinity. We find similarly strong positioning for most other centromeric nucleosomes (Supplementary Figures \ref{FigS14}-\ref{FigS13}), but not all (chrX, chrXII and chrXIII appear to be exceptions to this general strong footprinting pattern). 

We also illustrate the ability of SMAC-seq to capture accessibility at long-range multikilobase scales with an example from a more generic genomic neighborhood in Figure \ref{Fig2part1}b). This $\sim$6.6-kb span of chrIX contains five genes and three open chromatin regions, one of them fairly large and diffuse. In contrast to the more localized accessibility observed elsewhere, this region exhibits considerable heterogeneity in its accessibility suggesting a complex landscape of protein occupancy. 

We next asked if SMAC-seq could reveal binary states of chromatin accessibility. To approach this question we investigated SMAC-seq profiles at ribosomal DNA (rDNA). In yeast, rDNA is organized into multicopy arrays, each unit of which contains a copy of the 35S precursor pre-rRNA, transcribed by Pol I and later processed into mature 18S, 5.8S and 25S rRNAs, as well as a copy of the 5S RNA gene, transcribed by Pol III, and an origin of replication (ARS element in yeast) located in non-transcribed (NTS) regions of the array. Each array unit is $\sim$9.1 kb in length, and each cell's genome has an array of $\sim$150 copies of this unit\cite{Merz2008}. However, this number can vary from cell to cell, and the widely used \verb|sacCer3| \textit{S. cerevisiae} genome assembly only contains a single locus with two array copies. Chromatin structure at the rDNA locus has long been known to adopt two distinct conformations\cite{Conconi1989,French2003,Goetze2010}, depending on whether or not an individual unit is being transcribed. The active state is thought to be largely devoid of nucleosomes due to the extremely high levels of active transcription\cite{French2003}; the high-mobility group protein Hmo1 is proposed to replace nucleosomes\cite{Merz2008,Panday2016}. However, other studies have alternatively suggested that nucleosomes are found over actively transcribed rDNA arrays\cite{Jones2007}. Of note, rDNA indeed appears to be extremely accessible in short-read assays -- around half of reads in a typical ATAC-seq dataset in \textit{S. cerevisiae} are not uniquely mappable even though only a small fraction of the yeast genome consists of repetitive elements; these reads originate primarily from rDNA arrays (Supplementary Figure \ref{FigS54}a-e). The same phenomenon is also observed in other yeast species, such as \textit{Schizosaccharomyces pombe}  (Supplementary Figure \ref{FigS54}f). Active and inactive rDNA arrays are usually estimated to exist in roughly equal proportions in untreated, normally growing cells\cite{Jones2007}. However, methods to precisely observe these alternative states at the population level and relate them to sequence in fine detail have not been available. 

% We are able to address these questions using SMAC-seq, and are greatly aided by the fact that rDNA arrays exist in multiple copies, as this enabled us to achieve very high coverage even with reads fully spanning the whole 9.1-kb length of each of the arrays ($\geq$2000$\times$ in some samples). 

SMAC-seq reveals a striking picture of the two alternative mutually exclusive rDNA states at the single molecule level (Figure \ref{Fig2part2}a). About a quarter of full-length molecules exhibit near-full accessibility in the region spanning the 35S transcript, but not in the non-transcribed sequence (NTS) between the 35S and the 5S gene. The rest of the molecules show a typical nucleosomal state with several clearly accessible regions. We observe a broadly similar picture in all of our high-quality SMAC-seq samples (Supplementary Figures \ref{FigS56}--\ref{FigS80}). We note that while this picture is in contrast with the usually reported $\sim$50\% of arrays being fully accessible, it is possible that the fully accessible and nucleosomal states fragment differentially during DNA isolation; the fully accessible fraction appears to be in the 50\% range when a shorter window around the 35S promoter is examined (Figure \ref{Fig2part2}c). We also observe a localized accessible region just upstream of the 35S transcriptional unit that is present in the nucleosomal subpopulation but is not in an open state in the fully accessible population (in addition to a nearby open chromatin region present in all molecules), suggesting the possibility of a regulatory switch associated with that element. Finally, we also observe at least two (previously unreported) accessible regions located within the 35S transcriptional unit that exhibit strong accessibility in the nucleosome-protected fraction (Figure \ref{Fig2part2}a). 

To quantify the extent of (anti-)correlation between chromatin states, we developed a modified normalized mutual information (NMI) metric for assessing the degree of correlation between segments of the genome (see the Methods section for further details). NMI analysis of the rDNA arrays confirmed our observations of the inverse correlation between the 35S  open-chromatin state and the accessibility of the upstream element (Figure \ref{Fig2part2}b). 

What factors might be behind the observed chromatin state switch? Silencing of rDNA in yeast is thought to be mediated by a Sir2-containing complex called RENT \cite{Huang2003}, and a Reb1 binding site in NTS1 has been suggested to recruit corepressors to rDNA repeats\cite{Jones2007}. We took a higher-resolution view of NTS1 by integrating SMAC-seq data with available ChIP-exo data for Reb1 and transcription factor motif maps in the region. We find a clear pattern of protection from methylation around the Reb1 motif, which is concordant with ChIP-exo data (Figure \ref{Fig2part2}c), and we also observe signal consistent with footprinting from several other TF motifs. However, the anti-correlated accessibility profile seems to not be exclusively associated with Reb1 binding but rather with the region closer to the 35S TSS. It is likely that other proteins are responsible for establishing this state, but no currently annotated transcription factor recognition motifs are found in the underlying sequence.

\subsection*{SMAC-seq provides a high-resolution strand-specific view of protein occupancy on DNA}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig3-V3.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq provides a high-resolution strand-specific view of genomic occupancy by DNA-binding proteins and complexes}. 
(a-b) SMAC-seq allows for precise footprinting of transcription factor binding events. Shown is aggregate genome-wide SMAC-seq signal around occupied (as measured by ChIP-exo) Reb1 (a), and Rap1 (b) sequence recognition motifs. 
(c) SMAC-seq profiles around positioned nucleosome dyads reveal increased accessibility in the dyad and increased protection at the points of contact with the nucleosome (see Supplementary Figure \ref{FigS9} for additional details.)
(d) SMAC-seq provides a strand-specific view of nucleosome occupancy and reveals differential accessibility between the two DNA strands depending on their position on the nucleosomal particle.
(e-f) Coordination between the positions of individual nucleosomes at the level of single chromatin fibers. 
(e) Shown is the average normalized mutual information between each strongly or poorly positioned nucleosome in the yeast genome and its immediate genomic neighborhood (measured for windows of 10 bp length tiling at every genomic position centered on the nucleosome dyad).
(f) Shown is the average normalized mutual information between each $+1$ nucleosome and its immediate genomic neighborhood in highly expressed and in mostly silent genes (measured for windows of 10-bp length tiling at every genomic position centered on the $+1$ nucleosome dyad).
}
\label{Fig3}
\end{figure*}

We next asked if SMAC-seq can generally identify transcription factor footprints. We anticipated that the use of m$^6$A methylation ought to provide sufficient resolution to footprint most yeast TFs, which is confirmed by analysis of their consensus recognition motifs (Supplementary Figure \ref{FigS53}). Averaging genome-wide SMAC-seq profiles over occupied motifs indeed revealed strong protection footprints for several factors, such as Reb1, Rap1 and ORC1 (Origin of Replication Complex) as shown in Figure \ref{Fig3}b-c and Supplementary Figure \ref{FigS36}. Examination of individual sites confirmed these observations (Supplementary Figures \ref{FigS18}--\ref{FigS66}). We note that we did not observe strong footprinting for all TFs (e.g. Abf1 and Cbf1; Supplementary Figure \ref{FigS36}), suggesting that TF footprinting is probably dependent on the biophysical properties of individual factors.

We then explored global SMAC-seq signal around positioned nucleosomes (Figure \ref{Fig3}e-f and Supplementary Figure \ref{FigS9}). We find higher accessibility immediately at the dyad point, in contrast to the points on the nucleosome located two DNA helical turns away in each direction (Supplementary Figure \ref{FigS9}). The same pattern was observed for all nucleosomes irrespective of the strength of positioning; the difference between strongly positioned nucleosomes is the overall strength of protection from methylation (Figure \ref{Fig3}e). 

To further explore the limits of SMAC-seq's resolution, we studied methylation patterns around positioned nucleosomes in more detail (Figure \ref{Fig3}e-f and Supplementary Figure \ref{FigS9}). We find remarkably higher accessibility immediately at the dyad point, in contrast to the points of contact of the nucleosome with the DNA two helical turns away in each direction (Supplementary Figure \ref{FigS9}). The same pattern was observed for all nucleosomes irrespective of the strength of positioning (Figure \ref{Fig3}e). 

Because SMAC-seq directly maps accessibility independently on individual DNA strands, we next aimed to quantify strand-specificity in DNA accessibility within well-positioned nucleosomes. We observe a striking strand asymmetry in DNA accessibility around the nucleosome particle (Figure \ref{Fig3}f), especially within the dyad and at the points two helical turns away from the dyad. The magnitude of these differences in average methylation levels are similar to those observed between nucleosomes and flanking linker regions. Thus SMAC-seq reveals significant heterogeneity in DNA's accessibility potential within the nucleosomal particle, which has important implications for understanding how transcription factors interact with the genome in vivo, in particular in light of recent studies demonstrating that certain classes of TFs may preferentially occupy nucleosomal DNA\cite{Zhu2018}.

\subsection*{SMAC-seq reveals distal co-accessibility patterns in the genome}

We next examined co-accessibility patterns in the yeast genome. We first aimed to measure correlation between the positions of individual nucleosomes. Average NMI profiles centered on positioned nucleosomes reveal detectable correlation between nucleosome positions up to three to four nucleosomes away from an individual positioned nucleosome (Figure \ref{Fig3}g; Supplementary Figure \ref{FigS63}), with strongly positioned nucleosomal particles exhibiting stronger overall correlation over larger distances. These observations are consistent with a model whereby the restrictions on positioning that nucleosomes impose on each result in correlation between their protection footprints on short distances until random positional fluctuations of individual nucleosomes in the chromatin fiber eventually dephase this correlation signal on longer scales. 

We then examined co-accessibility patterns in the vicinity of promoters (Figure \ref{Fig3}h). Expressed yeast genes are characterized by a nucleosome-depleted/free region (NFR) upstream of the TSS and a well positioned $+1$ nucleosome. NMI profiles centered on the $+1$ nucleosome show significant differences between highly expressed and silent genes. While correlation patterns decay downstream of the TSS similarly for both groups of genes, highly expressed genes exhibit an inverse accessibility correlation pattern upstream of the TSS. The NMI profile for highly expressed genes also exhibits inverse accessibility correlation not just with the NFR but at a distance of at least one nucleosome beyond it.

Actively transcribed yeast genes often exist in a looped conformation, in which the promoter and termination regions are brought in physical proximity, potentially helping to enforce transcriptional directionality\cite{OSullivan2004,Tan2012}. Given this physical coupling, we wondered if correlation also exists between chromatin accessibility around gene ends. SMAC-seq data reveals low levels of correlation between the NFR and the accessible region in the 3$^\prime$ end of genes, and stronger correlation between positioned nucleosomes in these locations (Supplementary Figure \ref{FigS67}). The correlation between the accessibility in the NFR and the 3$^\prime$ end is increased for highly expressed genes and decreased for silent genes, suggesting that transcriptional activity and looping may help more strongly position nucleosomes at the beginning and end of transcribed regions; the decreased correlation between the nucleosome-depleted areas could be explained by transcriptional activity and dynamic regulatory occupancy leading to less stable protection patterns at each end that are accordingly less well correlated with one other.

We next assessed coordinated accessibility between yeast TSSs. To this end, we devised an explicit test of coordinated coaccessibility based on splitting reads into separate pieces, randomly reassembling them, then deriving an empirical coaccessibility distribution (see the Methods for details). Using this approach we identified 1,115 TSS pairs as significantly correlated out of 19,578 pairs covered with $\geq$100 reads in our initial sample (Supplementary Figure \ref{FigS87}). Of these, 560 were located a distance $\geq$ 1 kb from each other. An example of significantly coordinated accessibility (between the \textit{GAL10} and \textit{GAL1} genes) is shown in (Supplementary Figure \ref{FigS88}). As one possible mechanism for generating correlated accessibility is increased frequency of physical association in 3D space, we used publicly available Micro-C\cite{Hsieh2015} data to assess whether promoters exhibiting coordinated accessibility are more often physically interacting with each other. Indeed, we observe that significantly coaccessible promoters interact more frequently than non-coaccessible promoters at a similar distance (Supplementary Figure \ref{FigS86}).

% Finally, we note that when surveying NMI patterns genome-wide, we found several cases of striking switch-like linked accessibility patterns, in which a sharp discontinuity between a fully accessible state and a conventional nucleosomal state is observed (examples shown in Supplementary Figures \ref{FigS43}-\ref{FigS44}). These are primarily found in subtelomeric regions as well as around some repetitive elements. However, there is usually a clear pattern of strand-specificity to the switch-like behavior, thus at present the best explanation for these observations is that they are due to technical artifacts during read mapping and base calling, and we exclude them from further analysis.

\subsection*{SMAC-seq charts dynamic coordinated changes in chromatin accessibility in the course of the yeast stress response}

To monitor chromatin accessibility at long-range scales during dynamic changes in gene regulatory activity brought about by environmental stimuli, we carried out SMAC-seq experiments during a time course of diamide treatment. Diamide oxidizes thiols in proteins, resulting in the generation of disulfides and activation of the stress response pathway, leading to changes in the expression of several hundred genes\cite{Weiner2015}. Stress response in yeast is mediated through the action of the Hsf1 transcription factor (with the Msn2/4 proteins also playing a key role). Hsf1 binds to its cognate genes when activated upon stress leading to their transcriptional activation\cite{Morano2012}.

We performed SMAC-seq at 0, 30 and 60 minutes after diamide treatment of yeast cells, as well as RNA-seq, ATAC-seq, and ChIP-seq for RNA Pol2, the elongating version of Pol2 (Pol2pS2), and a V5-tagged version of HSF1 (RNA-seq and ATAC-seq data was also collected at 15 and 45 minutes; Figure \ref{Fig6}a). We observe several hundred genes exhibiting strong changes in gene expression during the time course (Supplementary Figure \ref{FigS58}), and strong induction of Hsf1 occupancy at hundreds of sites in the genome (Supplementary Figure \ref{FigS68}). SMAC-seq data at 30 minutes post-diamide treatment shows strong footprinting over Hsf1 motifs within induced Hsf1 binding sites.

\filllastline{We illustrate the dynamic patterns of chromatin accessibility that we observe upon diamide treatment using the \textit{TMA10} and \textit{HSP12} genes as examples in Figure \ref{Fig6}d-e, and multiple others in Supplementary Figures \ref{FigS69}-\ref{FigS85}. The \textit{TMA10} and \textit{HSP12} genes are strongly upregulated at 15 minutes after diamide treatment; TMA's expression subsequently declines somewhat and stabilizes (Figure \ref{Fig6}c) while that of \textit{HSP12} continues to increase up to the 45-minute time point. SMAC-seq reveals a relatively modest level of accessibility usptream of these genes before diamide treatment. However, at 30 minutes and upon Hsf1 binding, dramatic changes are evident. Nearby nucleosomes appear to be  evicted in many cells, and nucleosome depletion is also  observed at increased levels within the gene bodies, where ChIP-seq data for Pol2 and Pol2pS2 shows highly active transcription. At 60 minutes we observe dampening of this response, with the accessible fraction of reads decreasing, more strongly so around the TMA10 gene, whose expression decreases earlier than that of HSP12. Examination of NMI co-accessibility maps (Supplementary Figures \ref{FigS69}-\ref{FigS85}) frequently shows loss of correlation between positioned nucleosomes within and upstream of activated gene bodies as a result of response to diamide treatment, consistent with an increased movement of nucleosomal particles due to the} 

\end{multicols}

\begin{FPfigure}
\begin{center}
\includegraphics[width=18.5cm]{Fig6V3.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. 
(a) Experimental outline. Yeast cells were treated with diamide, then SMAC-seq and other functional genomic assays where carried out at 15- or 30-minute intervals. 
(b) Sites occupied by the HSF1 transcription factor upon its activation by the stress response pathway exhibit strong footprints in SMAC-seq data.
(c) Changes in the expression of the \textit{TMA10} gene upon diamide treatment
(d) Changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), and of chromatin accessibility at the single molecule level in the vicinity of the \textit{TMA10} gene during the diamide time course.
(e) Changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), and of chromatin accessibility at the single molecule level in the vicinity of the \textit{HSP12} gene during the diamide time course.
}
\label{Fig6}
\end{FPfigure}

\begin{multicols}{2}

\noindent activities of transcribing polymerase molecules and chromatin remodelers.

\section*{Discussion}

SMAC-seq is a novel single-molecule method for profiling chromatin accessibility within individual chromatin fibers at a multikilobase resolution that leverages the power of nanopore sequencing to detect DNA modifications and the preference of DNA methyltransferases for chemically modifying accessible DNA. We show that SMAC-seq generates accessibility signals similar to those of widely used methods such as DNase-seq and ATAC-seq while also opening new windows into the long-range structure of eukaryotic chromatin. SMAC-seq enables the simultaneous profiling of nucleosome positioning and accessible chromatin on a truly genome-wide scale (including repetitive regions that are poorly mappable with short reads), the assessment of the absolute distribution of chromatin accessibility states within a population of cells, and the identification of pairs of genomic loci that exhibit significant correlations in co-accessibility. 

Our initial work on SMAC-seq focused on applications in \textit{S. cerevisiae} because of its modest genome size. Extending SMAC-seq to larger genomes will require significantly increased sequencing throughput, or a method for selective enrichment of a collection of individual loci. Fortunately, both of these approaches appear feasible in the near term, as nanopore sequencing throughput is increasing rapidly, while methods for selective enrichment of genomic regions for nanopore sequencing are also becoming available\cite{Gabrieli2018}. Read length increases will also be useful, in particular for correlating the activity of distal regulatory elements to their cognate promoters, which can often be located many tens of kilobases apart in mammalian genomes.

Improvements to base calling accuracy constitute another major area of potential future advances. At the time of our analysis, Tombo was the only readily available algorithm for calling m$^6$A, but, as discussed above, its base calls exhibit on the order of 20-25\% error rates for all three bases combined. Results obtained with the other widely used methylation-aware base caller, Nanopolish, which is currently capable of calling only CpG and GpC methylation, show $\sim$15\% error on CpG and GpC calls. The major barrier to base-calling improvements is the lack of sets of ground truth controls that can be used to train base calling algorithms. Pools of DNA templates with individual modifications in well-defined yet highly diverse base pair contexts are the ideal training set for developing accurate such models. Alternatively, the introduction of tags bulkier than a simple methyl group\cite{TOPSeq} may provide a much stronger modulation of electric current through the nanopore than does simple methylation, enabling much more reliable modified base calls and accessibility evaluation.

We also anticipate an increase in the diversity of the DNA modifying enzymes available to carry out variations of SMAC-seq and related assays. In mammalian systems, endogenous 5mC methylation occurs primarily in the CpG context, thus it is not possible to use M.SssI for mapping accessibility, leaving only GpC as an option for traditional NOME-seq, bisulfite-based assays. Fortunately, the bulk of SMAC-seq's increase in resolution is derived from m$^6$A (Supplementary Figure \ref{FigS1}), so this will not be a significant obstacle to its widespread application in mammalian systems. However, there are species where m$^6$A occurs endogenously and is strongly correlated with patterns of accessibility and nucleosome positioning (e.g. \textit{Chlamydomonas}\cite{Fu2015} and \textit{Tetrahymena}\cite{Wang2017,Luo2018}). Alternative methylation strategies such as 4mC methyltransferases\cite{Timinskas1995}, cytidine deaminases\cite{Salter2016}, or the conversion of thymine to modified bases such as 5-hydroxymethyluracil (dhmU) \cite{Kawasaki2017} are areas of potential future exploration. 

Finally, we note that there is considerable scope for integration of SMAC-seq-type approaches with other measurements of the physical genome and the epigenome, especially once the improvements in accuracy outlined above are achieved. We anticipate the potential for obtaining simultaneous measurements of accessibility and nucleosomal positioning together with endogenous DNA methylation, general and specific protein occupancy, chromatin interactions, DNA replication, and other features on a multikilobase scale and within single molecules. In principle, similar approaches may also be applicable to individual RNA molecules. We expect long-read single-molecule approaches to provide an important new class of tools for the study of the functional and physical organization of genomes in the coming years.

\section*{Materials and Methods}

Except for when explicitly stated otherwise, all analyses were carried out using custom-written Python or R scripts. 

\subsection*{Cell lines and cell culture}

The BY4741 \textit{S. cerevisiae} strain (a kind gift from Ji-Ping Wang and Xiaozhong Wang) was used for all experiments except for Hsf1 ChIP-seq experiments where MS143 (H4S47C\_Hsf1-V5::HphMX6, this study) was used. MS143 was generated by PCR-based C-terminal tagging of Hsf1 with the V5 epitope. Hsf1-V5 tagging was confirmed by colony PCR and western blotting. For all experiments, except the initial one (``Sample 1''), cells were grown in YPD media (30$\,^{\circ}\mathrm{C}$) to OD$\sim$0.8 before collection. 

\subsection*{SMAC-seq experiments}

\subsubsection*{Enzymatic treatment of chromatin}

We developed and optimized SMAC-seq using the equivalent of $1 \times 10^6$ human cells, which in the case of \textit{S. cerevisiae} translates in to $2.5 \times 10^8$ (the size of the haploid human genome is $\sim$$3 \times 10^9$bp while that of \textit{S. cerevisiae} is $1.2 \times 10^6$bp). As yeast cells have a cell wall, we adapted the spheroplasting protocol previously used for carrying out ATAC-seq in yeast cells\cite{Schep2015} for our SMAC-seq experiments.

Yeast cells in log phase (OD$_{660} \leq 1.0$) were first centrifuged at 13,000 rpm for 1 minute, then washed with 100 $\mu$L Sorbitol Buffer(1.4 M Sorbitol, 40 mM HEPES-KOH pH 7.5, 0.5 mM MgCl$_2$), and centrifuged again at 13,000 rpm for 1 minute. Cells were then spheroplasted by resuspending in 200 $\mu$L Sorbitol Buffer with DTT added at a final concentration of 10 mM and 0.5 mg/mL 100T Zymolase, followed by incubating for 5 minutes at $30\,^{\circ}\mathrm{C}$ at 300 rpm in a Thermomixer. The pellet was centrifuged for 2 minutes at 5,000 rpm, washed in 100 $\mu$L Sorbitol Buffer, and centrifuged again at 5,000 rpm for 2 minutes.

Cells were then resuspended in 100 $\mu$L ice-cold Nuclei Lysis Buffer (10 mM Tris pH 7.4, 10 mM NaCl, 3 mM MgCl$_2$, 0.1 mM EDTA, 0.5\% NP-40) and incubated on ice for 10 minutes. Nuclei were then centrifuged at 5000 rpm for 5 min at $4\,^{\circ}\mathrm{C}$, resuspended in 100 $\mu$L cold Nuclei Wash Buffer (10 mM Tris pH 7.4, 10 mM NaCl, 3 mM MgCl$_2$, 0.1 mM EDTA), and centrifuged again at 5,000 rpm for 5 min at $4\,^{\circ}\mathrm{C}$. Finally, nuclei were resuspended in 100 $\mu$L M.CviPI Reaction Buffer (50 mM Tris-HCl pH 8.5, 50 mM NaCl, 10 mM DTT). 

Nuclei were then first treated with M.CviPI + EcoGII by adding 200 U of M.CviPI (NEB) and 200 units of EcoGII (NEB), SAM at 0.6 mM and sucrose at 300 mM, and incubating at $30\,^{\circ}\mathrm{C}$ for 7.5 min. After this incubation, 128 pmol SAM and another 100 U of enzymes were added, and a further incubation at $30\,^{\circ}\mathrm{C}$ for 7.5 min was carried out. Immediately after, M.SssI treatment followed, by adding 60 U of M.SssI (NEB), 128 pmol SAM, MgCl$_2$ at 10 mM and incubation at $30\,^{\circ}\mathrm{C}$ for 7.5 min. 

The reaction was stopped by adding an equal volume of Stop Buffer (20 mM Tris-HCl pH 8.5, 600 mM NaCl, 1\% SDS, 10 mM EDTA). 

\subsubsection*{High-molecular weight DNA isolation}

HMW DNA was isolated using the MagAttract HMW DNA Kit (Qiagen; cat $\#$ 67563) following the manufacturer's instructions.

\subsubsection*{Enzymatic treatment of naked DNA}

Naked DNA was treated under exactly the same conditions as chromatin except that the reaction volume and enzyme amounts were reduced in half. HMW DNA was purified as described above

\subsection*{SMAC-seq analysis}

\subsubsection*{Nanopore sequencing}

HMW DNA was converted into libraries using the Ligation Sequencing Kit 1D (Oxford Nanopore Technologies, SQK-LSK108) following the manufacturer's instructions. Nanopore sequencing was carried out on R9.4 MinION flowcells (Oxford Nanopore Technologies) for up to 48 hours. 

\subsubsection*{Nanopore base calling}

Nanopore events were converted to DNA sequence using Albacore (V2.3.3) using default settings. Reads were resquiggled using Tombo\cite{Stoiber2017}, version 1.3, using the \verb|sacCer3| reference genome. Methylated bases were identified using Tombo in the ``de novo model'' mode. 

\subsubsection*{Aggregation of accessibility information over multibasepair windows}

Even with the addition of m$^6$A methylation, the resolution of SMAC-seq still does not cover every nucleotide in the genome, and it varies substantially between different locations depending on local sequence content differences. In addition to that, nanopore base calling is still far from being a fully resolved problem, and even more so in methylation-aware mode. For these reasons, for many of the analyses described in this study we aimed at assigning aggregate accessibility scores over windows, taking the totality of the available evidence into account, thus obtaining more reliable, if coarser-grained, views of accessibility patterns along the genome. We used a Bayesian approach to carry out aggregation, as follows.

For a given window of width $w$ in the genome, specified by coordinates $c,i,i+w$ (where $c$ denotes the chromosome, and $i$ the leftmost coordinate of the window), and for all reads $r \in R_{c,i,i+w}$ fully spanning the window, we obtain all Tombo probabilities $p_{r,(c,j)}$ such that $j \in [i,i+w)$ for sequence contexts CpG, GpC and A on the corresponding genomic strand. We use a Beta prior $B$($\alpha$,$\beta$), with $\alpha = \beta = 10$, which we then updated based on each probability $p_{r,(c,j)}$ for all $j \in [i,i+w)$. The final binary accessibility score $p_{r,(c,i,i+w)}$ for read $r$ and window $c,i,i+w$ is determined by the final state of the prior. 

\subsubsection*{Read filtering}

As discussed above, we sometimes observe a population of reads that are fully methylated across their whole length or over large segments of it. There reads most likely derive from dead cells, as our initial experiment, which was carried out on a very dense yeast population containing a substantial number of dead cells, exhibited much higher proportion of such reads compared to subsequent experiments using early log-phase cells. In order to remove such potentially artifactual reads, ``filtered'' sets of reads were obtained by removing all reads containing a $\geq$1-kbp stretch that is $\geq$75\% methylated (while also filtering out reads shorter than 1 kb).

\subsubsection*{Read clustering}

For most analyses presented in this manuscript, the \verb|tglkmeans| package was used to cluster SMAC-seq reads (implemented in \verb|R|, \burl{https://bitbucket.org/tanaylab/tglkmeans}. In addition, the hierarchical clustering implementation in \verb|scipy| was also used in certain cases. 

\subsubsection*{Co-accessibility assessment using Normalized Mutual Information}

To evaluate co-accessibility patterns along the genome, we applied a Normalized Mutual Information as follows. Each chromosome in the genome $c$ was split into windows of size $w$. For each such window $(c,i,i+w)$, we identified the maximum range to the right of it, $(c,j,j+w)$ such that the span $(c,i,j+w)$ is covered by $\geq M$ reads. All reads spanning $(c,i,j+w)$ were then extracted and subsampled down to $M$ reads (usually $M = 100$, unless specified otherwise). Accessibility scores were then aggregated and binarized as described above for all windows located in the span $(c,i,j+w)$, and for all $M$ reads fully spanning it, resulting in a local co-accessibility matrix $LCM$ of size $M \times (j+w-i)/w$. We then calculated Normalized Mutual Information scores for each pair of columns $LCM_k$ and $LCM_l$ as follows:

\begin{equation}
\begin{split}
MI(LCM_k,LCM_l) & = p(0,0) \log_2{ \left(\cfrac{p(0,0)}{p_k(0)\,p_l(0)}\right) } \\
                & + p(1,1) \log_2{ \left(\cfrac{p(1,1)}{p_k(1)\,p_l(1)}\right) } \\
                & + p(0,1) \log_2{ \left(\cfrac{p(0,1)}{p_k(0)\,p_l(1)}\right) } \\
                & + p(1,0) \log_2{ \left(\cfrac{p(1,0)}{p_k(1)\,p_l(0)}\right) } \\
\end{split}
\end{equation}

MI scores were then normalized and rescaled in the interval $(-1,1)$:

\end{multicols}
\par\noindent\rule{\dimexpr(0.5\textwidth-0.5\columnsep-0.4pt)}{0.4pt}%
\rule{0.4pt}{6pt}

\begin{equation}
NMI(LCM_k,LCM_l) = \begin{cases}
\cfrac{MI(LCM_k,LCM_l)}{\sqrt{H(LCM_k)H(LCM_l)}} & \text{ for } p(0,0) + p(1,1) \geq 0.5 \\
-\cfrac{MI(LCM_k,LCM_l)}{\sqrt{H(LCM_k)H(LCM_l)}} & \text{ for } p(0,0) + p(1,1) < 0.5
\end{cases}
\end{equation}

\vspace{\belowdisplayskip}\hfill\rule[-6pt]{0.4pt}{6.4pt}%
\rule{\dimexpr(0.5\textwidth-0.5\columnsep-1pt)}{0.4pt}
\begin{multicols}{2}

Where $H$ refers to the entropy of each individual distribution. 

For computational efficiency, local NMI matrices were calculated for even-sized (50kb) evenly spaced (every 10kb) tiles of the genome. The entries of the general genome-wide NMI matrix were then calculated as the average of all local NMI matrices containing each entry.	

\subsubsection*{Testing for coordinated accessibility}

Coordinated accessibility was evaluated as follows. For each pair of locations $(c,i_1,i_1 + r_1)$ and $(c,i_2,i_2 + r_2)$ (usually $r_1 = r_2$), a minimum number of reads $N$ was required that fully spans the $(c,i_1,i_2 + r)$ interval. All such reads were then obtained for each pair, and then subsampled multiple times down to $N$ reads (in order not to introduce bias in coordinated accessibility tests arising due to differential read coverage between locations closer/further apart). For each subsampling, the fraction of accessible regions $p_1$ and $p_2$ was estimated for each of the two locations using the Bayesian procedure described above, as well as the distribution of joint accessibilities over the four states $(0,0)$, $(1,0)$, $(0,1)$, and $(1,1)$. The two halves of the reads were then virtually split in half and recombined for a total of $10^3$ random combinations. The empirical distribution $\mathcal{N}(\mu,\sigma)$ of the four states was then estimated from these random combinations, where $\mu = N*(p_{(0,0)} + p_{(1,1)})$  if $p_{(0,0)} + p_{(1,1)} > 0.5$ and $\mu = N*(p_{(1,0)} + p_{(1,0)})$  if $p_{(0,0)} + p_{(1,1)} \leq 0.5$. Empirical coordinated accessibility $p$-values were then estimated based on the observed counts $|(0,0)| + |(1,1)|$ if $|(0,0)| + |(1,1)| > 0.5*N$ or $|(0,1)| + |(1,0)|$ if $|(0,0)| + |(1,1)| \leq 0.5*N$. Bonferroni correction was applied to account for multiple hypothesis testing.

\subsection*{dSMF and Bisulfite sequencing}

Illumina measurements of CpG and GpC methylation levels were carried out using the PBAT\cite{Miura2012} with modifications. HMW DNA ($\sim$500 ng) was bisulfite converted using the EZ DNA Methylation-Lightning Kit (Zymo, Cat $\#$ D5030) by mixing 20 $\mu$L of purified DNA ($\sim$500 ng) with 130 $\mu$L DNA Methylation Lightning Conversion reagent and incubating at $98\,^{\circ}\mathrm{C}$ for 8 minutes and then at $64\,^{\circ}\mathrm{C}$ for 60 minutes. Bisulfite converted DNA was then cleaned up using the EZ DNA Methylation-Lightning Kit following manufacturer's instructions. 

First strand synthesis was carried out by mixing 20 $\mu$L bisulfite converted DNA, 19.75 $\mu$L H$_2$O, 5 $\mu$L 10$\times$ Blue Buffer (ThermoFisher), 1.25 $\mu$L 10 mM dNTP (NEB), and 4 $\mu$L custom-designed biotinylated adapter. Samples were then incubated at $94\,^{\circ}\mathrm{C}$ for 5 minutes, and at $4\,^{\circ}\mathrm{C}$ for 5 minutes, after which 1.5 $\mu$L Klenow (3$^\prime$ $\rightarrow$ 5$^\prime$ exo minus; MCLab) were added, and the reaction was incubated at $4\,^{\circ}\mathrm{C}$ for 15 minutes, at $37\,^{\circ}\mathrm{C}$ for 90 minutes, and at $70\,^{\circ}\mathrm{C}$ for 5 minutes. First-strand reaction cleanup was carried out using 50 $\mu$L AMPure XP beads (Beckman Coulter); DNA was eluted 50 $\mu$L EB buffer. 

Biotinylated DNA was captured on streptavidin beads. A total of 20 $\mu$L streptavidin Dybaneads M-280 (ThermoFisher) per sample were added to a PCR tube, separated on a magnet and then resuspended in 50 $\mu$L 2$\times$ BW(Li) buffer (6.3 g LiCl, 0.5 mL Tris-HCL pH 8.0, and 0.1 mL 500 mM EDTA for 50 mL total volume), to which the 50 $\mu$L of eluted first-strand reaction DNA was added. Beads were then incubated at room temperature for 30 minutes, washed with 180 $\mu$L 2$\times$ BW(Li) buffer, twice with 0.1 N NaOH (by resuspending well and incubating at room temperature for 2 minutes), washed again with 180 $\mu$L 2$\times$ BW(Li) buffer, then with 180 $\mu$L 10 mM Tris-HCL pH 7.5. 

Second-strand synthesis was carried out by resuspending streptavidin beads in the following reaction mix: 5 $\mu$L 10 $\times$ Blue Buffer, 1.25 $\mu$L 10 mM dNTPs, 39.75 $\mu$L H$_2$O, 4 $\mu$L custom-designed second-strand adapter. Samples were then incubated at $94\,^{\circ}\mathrm{C}$ for 5 minutes, and at $4\,^{\circ}\mathrm{C}$ for 5 minutes, after which 1.5 $\mu$L Klenow (3$^\prime$ $\rightarrow$ 5$^\prime$ exo minus) were added, followed by further incubation at $4\,^{\circ}\mathrm{C}$ for 15 minutes, at $37\,^{\circ}\mathrm{C}$ for 30 minutes, and at $70\,^{\circ}\mathrm{C}$ for 5 minutes. 

Beads were separated on magnet and the chase reaction was carried out by resuspending in a mix of 5 $\mu$L 10$\times$ Thermo Pol Buffer, 1.25 $\mu$L 10 mM dNTPs, 43.5  $\mu$L H$_2$O, and 1  $\mu$L \textit{Bst} DNA Polymerase Large Fragment (NEB). Samples were incubated at $65\,^{\circ}\mathrm{C}$ for 30 minutes, then again separated on magnet.

PCR was performed on beads in 50 $\mu$L reactions composed of 25 $\mu$L 2$\times$ NEB Next PCR Master Mix, 20 $\mu$L H$_2$O, 2.5 $\mu$L i7 and 2.5 $\mu$L i5 primers (both custom-designed), with initial extension at 72$\,^{\circ}\mathrm{C}$ for 3 min, denaturation at 98$\,^{\circ}\mathrm{C}$ for 30 sec, 15 cycles of 98$\,^{\circ}\mathrm{C}$ for 10 sec, 63$\,^{\circ}\mathrm{C}$ for 30 sec, and 72$\,^{\circ}\mathrm{C}$ for 30 sec, and final extension at 72$\,^{\circ}\mathrm{C}$ for 5 min. PCR reactions were cleaned up and size-selected using AMPure XP beads. 

Libraries were sequenced on Illumina NextSeq or MiSeq instruments, as 2$\times$75mers or 2$\times$300mers, respectively. 

\subsection*{dSMF data processing}

Bisulfite reads were trimmed using cutadapt (version 0.16) and Trim Galore (version 0.4.4), using the following settings (taking into account that the bisulfite sequencing libraries are generated with the PBAT protocol): \verb|--clip_R1 9| \verb|--clip_R2 9| \verb|--three_prime_clip_r1 6| \verb|--three_prime_clip_r2 6| \verb|--paired|. Trimmed reads were the mapped to the \verb|sacCer3| version of the yeast genome using Bismark\cite{Bismark2011} (version 0.19.0) with the following settings: \verb|--bowtie2| \verb|--pbat|. Methylation calls were extract using the \verb|bismark_methylation_extractor| program within Bismark and the following settings: \verb|-s| \verb|--no_overlap| \verb|--comprehensive| \verb|--merge_non_CpG| \verb|--cytosine_report| \verb|--CX|. 

\subsection*{ATAC-seq}

% ATAC-seq was carried out on the same nuclei isolated for SMAC-seq as described above (before resuspension in M.CviPI Reaction Buffer), by resuspending nuclei with 25 $\mu$L 2$\times$ TD buffer (20 mM Tris-HCl pH 7.6, 10 mM MgCl$_2$, 20\% Dimethyl Formamide), 2.5 $\mu$L transposase (custom produced) and 22.5 $\mu$L nuclease-free H$_2$0, and incubating at $37\,^{\circ}\mathrm{C}$ for 30 min in a Thermomixer at 1000 RPM. Transposed DNA was isolated using the DNA Clean \& Concentrator Kit (Zymo, cat $\#$ D4014) and PCR amplified as described before\cite{Corces2017}.

ATAC-seq was carried out on the same nuclei isolated for SMAC-seq as described above (before resuspension in M.CviPI Reaction Buffer), by resuspending nuclei with 25 $\mu$L 2$\times$ TD buffer (20 mM Tris-HCl pH 7.6, 10 mM MgCl$_2$, 20\% Dimethyl Formamide), 2.5 $\mu$L transposase (custom produced) and 22.5 $\mu$L nuclease-free H$_2$0, and incubating at $37\,^{\circ}\mathrm{C}$ for 30 min in a Thermomixer at 1000 RPM. Transposed DNA was isolated using the DNA Clean \& Concentrator Kit (Zymo, cat $\#$ D4014) and PCR amplified as described before\cite{Corces2017}. Libraries were then sequenced on a Illumina NextSeq instrument as 2$\times$36mers or  as 2$\times$75mers. 

\subsection*{ATAC-seq data processing}

Demultipexed fastq files were mapped to the \verb|sacCer3| assembly of the \textit{S. cerevisiae} genome as 2$\times$36mers using Bowtie\cite{Bowtie2009} with the following settings: \verb|-v 2| \verb|-k 2| \verb|-m 1| \verb|--best| \verb|--strata|. Duplicate reads were removed using \verb|picard|\verb|-tools| (version 1.99). 

\subsection*{ChIP-seq experiments}

Cell lysis and ChIP reactions were performed as previously described\cite{Hu2015} with minor modifications. Cells were fixed with 1\% formaldehyde for 20 minutes (Rpb1-CTD and Rbp1-CTD-S2P ChIP) or 30 minutes (Hsf1-V5 ChIP) and quenched with 0.125 M glycine for 5 minutes. A total of $\sim$50 ODs of cells were used per Rpb1-CTD or Rpb1-CTD-S2P ChIP and $\sim$300 ODs per Hsf1-V5 ChIP. Fixed cell were washed 2$\times$ in cold 1$\times$ PBS, pelleted and stored at $-80\,^{\circ}\mathrm{C}$. Pellets were lysed in 300 $\mu$L FA lysis buffer (50 mM HEPES--KOH pH 8.0, 150 mM NaCl, 1 mM EDTA, 1\% Triton X-100, 0.1\% sodium deoxycholate, 1 mM PMSF, Roche protease inhibitor) with $\sim$1 mL ceramic beads on a Fastprep-24 (MP Biomedicals). The entire lysate was then collected and adjusted to 1 mL with FA lysis buffer before sonication with a 1/8' microtip on a Q500 sonicator (Qsonica) for 14 minutes (10 seconds on, 20 seconds off). The sample tube was held in a $-20\,^{\circ}\mathrm{C}$ 80\% ethanol bath throughout sonication to prevent sample heating. After sonication, cell debris was pelleted and the supernatant was retained for ChIP. For each ChIP reaction, 30 $\mu$L Protein G Dynabeads (Invitrogen) were blocked (PBS + 0.5\% BSA), prebound with 5-10 $\mu$L antibody (8wG16 Rpb1-CTD, Abcam cat $\#$ ab817); 3E10 Rpb1-CTD-S2P, Milipore cat $\#$ 04-1571-1) or SV5-Pk1 (anti-V5, BioRad cat $\#$ MCA1360G)) and washed 1$\times$ with PBS before incubation with supernatant ($4\,^{\circ}\mathrm{C}$, overnight). Dynabeads were then washed (5 minutes per wash) 3$\times$ in FA lysis buffer, 3$\times$ in high-salt FA lysis buffer (50 mM Hepes–KOH pH 8.0, 500 mM NaCl, 1 mM EDTA, 1\% Triton X-100, 0.1\% sodium deoxycholate, 1 mM PMSF), 1$\times$ in ChIP wash buffer (10 mM Tris–HCl pH 7.5, 0.25 M LiCl, 0.5\% NP-40, 0.5\% sodium deoxycholate, 1 mM EDTA, 1 mM PMSF) and 1$\times$ in TE wash buffer (10 mM Tris–HCl pH 7.5, 1 mM EDTA, 50 mM NaCl). DNA was eluted from the beads in ChIP elution buffer (50 mM Tris–HCl pH 7.5; 10 mM EDTA; 1\% SDS) at $65\,^{\circ}\mathrm{C}$ for 20 min. Eluted DNA was incubated at $65\,^{\circ}\mathrm{C}$ overnight to reverse crosslinks, before treatment with RNAse A ($37\,^{\circ}\mathrm{C}$, 1 hour) and then Proteinase K ($65\,^{\circ}\mathrm{C}$, 2 hours). DNA was purified using the ChIP DNA Clean \& Concentrator kit (Zymo Research). Sequencing libraries were generated using the NEBNext Ultra II DNA Library Prep kit  (NEB Cat $\#$ E7645) and sequenced on a Illumina NextSeq instrument as 2$\times$36mers or as 2$\times$75mers. 

\subsection*{ChIP-seq data processing}

Demultipexed fastq files were mapped to the \verb|sacCer3| assembly of the \textit{S. cerevisiae} genome as 2$\times$36mers using Bowtie\cite{Bowtie2009} with the following settings: \verb|-v 2| \verb|-k 2| \verb|-m 1| \verb|--best| \verb|--strata|. Duplicate reads were removed using \verb|picard|\verb|-tools| (version 1.99). Hsf1 peaks were called using MACS2\cite{Feng2012} (version 2.1.0) with the following settings: \verb|-g 12000000| \verb|-f BAMPE|.

\subsection*{Multiread-preserving alignment and normalization}

Multiread-preserving alignment and track generation was carried out by mapping reads to the \verb|sacCer3| assembly of the \textit{S. cerevisiae} genome using Bowtie\cite{Bowtie2009} with the following settings: \verb|-v 2| \verb|-a| \verb|--best| \verb|--strata|. Each alignment was then given a weight inversely proportional to the number of locations that the read maps to i.e. each position's score was normalized to RPMs as follows:

\begin{equation}\label{multireads}
S_{c,i} = \cfrac{\displaystyle \sum_{R \in R_{c,i}}{\cfrac{1}{NH_R}}}{\cfrac{|R|}{10^6}}
\end{equation}

Where $NH_R$ is the number of locations in the genome a read maps to.

\subsection*{RNA-seq experiments}

Cells (1 mL) were pelleted and flash frozen in liquid N$_2$. Pellets were resuspended in 300 $\mu$L TRIzol and lysed with $\sim$1 mL ceramic beads on a Fastprep-24 (MP Biomedicals). Cell debris were pelleted and RNA was extracted from the supernatant using the Direct-Zol RNA Microprep Kit (Zymo Research). RNA-seq libraries were generated using the NEBNext Ultra Directional RNA Library Prep Kit (NEB Cat $\#$ E7420)

\subsection*{RNA-seq processing and gene expression quantification}

RNA-seq reads were mapped to the yeast genome as 1$\times$50mers (external datasets) or 2$\times$75mers (diamide experiments) using TopHat version 2.0.8\cite{Trapnell2012}. Gene-level  quantifications in FPKMs (Fragments Per Kilobase per Million mapped reads) were generated using Cufflinks version 2.0.2\cite{Trapnell2012}. The mean from all replicates was taken as the expression level for each gene for subsequent analyses. 

\subsection*{External sequencing datasets}

A number of previously published \textit{S. cerevisiae} genomics datasets were used in this study. ChIP-exo reads and called peaks for Abf1, Cbf1, Rap1 and Reb1 were downloaded from GEO accessions GSE93662 and GSE72106. ChIP-seq data for centromeric proteins was downloaded from GEO accessions GSE31466 and GSE51949. PRO-seq and PRO-CAP data was obtained from GEO accession GSE76142. ORC ChIP-seq data was downloaded from GEO accession GSE16926. DNase-seq was downloaded from GEO accession GSE69651 while DGF data was downloaded from SRA accession SRP000620. MNase-seq data was obtained from GEO accessions GSE26493 and GSE29292, TBP ChIP-seq from GSE44200, Rpb1 ChIP-seq from GSE93190, Rpb3 ChIP-seq from GSE74787, RPC128 ChIP-seq from GSE39566, and Mediator subunits ChIP-seq from GSE95051. RNA-seq data from accession GSE85590 was also used. Except where otherwise stated, raw reads were aligned using Bowtie\cite{Langmead2009} with the following settings: \verb|-v 2| \verb|--k 2| \verb|-m 1| \verb|--best| \verb|--strata|, with the addition of \verb|-X 1000| for paired-end reads. Paired-end reads were aligned as 2$\times$25mers), while single-end reads were aligned as 1$\times$36mers. PRO-seq and PRO-CAP data was aligned as 1$\times$16mers. 

\subsection*{Micro-C data and processing}

Micro-C data was downloaded from GEO accession GSE68016 and processed as described in the original publication\cite{Hsieh2015}.

\subsection*{Transcription factor motif mapping}

Transcription factor motif recognition sequences were mapped genome-wide using FIMO\cite{Grant2011} (version 4.11.2 of the MEME-Suite\cite{Bailey2009} using the CIS-BP database\cite{Weirauch2014} as a reference set of position weight matrices. 

\subsection*{Gene annotation update}

Publicly available gene models for \textit{S. cerevisiae} do not contain TSS and TTS information for a major fraction of genes in the genome, only including the coding (``CDS'') portions instead. As the omission of UTRs presents a problem for TSS- and TTS-centered analyses, we updated the existing gene models following the approach described previously\cite{Schep2015} and the \textit{S. cerevisiae} TIF-seq dataset from GEO accession GSE39128\cite{Pelechano2013}. New TSS and TTS positions were assigned to each gene for which such information was available based on the median UTR length as measured by TIF-seq.

\subsection*{Nucleosome positioning information}

H4S47C\cite{Brogaard2012,Ramachandran2015} chemical mapping data was downloaded form GEO accessions GSE59523 and GSE36063. H3Q85C\cite{Chereji2018} chemical mapping data was downloaded from GEO accession GSE97290. We used the nucleosome positioning calls obtained from the original 2012 Brogaard et al. study for our analyses, after transforming them from coordinates in the \verb|sacCer2| version of the \textit{S. cerevisiae} genome assembly to \verb|sacCer3| using the \verb|liftOver| function in the UCSC Genome Browser utilities toolkit.

\subsection*{Mappability tracks generation}

To evaluate unique read mappability, the whole genome was tiled with reads of given length at every position. The reads were then mapped back to the genome using the same settings used to map single-end ChIP-seq reads. For every position coverage by mapped reads was calculated, and mappability was scored as the ratio between read coverage and the read length used to tile the genome.

\section*{Author contributions}

Z.S, G.K.M and N.A.S.A. conceived and designed the study. Z.S, G.K.M and N.A.S.A. performed initial experiments. Z.S., M.P.S. and G.K.M. performed diamide time course experiments. G.K.M. and Z.S. analyzed data. W.J.G. and A.K. supervised the study. G.K.M., Z.S., W.J.G. and A.K. wrote the manuscript with input from all authors. 

\section*{Acknowledgments}

This work was supported by NIH grants (P50HG007735, RO1 HG008140, U19AI057266 and UM1HG009442 to W.J.G., 1UM1HG009436 to W.J.G. and A.K., 1DP2OD022870-01 and 1U01HG009431 to A.K.), the Rita Allen Foundation (to W.J.G.), the Baxter Foundation Faculty Scholar Grant, and the Human Frontiers Science Program grant RGY006S (to W.J.G). W.J.G is a Chan Zuckerberg Biohub investigator and acknowledges grants 2017-174468 and 2018-182817 from the Chan Zuckerberg Initiative. Z.S. is supported by EMBO Long-Term Fellowship EMBO ALTF 1119-2016 and by Human Frontier Science Program Long-Term Fellowship HFSP LT 000835/2017-L G.K.M. is supported by the Stanford School of Medicine Dean's Fellowship. N.A.S.A. is funded by the Department of Defense through a National Defense Science and Engineering Grant and by a Stanford Graduate Fellowship. The authors would also like to thank members of the Greenleaf and Kundaje labs for their helpful suggestions and discussions on the subject over the course of the study.

\begin{thebibliography}{100}

% \section*{References}

\input{references-V4}

\end{thebibliography}

\end{multicols}

% \end{document}

\clearpage

\setcounter{table}{0}
\renewcommand{\tablename}{Supplementary Table}
\setcounter{figure}{0}
\renewcommand{\figurename}{Supplementary Figure}

\setcounter{page}{1}
\renewcommand\thepage{{SM }\arabic{page}}

\begin{center}
% {\LARGE \textbf{\begin{spacing}{1.1}The Cost of a Gene. \\ Supplementary Materials\end{spacing} }}
{\LARGE \textbf{Supplementary Materials}}
\end{center}


\section*{Supplementary Tables}

% \begin{scriptsize}
\begin{center}
\begin{longtable}{m{4cm}m{2.5cm}m{2.5cm}m{2.5cm}m{2.5cm}}
\caption{Mapping statistics for SMAC-seq datasets in this study}\\
\hline
sample & number reads & total bases & mean read length & median read length \\
\hline
\endfirsthead
\multicolumn{5}{c}%
{\tablename\ \thetable\ -- \textit{Continued from previous page}} \\
\hline
sample & number reads & total bases & mean read length & median read length \\
\hline
\endhead
\hline \multicolumn{5}{r}{\textit{Continued on next page}} \\
\endfoot
\hline
\endlastfoot
Sample 1 & 847,262 & 3,475,258,237 & 4,102 & 1,474 \\
diamide 0 min rep1 & 2,462,311 & 4,018,221,647 & 1,632 & 722 \\
diamide 30 min rep1 & 1,127,937 & 3,376,196,280 & 2,993 & 1,254 \\
diamide 60 min rep1 & 1,543,217 & 4,187,296,821 & 2,713 & 1,171 \\
diamide 0 min rep2 & 1,375,008 & 2,232,750,288 & 1,624 & 678 \\
diamide 30 min rep2 & 2,266,191 & 6,046,929,720 & 2,668 & 977
\label{TableS1}
\end{longtable}
\end{center}
% \end{scriptsize}

\begin{center}
\begin{longtable}{m{4.5cm}m{1.5cm}m{1cm}m{1cm}m{2cm}m{2cm}}
\caption{Mapping and QC statistics for ATAC-seq datasets used in this study}\\
\hline
Dataset & Complexity & TSS ratio & Read Length & uniquely mapped deduplicated reads & Raw fragments \\
\hline
\endfirsthead
\multicolumn{6}{c}%
{\tablename\ \thetable\ -- \textit{Continued from previous page}} \\
\hline
Dataset & Library complexity & TSS ratio & Read Length & uniquely mapped deduplicated reads & Raw fragments \\
\hline
\endhead
\hline \multicolumn{6}{r}{\textit{Continued on next page}} \\
\endfoot
\hline
\endlastfoot
L464 Diamide 0 min rep1 & 0.79 & 1.38 & 2 $\times$ 36 & 2,284,992 & 2,363,608 \\
L465 Diamide 0 min rep2 & 0.79 & 1.36 & 2 $\times$ 36 & 2,383,094 & 2,409,446 \\
L466 Diamide 15 min rep1 & 0.77 & 1.49 & 2 $\times$ 36 & 1,907,268 & 1,760,891 \\
L467 Diamide 15 min rep2 & 0.75 & 1.42 & 2 $\times$ 36 & 3,415,058 & 3,058,834 \\
L468 Diamide 30 min rep1 & 0.72 & 1.49 & 2 $\times$ 36 & 3,223,114 & 3,414,835 \\
L469 Diamide 30 min rep2 & 0.74 & 1.41 & 2 $\times$ 36 & 3,081,970 & 2,846,193 \\
L470 Diamide 45 min rep1 & 0.79 & 1.37 & 2 $\times$ 36 & 2,411,938 & 2,457,100 \\
L471 Diamide 45 min rep2 & 0.75 & 1.40 & 2 $\times$ 36 & 3,135,316 & 2,885,651 \\
L472 Diamide 60 min rep1 & 0.74 & 1.45 & 2 $\times$ 36 & 3,205,726 & 3,120,294 \\
L473 Diamide 60 min rep2 & 0.76 & 1.43 & 2 $\times$ 36 & 2,387,398 & 2,244,141
\label{TableS2}
\end{longtable}
\end{center}

\clearpage

\begin{center}
\begin{longtable}{m{6cm}m{1.5cm}m{1cm}m{2cm}m{1.75cm}}
\caption{Mapping and QC statistics for ChIP-seq datasets used in this study}\\
\hline
Dataset & Complexity & Read Length & Uniquely mapped deduplicated reads & Raw fragments \\
\hline
\endfirsthead
\multicolumn{5}{c}%
{\tablename\ \thetable\ -- \textit{Continued from previous page}} \\
\hline
Dataset & Complexity & Read Length & Uniquely mapped deduplicated reads & Raw fragments \\
\hline
\endhead
\hline \multicolumn{5}{r}{\textit{Continued on next page}} \\
\endfoot
\hline
\endlastfoot
L482 Diamide 0 min Input & 0.88 & 2 $\times$ 36 & 8,776,410 & 5,488,322 \\
L483 Diamide 0 min Pol2 & 0.87 & 2 $\times$ 36 & 4,944,766 & 2,894,414 \\
L484 Diamide 0 min Pol2pS2 & 0.80 & 2 $\times$ 36 & 5,834,864 & 3,434,723 \\
L485 Diamide 0 min HSF1-V5 Input & 0.90 & 2 $\times$ 36 & 6,089,572 & 3,850,484 \\
L486 Diamide 0 min HSF1-V5 & 0.93 & 2 $\times$ 36 & 2,540,876 & 1,726,223 \\
L487 Diamide 30 min Input & 0.86 & 2 $\times$ 36 & 10,052,178 & 6,212,240 \\
L488 Diamide 30 min Pol2 & 0.84 & 2 $\times$ 36 & 6,763,128 & 3,961,384 \\
L489 Diamide 30 min Pol2pS2 & 0.84 & 2 $\times$ 36 & 6,332,462 & 3,901,689 \\
L490 Diamide 30 min HSF1-V5 Input & 0.92 & 2 $\times$ 36 & 4,587,466 & 2,831,871 \\
L491 Diamide 30 min HSF1-V5 & 0.89 & 2 $\times$ 36 & 3,324,498 & 2,054,869 \\
L492 Diamide 60 min Input & 0.90 & 2 $\times$ 36 & 5,812,736 & 3,539,630 \\
L493 Diamide 60 min Pol2 & 0.90 & 2 $\times$ 36 & 3,774,106 & 2,204,126 \\
L494 Diamide 60 min Pol2pS2 & 0.85 & 2 $\times$ 36 & 4,873,244 & 2,924,345 \\
L495 Diamide 60 min HSF1-V5 Input & 0.87 & 2 $\times$ 36 & 7,683,586 & 4,679,564 \\
L496 Diamide 60 min HSF1-V5 & 0.91 & 2 $\times$ 36 & 1,094,048 & 698,664
\label{TableS3}
\end{longtable}
\end{center}

\begin{center}
\begin{longtable}{m{3.5cm}m{1.5cm}m{1cm}m{1.5cm}m{1.5cm}m{1.5cm}m{1.5cm}m{1.75cm}}
\caption{Mapping and QC statistics for RNA-seq datasets used in this study}\\
\hline
Dataset & Complexity & Read Length & Unique & Unique Splices & Multi & Multi Splices & Raw fragments \\
\hline
\endfirsthead
\multicolumn{8}{c}%
{\tablename\ \thetable\ -- \textit{Continued from previous page}} \\
\hline
Dataset & Complexity & Read Length & Unique & Unique Splices & Multi & Multi Splices & Raw fragments \\
\hline
\endhead
\hline \multicolumn{8}{r}{\textit{Continued on next page}} \\
\endfoot
\hline
\endlastfoot
Diamide 0 min & 0.20 & 2 $\times$ 75 & 30,692,672 & 461,804 & 3,974,211 & 38,859 & 18,719,790 \\
Diamide 15 min & 0.21 & 2 $\times$ 75 & 26,991,043 & 182,208 & 3,217,924 & 18,788 & 15,960,277 \\
Diamide 30 min & 0.19 & 2 $\times$ 75 & 34,711,193 & 201,334 & 3,716,196 & 21,450 & 21,717,800 \\
Diamide 45 min & 0.22 & 2 $\times$ 75 & 28,668,901 & 116,773 & 3,103,601 & 19,301 & 17,832,222 \\
Diamide 60 min & 0.31 & 2 $\times$ 75 & 14,991,437 & 71,198 & 1,759,048 & 10,455 & 13,548,619 \\
Diamide Hsf-V5 0min & 0.19 & 2 $\times$ 75 & 32,981,363 & 519,649 & 4,320,132 & 35,513 & 21,535,487 \\
Diamide Hsf-V5 30min & 0.19 & 2 $\times$ 75 & 33,505,595 & 211,992 & 3,973,770 & 23,016 & 22,974,020 \\
Diamide Hsf-V5 60min & 0.19 & 2 $\times$ 75 & 40,290,524 & 196,082 & 4,698,437 & 29,628 & 28,894,301
\label{TableS4}
\end{longtable}
\end{center}

\clearpage

\section*{Supplementary Figures}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS1-resolution.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Resolution of the SMAC-seq assay in its current form and some potential other versions of it in the main model organism systems}. Note that where endogenous methylation confounds readout of accessibility, the corresponding combination of sequence contexts has been omitted from the plot. Also note that while m$^6$A has been reported to be present in the genomes of \textit{Arabidopsis}\cite{Liang2018} and \textit{C. elegans}\cite{Greer2015}, it is generally found at low levels ($\leq 1$\%), and is not strongly as strongly correlated with open chromatin and nucleosome positioning as it is in some other eukaryotes such as \textit{Chlamydomonas}\cite{Fu2015}, thus its utility for accessibility profiling is not altered significantly. Nevertheless, a universally applicable version of SMAC-seq that is minimally confounded by endogenous methylation status in all eukaryotes will probably require the use of different methyltransferases (once they become available as efficient recombinant enzymes), for example, ones depositing the 4mC mark, which is what the ``C'' sequence context shown here corresponds to.
(a) \textit{Saccharomyces cerevisiae} (complete absence of endogenous methylation);
(b) \textit{Caenorhabditis elegans} (no endogenous 5mC, small amounts of endogenous m$^6$A); 
(c) \textit{Drosophila melanogaster} (no significant endogenous 5mC or m$^6$A methylation); 
(d) \textit{Arabidopsis thaliana} (endogenous 5mC in CpG, CHG and CHH contexts, small amounts of endogenous m$^6$A);
(e) \textit{Homo sapiens} (endogenous 5mC in CpG contexts, small amounts of endogenous 5mC in CHG and CHH contexts);
(f) \textit{Mus musculus} (endogenous 5mC in CpG contexts, small amounts of endogenous 5mC in CHG and CHH contexts).
}
\label{FigS1}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=12cm]{FigS61-Distance.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Distance between regulatory elements (i.e. promoters in the case of \textit{S. cerevisiae}) and distribution of obtained nanopore read lengths}. 
(a) Distance between annotated promoters in \textit{Saccharomyces cerevisiae};
(b) Distribution of nanopore read lengths (in the ``Sample 1'' experiment).
}
\label{FigS61}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=8cm]{FigS62-Controls.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Global methylation levels in yeast dSMF experiments and in positive and negative controls (measured by bisulfite sequencing)}. DNA from SMAC-seq experiments was subjected to Illumina bisulfite sequencing using the PBAT protocol. In parallel, naked genomic DNA (gDNA) was either treated with all three enzymes under the same conditions as in the SMAC-seq protocol or it was left untreated. These samples were also subjected to Illumina PBAT bisulfite sequencing. Shown are the global methylation levels for each sample. 
}
\label{FigS62}
\end{figure*}

\clearpage 

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=17cm]{FigS26-different-enzymes.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Impact of the addition of m$^6$A on assay resolution}. Shown is the raw unfiltered nanopore read coverage around a divergent promoter region on chrXII, considering only CG, only GC, only m$^6$A, and all bases at 1-bp resolution as well as all bases at averaged and aggregated 10-bp resolution.
}
\label{FigS26}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=16.5cm]{FigS27-different-enzymes-nucleosome.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Impact of the addition of m$^6$A on assay resolution}. Shown is the raw unfiltered nanopore read coverage around a strongly positioned +1 nucleosome, considering only CG, only GC, only m$^6$A, and all bases at 1-bp resolution as well as all bases at averaged and aggregated 10-bp resolution.
}
\label{FigS27}
\end{figure*}

% \end{document}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS42-probabilities-distribution.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Distribution of methylation call probabilities}. Shown is the distribution of Tombo ``alternative model'' probabilities for all positions, and each of the three sequence contexts, as well as the distribution after filtering potential poor quality/non-chromatinized reads (see further below for more details).
}
\label{FigS42}
\end{figure*}

\clearpage

\begin{center}
\includegraphics[width=15cm]{FigS40-binarization-ab.png}
\newpage
\end{center}
\begin{figure*}
\begin{center}
\includegraphics[width=15cm]{FigS40-binarization-cd.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Transformation of raw methylation probabilities into binary methylation calls}. Shown is raw unfiltered SMAC-seq single-molecule data over a strongly positioned nucleosome on chrXVI (1-bp resolution).
(a) raw Tombo alternative model methylation probabilities; 
(b) $p < 0.5$ thresholding; 
(d) $p < 0.3$ thresholding; 
(d) $p < 0.1$ thresholding.
}
\label{FigS40}
\end{figure*}

\clearpage

\begin{center}
\includegraphics[width=15cm]{FigS41-binarization-ab.png}
\end{center}
\newpage
\begin{figure*}
\begin{center}
\includegraphics[width=15cm]{FigS41-binarization-cd.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Transformation of raw methylation probabilities into binary methylation calls}. Shown is raw unfiltered SMAC-seq single-molecule data over a $\sim$2.8kb locus on chrIV (10-bp average for all panels).
(a) raw Tombo alternative model methylation probabilities; 
(b) $p < 0.5$ thresholding; 
(d) $p < 0.3$ thresholding; 
(d) $p < 0.1$ thresholding.
}
\label{FigS41}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS2-read_length_vs_methylation.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Absence of strong correlation between nanopore sequencing read length and methylation status}. Shown is the fraction of bases within each read that is scored as methylated. 
(a) CG positions only. 
(b) GC positions only. 
(c) m$^6$A positions only. 
(d) All positions.}
\label{FigS2}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS3-seq-bias.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Examination of enzymatic/methylated base calling bias}. Shown is the ratio of observed versus expected fraction of methylated bases for each sequence context of size $k$, calculated as follows:
\newline
\centerline{$f_{\mbox{obs/exp},k} = \displaystyle\cfrac{k_m/k_u}{\sum_k{k_m}/\sum_k{k_u}}$}
\newline
where $k_m$ refers to the number of bases called as methylated across all reads and $k_m$ refers to the number of bases called as unmethylated .
(a) CG positions only, $k=4$. 
(b) GC positions only, $k=4$. 
(c) A positions only, $k=4$. 
(d) CG positions only, $k=6$. 
(e) GC positions only, $k=6$. 
(f) A positions only, $k=6$. 
(g) CG positions only, $k=8$. 
(h) GC positions only, $k=8$. 
(i) A positions only, $k=8$.
}
\label{FigS3}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS4-seq-bias-vs-GC-content.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Relationship between local GC content and enzymatic/methylated base calling bias}. Shown is the ratio of observed versus expected fraction of methylated bases for each sequence context of size $k$, calculated as follows:
\newline
\centerline{$f_{\mbox{obs/exp},k} = \displaystyle\cfrac{k_m/k_u}{\sum_k{k_m}/\sum_k{k_u}}$}
\newline
where $k_m$ refers to the number of bases called as methylated across all reads and $k_m$ refers to the number of bases called as unmethylated .
(a) CG positions only, $k=4$. 
(b) GC positions only, $k=4$. 
(c) A positions only, $k=4$. 
(d) CG positions only, $k=6$. 
(e) GC positions only, $k=6$. 
(f) A positions only, $k=6$. 
(g) CG positions only, $k=8$. 
(h) GC positions only, $k=8$. 
(i) A positions only, $k=8$.
}
\label{FigS4}
\end{figure*}

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=14cm]{FigS5-neighboring-sequences.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf \hl{CAPTION}}. .
% (a) CG positions only;
% (b) GC positions only; 
% (c) A positions only. 
% }
% \label{FigS5}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS30-plus-minus-strand-correlation-10bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 10 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS30}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS31-plus-minus-strand-correlation-20bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 20 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS31}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS32-plus-minus-strand-correlation-30bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 30 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS32}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS33-plus-minus-strand-correlation-50bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 50 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS33}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS34-plus-minus-strand-correlation-75bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 75 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS34}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.70\linewidth}
% \includegraphics[width=15cm]{FigS35-plus-minus-strand-correlation-100bp.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.30\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Similarity between methylation percentages observed on the forwards/plus and reverse/minus strand of the genome}. The nanopore sequencing reads generated in this study are stranded, and provide information about methylation status only on one of the two strands. In addition, unlike CpG and GpC methylation events, which are inherently symmetric due to the palindromic nature of the two sequence contexts, m$^6$A methylation is asymmetric, and methylation estimates relying on it are expected to be more sensitive to local variations in sequence composition. However, note that because different molecules are sequenced for each strand, the estimates for the two strands are not expected to be exactly the same even for the CpG and GpC sequence contexts. Shown in this figure is the correspondence between methylation fractions calculated over windows of size 100 bp. The plots on the left show all windows in the genome, the plots on the right show only windows with an estimated methylated fraction of reads $\geq$0.5.
% (a) CG positions only;
% (b) GC positions only;
% (c) m$^6$A positions only;
% (d) all positions.
% }
% \label{FigS35}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS60-dSMF-vs-SMAC.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Comparison of dSMF results and different approaches to methylation-aware base calling on SMAC-seq data}. Shown is the inverse of the methylated fraction of nucleotides around TSSs of all, highly express (top quantile) and low expression-level (bottom quantiles) yeast genes (unfiltered ``Sample 1'' dataset). Note that the different panels are not drawn to the same scale.
(a) dSMF; 
(b) SMAC-seq data using Nanopolish methylation base-calling on CG and GC nucleotides;
(c) SMAC-seq data using Tombo methylation base-calling, CG positions only;
(d) SMAC-seq data using Tombo methylation base-calling, GC positions only;
(e) SMAC-seq data using Tombo methylation base-calling, m6A positions only.
}
\label{FigS60}
\end{figure}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=13cm]{FigS8-positioned-nucleosomes.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation of SMAC-seq (Figure \ref{Fig1}e) with other measures of chromatin structure around the dyad centers of positioned nucleosomes in the \textit{S. cerevisiae} genome}.
(a) H4S47C and H3Q85C nucleosome chemical mapping; 
(b) MNaseq-seq;
(c) DNase-seq;
(d) Digital Genomic Footprinting (DGF, 5$^\prime$ ends of deeply sequenced DNAseq-seq data).
}
\label{FigS8}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS7-TSSs.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation of SMAC-seq (Figure \ref{Fig1}f and g) with other types of functional genomic measurements of chromatin structure, protein occupancy and transcriptional activity around TSSs}. Shown is average coverage over all \textit{S. cerevisiae} genes, for the most highly expressed 20\% of genes (``quantile1''), and for the bottom 20\% of genes (``quantile5'').
(a) DNase-seq;
(b) H4S47C nucleosome chemical mapping; 
(c) H3Q85C nucleosome chemical mapping;
(d) TBP ChIP-seq;
(e) Rpb1 ChIP-seq;
(f) Rpb3 ChIP-seq;
(g) Med3 ChIP-seq;
(h) PRO-cap;
(i) PRO-seq.
}
\label{FigS7}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18cm]{FigS59-DNAse-ATAC-vs-SMAC.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation of SMAC-seq signal with ATAC-seq and DNA-seq signal over yeast promoters}. Shown is the average methylation over the TSS $\pm$ 200 bp for SMAC-seq and RPM (Reads Per Million mapped reads) values for DNase-seq (a) and ATAC-seq (b). Note that the DNase-seq dataset is obtained from an external study while the SMAC-seq and ATAC-seq ones are from the same sample (``diamide 0 min rep1'').
}
\label{FigS59}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=16cm]{FigS9-nucleosome-structure-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq data displays higher methylation propensity in more exposed parts of the nucleosomal particle}.
(a) Structure of the eukaryotic nucleosome; 
(b) High-resolution (50-bp radius) view of chemical nucleosome mapping data relative to nucleosome dyads;
(c) High-resolution (50-bp radius) view of SMAC-seq data relative to nucleosome dyads;
(d) Strand-specific  (100-bp radius) view of SMAC-seq data relative to nucleosome dyads;
(e) High-resolution (50-bp radius) view of DGF cleavage profiles relative to nucleosome dyads.
}
\label{FigS9}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS28-unmappable-regions-coverage.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq provides coverage of areas of the genome that cannot be uniquely mapped using short reads}. Four different repetitive regions are shown (a, b, c and d) together with short read coverage for DNase-seq and chemical nucleosome mapping data and unique mappability tracks for 36mers and 100mers. 
}
\label{FigS28}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\begin{minipage}[c]{0.60\linewidth}
\includegraphics[width=11cm]{FigS29-mappability-vs-SMAC-seq.png}
\end{minipage}\hfill
\begin{minipage}[c]{0.40\linewidth}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq provides coverage of areas of the genome that cannot be uniquely mapped using short reads}. 
(a) Average short-read mappability around TSSs of annotated transposable elements in the \textit{S. cerevisiae} genome;
(b) SMAC-seq signal around TSSs of annotated transposable elements in the \textit{S. cerevisiae} genome.
}
\label{FigS29}
\end{minipage}
\end{center}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=17cm]{FigS38-filtering.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Removal of potentially artifactual high-methylation reads}. Shown is unfiltered SMAC-seq data and the same locus after removal of all reads containing a 1-kb $\geq$75\% methylated stretch (average accessibility over 10 bp).
}
\label{FigS38}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=17cm]{FigS39-filtering.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Removal of potentially artifactual high-methylation reads}. Shown is unfiltered SMAC-seq data and the same locus after removal of all reads containing a 1-kb $\geq$75\% methylated stretch (average accessibility over 10 bp).
}
\label{FigS39}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS14-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrI;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrII.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS14}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS15-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrIII;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrIV.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS15}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS16-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrV;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrVI.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS16}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS17-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrVII;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrVIII.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS17}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS10-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrIX;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrX.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS10}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS11-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXI;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXII.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS11}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS12-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXIII;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXIV.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS12}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS13-centromeres-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule long-read accessibility around well positioned centromeres}.
(a) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXV;
(b) Raw unfiltered nanopore reads fully spanning the 4-kilobase neighborhood of the centromere of \textit{S. cerevisiae} chrXVI.
In both cases, accessibility is shown at aggregated 10-bp resolution (see Methods section for details) for the single-molecule display, and aggregated over sliding (every 5 bases) 50-bp windows for the genome browser SMAC-seq track.
}
\label{FigS13}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS54-multimappers_fractions.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Ribosomal DNA arrays are highly enriched for chromatin accessibility as measured by ATAC-seq}.
(a) Unique read mappability of the \textit{Saccharomyces cerevisiae} genome as a function of read length
(b) Unique read mappability of the \textit{Schizosaccharomyces pombe} genome as a function of read length
(c and d) Enrichment of multimapping reads in \textit{Saccharomyces cerevisiae} ATAC-seq datasets
(e) ATAC-seq multireads are highly enriched for rDNA-mapping reads 
(f) Enrichment of multimapping reads in \textit{Schizosaccharomyces pombe} ATAC-seq datasets
ATAC-seq datasets were obtained from Schep et al. 2015\cite{Schep2015} 
}
\label{FigS54}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=16cm]{FigS56-rRNA-diamide-30min-rep1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq reveals the distribution of alternative chromatin states of rDNA arrays. Shown are all reads covering the \textit{RDN37-1} and \textit{RDN37-2} arrays in the \textit{RDN1} locus in the ``diamide 30 min rep1'' experiment (unfiltered reads, ``aggregate'' signal)}.
}
\label{FigS56}
\end{figure*}

\clearpage

\begin{FPfigure}
\begin{center}
\includegraphics[width=16cm]{FigS57-rRNA-diamide-60min-rep1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq reveals the distribution of alternative chromatin states of rDNA arrays. Shown are all reads covering the \textit{RDN37-1} and \textit{RDN37-2} arrays in the \textit{RDN1} locus in the ``diamide 60 min rep1'' experiment (unfiltered reads, ``aggregate'' signal)}.
}
\label{FigS57}
\end{FPfigure}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=16cm]{FigS55-rRNA-diamide-0min-rep1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf SMAC-seq reveals the distribution of alternative chromatin states of rDNA arrays. Shown are all reads covering the \textit{RDN37-1} and \textit{RDN37-2} arrays in the \textit{RDN1} locus in the ``diamide 0 min rep1'' experiment (unfiltered reads, ``aggregate'' signal)}.
}
\label{FigS55}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=16cm]{FigS80-RDN2-NMI.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf NMI profile for the \textit{RDN37-2} array, as in Figure \ref{Fig2part2}b}.
}
\label{FigS80}
\end{figure*}

% \begin{FPfigure}
% \begin{center}
% \includegraphics[width=15cm]{FigS37-nucleosome-asymmetry.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Strand asymmetry of accessibility over positioned nucleosomes}. Shown is the average methylation status (either a 1-bp average, or aggregated over 10bp for every position in the genome, as indicated in each panel) in the neighborhood of positioned nucleosomes (as defined by H4S47C chemical mapping): 
% (a) Average forward/plus and reverse/minus strand methylation status in the 500-bp neighborhood of positioned nucleosomes, m$6$A + CG + GC positions;
% (b) Aggregated forward/plus and reverse/minus strand methylation status in the 500-bp neighborhood of positioned nucleosomes, m$6$A + CG + GC positions;
% (c) Average forward/plus and reverse/minus strand methylation status in the 100-bp neighborhood of positioned nucleosomes, m$6$A + CG + GC positions;
% (d) Aggregated forward/plus and reverse/minus strand methylation status in the 100-bp neighborhood of positioned nucleosomes, m$6$A + CG + GC positions;
% (e) Average forward/plus and reverse/minus strand methylation status in the 100-bp neighborhood of positioned nucleosomes, CG positions only;
% (f) Average forward/plus and reverse/minus strand methylation status in the 100-bp neighborhood of positioned nucleosomes, GC positions only;
% Asymmetry between the strand is only observed for m$6$A; it is not observed for the CG and GC sequence contexts.
% }
% \label{FigS37}
% \end{FPfigure}

% \clearpage

\begin{FPfigure}
\begin{center}
\includegraphics[width=17cm]{FigS53-TF-footpints-resolution.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf The impact of the addition of m$6$A to SMAC-seq on assay resolution and the potential ability to footprint individual transcription factors}. Shown is the fraction of motifs in the genome for each transcription factor in the yeast genome containing the indicated number of informative positions using GC alone,  m$6$A alone, and GC + m$6$A methyltransferase.
}
\label{FigS53}
\end{FPfigure}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS36-TF-footpints-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting by transcription factors}. Shown is the average methylation status (averaged over 10bp) in the neighborhood of occupied (as measured by ChIP-exo or ChIP-seq) recognition motifs for several \textit{S. cerevisiae} DNA binding proteins: 
(a) Reb1;
(b) Rap1; 
(c) Abf1;
(d) Cbf1. 
(e) ORC1. 
Strong footprinting is observed for Reb1, Rap1, and ORC1, while Abf1 and Cbf1 occupancy does not appear to be strongly protective against methylation. 
}
\label{FigS36}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS18-Reb1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting by Reb1 binding sites}.
(a) Raw unfiltered nanopore reads fully spanning the 400-bp neighborhood of a Reb1 binding site on chrXVIII, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS18}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS19-Reb1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting by Reb1 binding sites}.
(a) Raw unfiltered nanopore reads fully spanning the 400-bp neighborhood of a Reb1 binding site on chrXVIII, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS19}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS20-Rap1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with Rap1 occupancy}.
(a) Raw unfiltered nanopore reads fully spanning a 400-bp neighborhood of the subtelomeric region of chrXIV, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS20}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS21-Rap1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with Rap1 occupancy}.
(a) Raw unfiltered nanopore reads fully spanning a 400-bp neighborhood of the subtelomeric region of chrXV, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS21}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS22-Rap1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with Rap1 occupancy}.
(a) Raw unfiltered nanopore reads fully spanning a 400-bp neighborhood of the subtelomeric region of chrVI, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS22}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS23-Rap1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with Rap1 occupancy}.
(a) Raw unfiltered nanopore reads fully spanning a 400-bp neighborhood of the subtelomeric region of chrXII, at single-bp resolution;
(b) Same as in (a), but at aggregated 10-bp resolution.
}
\label{FigS23}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=9cm]{FigS64-ORC-footprints-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with ORC occupancy}.
(a) Raw unfiltered nanopore reads fully spanning the neighborhood of an ARS site on chrII, at 5-bp aggregated resolution;
}
\label{FigS64}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=9cm]{FigS65-ORC-footprints-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with ORC occupancy}.
(a) Raw unfiltered nanopore reads fully spanning the neighborhood of an ARS site on chrII, at 5-bp aggregated resolution;
}
\label{FigS65}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=9cm]{FigS66-ORC-footprints-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Single-molecule footprinting associated with ORC occupancy}.
(a) Raw unfiltered nanopore reads fully spanning the neighborhood of an ARS site on chrII, at 5-bp aggregated resolution;
}
\label{FigS66}
\end{figure*}

\clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.75\linewidth}
% \includegraphics[width=16cm]{FigS24-rRNA.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.25\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Single-molecule chromatin accessibility landscape of \textit{S. cerevisiae} rRNA arrays}. Shown is unfiltered SMAC-seq data around the first rRNA array, ``RDN37-1''.
%}
% \label{FigS24}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \begin{minipage}[c]{0.75\linewidth}
% \includegraphics[width=18.5cm]{FigS25-rRNA.png}
% \end{minipage}\hfill
% \begin{minipage}[c]{0.25\linewidth}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Single-molecule chromatin accessibility landscape of \textit{S. cerevisiae} rRNA arrays}. Shown is unfiltered SMAC-seq data around the second rRNA array, ``RDN37-2''.
% }
% \label{FigS25}
% \end{minipage}
% \end{center}
% \end{figure*}

% \clearpage


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS43-switch1.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts. 
% }
% \label{FigS43}
% \end{figure*}

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS45-switch3.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS45}
% \end{figure*}

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=13.5cm]{FigS44-switch2.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed in the vicinity of a repetitive element on chromosome IV}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS44}
% \end{figure*}

\clearpage

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS46-switch4.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS46}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS47-switch5.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS47}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS48-switch6.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS48}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS49-switch7.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS49}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS50-switch8.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS50}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS51-switch9.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS51}
% \end{figure*}

% \clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS52-switch10.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Sharp discontinuities in methylation signal observed nearby telomeric regions of some yeast chromosomes}. As these patterns are not observed equally on each strand, they are currently best explained as alignment and base calling artifacts.
% }
% \label{FigS52}
% \end{figure*}

% \clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=13cm]{FigS63-metanucleosome.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Metanucleosome NMI profiles in the yeast genome}.
Shown are average NMI maps between all 20-bp segments centered on each positioned nucleosome in the genome (a), the top 10\% strongly positioned nucleosomes (b), or the top 10\% strongly positioned nucleosomes (c).
}
\label{FigS63}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS67-gene-ends-correlation.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Patterns of coaccessibility between the 5$^\prime$ and 3$^\prime$ ends of genes}.
Shown is the average NMI for the $\pm$500bp regions in the 5$^\prime$ and 3$^\prime$ end of all yeast genes as well as the top and bottom 20\% expression-ranked genes (calculated over 10-bp windows). Only genes $\geq$1000 bp in length are shown. Similar results are obtained using windows of size 20bp or 50bp (data not shown).
}
\label{FigS67}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS87-TSS-correlation-NMI.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Accessibility correlation between TSSs in the yeast genome}.
Shown are NMI values for each pair of significantly and non-significantly correlated TSSs (defined as the regions $\pm$100 bp around the TSS).
}
\label{FigS87}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS88-TSS-correlation-chrII-265000-280000.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Example of accessibility correlation between TSSs}.
Shown are accessibility correlations between all loci in the region with coordinates ``chrII:265000-280000''. The size and color of the dots corresponds to their NMI scores.
}
\label{FigS88}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS86-TSS-correlation-MicroC.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Accessibility correlation between TSSs and 3D interactions}.
Shown are significantly and non-significantly correlated TSSs (defined as the regions $\pm$100 bp around the TSS) split into distance bins and the number of 3D interactions between each group (measured by MicroC).
}
\label{FigS86}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=9cm]{FigS58-diamide-RNA-seq.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Gene expression changes upon diamide treatment}.
Shown is RNA-seq data (mean and unit-variance normalized across time points) for all genes expressed at $\geq$50 FPKM in at least one time point.
}
\label{FigS58}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=9cm]{FigS68-diamide-HSF1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Dynamic changes in HSF1 occupancy upon diamide treatment}.
Shown are ChIP-seq RPMs (mean and unit-variance normalized across time points) for all HSF1 peaks detected in least one time point.
}
\label{FigS68}
\end{figure*}

\clearpage

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=14cm]{FigS89-HSF1-footprint.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Example of single-molecule footprinting by HSF1 at 30 minutes after diamide treatment}.
% }
% \label{FigS89}
% \end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS69-TMA10-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{TMA10} gene. 
% NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{TMA10} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS69}
\end{figure*}

\clearpage

\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS70-AAD6-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{AAD6} gene. 
% NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{AAD6} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS70}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS71-HOR7-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HOR7} gene. 
% NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HOR7} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS71}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS72-HSP26-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HSP26} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HSP26} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS72}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS73-HSP12-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HSP12} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HSP12} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS73}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS74-HSP31-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HSP31} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HSP31} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS74}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS75-HSP82-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HSP82} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HSP82} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS75}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS76-HSP104-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{HSP104} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{HSP104} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS76}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS77-DCS1-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{DSC1} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{DCS1} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS77}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS78-CTT1-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{CTT1} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{CTT1} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS78}
\end{figure*}

\clearpage
\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS79-SSA1-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{SSA1} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{SSA1} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS79}
\end{figure*}

\clearpage

\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS83-SSA4-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{SSA4} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{SSA4} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS83}
\end{figure*}

\clearpage

\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS84-GRE2-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{GRE2} gene. 
 % NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{GRE2} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS84}
\end{figure*}

\clearpage

\begin{figure*}
\begin{center}
\includegraphics[width=16.5cm]{FigS85-SRX1-V2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Coordinated changes in chromatin accessibility and nucleosomal occupancy during the yeast stress response}. Shown are changes in RNA Polymerase and HSF1 occupancy (measured by ChIP-seq), SMAC-seq profiles (1-bp resolution, 10-bp aggregate scores) and NMI profiles in the vicinity of the \textit{SRX1} gene. 
% NMI profiles and nucleosome occupancy ``states'' in the vicinity of the \textit{SRX1} gene. Nucleosome positions were called as occupied or unoccupied in $\pm$60 bp windows around the dyad positions determined by H4S47C chemical mapping.
}
\label{FigS85}
\end{figure*}

\end{document}
