Homec4science
Diffusion PhD_thesis (master) Browse

R9524/phd_thesis.tex
2cb4179d12c2master

Last Change
Raw File

phd_thesis.tex
View Options

\documentclass[12pt]{report}
\usepackage{graphicx}
\usepackage[round]{natbib}
\usepackage[utf8]{inputenc}
\usepackage[labelfont=bf]{caption}
\usepackage[margin={2.5cm,2.5cm}]{geometry}
\usepackage[titletoc,title]{appendix}
\usepackage{amsmath} % for equations with braces and text inside equations
\usepackage{hyperref} % for url
\setlength{\columnsep}{20pt}
\renewcommand{\familydefault}{\sfdefault}
% document
\begin{document}
\title{Bioinformatic analysis of transcription factor binding sites}
\author{Romain Groux}
\date{\today}
\maketitle
\chapter{Plan}
\begin{enumerate}
\item in vitro data study
\begin{itemize}
\item the analysis of SMiLE-seq data
\end{itemize}
\item in vivo data study
\begin{itemize}
\item SPar-K method
\item Chromatin architecture around TFBS
\item Motif co-occurence
\item Using single-cell ATAC-seq data to create a lexicon of chromatin architectures
\end{itemize}
\end{enumerate}
\chapter{Introduction}
TODO \citep{alipanahi_predicting_2015}.
\chapter{A study of chromatin accessibility in human blood monocytes}
Digital genomic footprinting (DGF) methods are a power powerful mean to reveal protein occupancy, genome-wide, at once \citep{vierstra_genomic_2016}. These methods allow to identify open chromatin regions within a genome and to list its active regulatory sites. These technologies are based on a targeted degradation of the open regions of the genome, either by DNaseI \citep{neph_expansive_2012} or by a transposon-based system \citep{adey_rapid_2010,buenrostro_transposition_2013}.
DGF techologies encounter a yet ever-growing popularity because of the wealth of data produced in a single experiment. Indeed, instead of running thousands of chromatin immonuprecitpitation followed by sequencing (ChIP-seq) \citep{barski_high-resolution_2007}, one per transcription factor (TF), to know where each TF is binding, it is sufficient to run a single chromatin accessibility assay. The price to pay for this gain of simplicity, compared to ChIP-seq, is a loss of information. Indeed, chromatin accessibility assays allow to list any active regulatory region within a sample, at once but do not give any information about which TF or complex is bound at a given location
The transposon-based method - named assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) - is currently gaining a rapid popularity due to its technical affordability and cheap labor costs compared to DNaseI-based methods.
\section{ATAC-seq}
\begin{figure}[!htbp]
\begin{center}
\includegraphics[scale=0.5]{scATAC-seq/figures/ATAC-seq2.png}
\captionof{figure}{\textbf{ATAC-seq principle :} ATAC-seq uses a hyperactive Tn5 transposase to simultaneously cleave genomic DNA at accessible loci and ligate adaptors to accessible DNA. These adaptors can serve as sequencing barcodes. A subsequent step of ligation allows to add sequencing adaptors. The released fragments are subjected to massively parallel sequencing of fragment ends (transposition sites) to generate a digital readout of per-nucleotide insertion genome-wide. A footprint is present at the level of the protected nucleotides. Taken and adapted from \citep{vierstra_genomic_2016}.}
\label{atac_seq_atac_seq}
\end{center}
\end{figure}
ATAC-seq assays uses a modified Tn5 transposase enzyme to selectively fragment the accessible regions of the genome \citep{adey_rapid_2010,buenrostro_transposition_2013}. The enzyme inserts small double stranded barcodes inside the DNA wherever it is accessible resulting a the creation of double strand breaks (shown in Figure \ref{atac_seq_atac_seq}). This process, known as tagmentation, allows to i) fragment the genome and ii) inserts sequencing primers at once. It should be noted that the Tn5 acts as an homodimer and thus inserts two copy of the same adaptors separated from each other by 9bp \citep{adey_rapid_2010}.
For a given genomic locus, the number of insertions depends on several parameters. First, the Tn5 transposition rate itself depends on the enzyme and substrate concentrations. Second, it has been demonstrated that the Tn5 posses a slight sequence preference \citep{adey_rapid_2010}. Finally, the most important parameter is steric hindrance. Stretches of DNA occupied by other proteins, such a transcription factors (TF) binding sites (BS) are protected from degradation. This usually leads to the creation of a characteristic signal at the level of the BS. If we represent the density of transposition events along a given region of the genome, the density usually show a decrease at the level of the BS. This decrease is named “footprint”.
The experimental readout of ATAC-seq is produced by sequencing the DNA fragment edges – the transposition sites - by using the inserted tags as sequencing primers. Mapping these sequences against the genome allows to retrieve the insertion sites. If the sequencing is pair-ended, it is possible to know the size of each fragment.
\section{Monitoring TF binding}
As discusses above, DGF assays are able to highlight active regulatory elements from an entire genome, at once. However, this come with the price of an information loss. First, even if we can identify active loci likely to be bound, we have no direct idea about the identities of the bound TFs. Second, we have no idea about the function of those regions. The function of a region – for instance transcriptional activation or repression - is ultimately bared by the TF and other complexes bound there. Thus delineating a region function necessitate to identify the TF and complexes bound here.
Retrieving the identify of the TFs/complexes bound require to implement dedicated strategies.
First, it is possible to collect evidences about the identity of TF likely to bind at a given location through a motif analysis. TFs can bind DNA directly through their own DNA binding domain or indirectly, through an interaction with at least on other partner TF which binds DNA directly \citep{neph_expansive_2012}. For a given TF, direct binding events can be detected by monitoring the presence of a binding motif if a specificity model is available. Thus a footprint baring a motif is likely to reflect a direct binding event. However, this method has two important limitations : related TF often share a common DNA specificity. For indirect binding, however, nothing can be done to detect such events. Also, evidences about the presence of biggest complexes can be collected by studying the size of the footprint. Large complexes should leave large footprints. This approach, even if limited is able to pinpoint a handful of candidate TFs.
Second, deciphering the functions of the regulatory elements can be undertaken by looking at the footprint produced by a given factor. Indeed, previous studies have showed that activator and repressor TFs tend to produce different types of footprints \citep{berest_quantification_2018}. Also, the spatial positioning of TF motif within the footprint seemed to be linked with the factor functions \citep{grossman_positional_2018}. For instance, factors associated with the regulation of transcription tend to have a motif in the middle of the footprint whereas factors known to interact with chromatin remodeling factors tend to have a footprint at the edge of the footprint, in contact with the surrounding nucleosomes.
\section{The advent of single cell DGF}
Recently, the advent of single-cell (sc) sequencing technologies have been a real game changer in the field of life science. These technological advances allowed to measure gene expression and chromatin accessibility (scATAC-seq) at a yet unprecedented resolution. As bulk sequencing was providing an average overview of what was going on, single-cell sequencing allows to monitor what is happening in each cell of a population. This advance had a profound impact on genomics for two reasons.
First, for the really first time, the heterogeneity of a cell population became accessible and could be studied at the chromatin, transcriptional and protein levels.
Second, the possibility of collecting high dimensionality data from tenth of thousands of individuals allows genomics to fully enter in the modern big data era, making commonly used machine learning methods usable as the number of parameters to estimate in the models became smaller than the number of individuals in the data \citep{angerer_single_2017}.
\section{A quick overview of scATAC-seq data analysis}
So far, most of the single cell technologies are targeted at measuring gene expression through scRNA-seq. Naturally, dedicated algorithms and computational methods have been developed to analyze these data. Currently, the most common types of analyses made are i) data projections and dimensionality reduction such as principal component analysis (PCA), t-stochastic distributed neighbours embedding (t-SNE) or uniform manifold approximation and projection (UMAP) and ii) cell population detection by clustering the cells based on the expression of genes \citep{fan_characterizing_2016, kiselev_sc3:_2017}, by reconstructing gene regulation network \citep{aibar_scenic:_2017} or by identifying cellular states based on the accessible region motif content \citep{gonzalez-blas_cistopic:_2019}. In all cases, the use of scATAC-seq data is to determined whether a region is accessible or not. The downstream analyses characterizes the accessible region using i) the number of reads mapping in these regions as a measure of the accessibility or ii) the sequence content within these accessible regions to determine regulatory topics.
\section{Open questions}
\begin{figure}[!htbp]
\begin{center}
\includegraphics[scale=0.5]{scATAC-seq/figures/pipeline.png}
\captionof{figure}{\textbf{framework to identify chromatin organization and use them to annotate cellular state :} the scATAC-seq data available in each individual cell are aggregated and used a if it was a bulk sequencing experiment. Regions of interest are listed using peak calling on the the bulk data. The read densities in these regions (center of the peaks +/- a given offset) are measured. The regions are then clustered based on their signal shape to identify different chromatin architectures and create a catalog. These chromatin signatures can then be used to annotate each region of interest in each cell, based on the signal resemblance. The information can be stored as a matrix (M) that can be used for downstream analyses, such as sub-population identification.}
\label{atac_seq_pipeline}
\end{center}
\end{figure}
All these methods have shown good performances to identify know and new cell populations [REFERENCES]. However, some issues remains open. First, none of these methods uses DGF data to identify different types of footprints or chromatin architecture, in terms of signal shape, at the single cell level. Second, ATAC-seq measures chromatin accessibility but also provides information about the nucleosome occupancy at accessible genomic regions \citep{buenrostro_transposition_2013}. Thus counting the number of reads mapping at a given loci is, indeed, an indication of accessibility but it does use only a small fraction of the available information. Finally, to date, no study has tried to determine whether what is observed at the bulk level can also be seen at the individual cell level and whether this can be used to infer the molecular state of the cells.
In this project, I designed and developed the basements of a computational framework to construct a catalog of prototypical chromatin architectures from single-cell data that can later on be used to annotate individual regions, in single cell. Such a method can be useful to determine cellular molecular state and to group cells accordingly. The entire pipeline is illustrate in Figure \ref{atac_seq_pipeline}.
\section{Data}
To this end, I choose to work with a publicly available single-cell ATAC-seq dataset from 5'000 human blood monocytes from a healthy donor. These data have been produced by 10xGenomics and are available as part of the advertisement for their single-cell library preparation and data analysis solutions.
For this reason, several data pre-processing steps have already been performed. Thus working with these data require minimum handling in order to be used. Cell demultiplexing has already been performed. Sequencing adapters sequences have been trimmed. Quality control checks have been performed and show good values. Mapping has been performed on the hg19 genome, peaks have been called and many more downstream analyses results and intermediate data are available.
Because of these reasons, this dataset offers all the conditions to be used as a standard to develop and benchmark new analyses methods.
The reads were downloaded in bam format from \url{http://s3-us-west-2.amazonaws.com/10x.files/samples/cell-atac/1.1.0/atac_v1_pbmc_5k/atac_v1_pbmc_5k_possorted_bam.bam} and the peaks called on the aggregated data were downloaded in bed format from \url{http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_pbmc_5k/atac_v1_pbmc_5k_peaks.bed}.
\section{Identification of catalog of chromatin architectures}
The study of signal shape (distribution) has been a quite active field for bulk sequencing experiments during the last decade. Dedicated algorithms \citep{hon_chromasig:_2008} \citep{nielsen_catchprofiles:_2012} \citep{kundaje_ubiquitous_2012} \citep{nair_probabilistic_2014} \citep{groux_spar-k:_2019} have been developed to cluster genomic regions based on their distribution of reads, which reflects their function. The major issue faced where i) to assess whether two regions had the same signal, they had to be properly aligned, ii) even if the regions were properly aligned, they had to be properly oriented and iii) the signal may be sparse due to an sub-optimal sequencing depth.
\subsection{EMRead : an algorithm to identify over-represented chromatin architecture}
\begin{figure}[!htbp]
\begin{center}
\includegraphics[scale=0.15]{scATAC-seq/figures/em.png}
\captionof{figure}{\textbf{Expectation-maximization algorithms :} illustration of the algorithms procedure. \textbf{A} illustration of EMRead, an algorithm dedicated to the discovery of over-represented chromatin patterns, as described in \citep{nair_probabilistic_2014}. \textbf{B} illustration of EMSequence, an algorithm to discover over-represented DNA motifs. The overall design is the same. Both algorithms model the data has having being sampled from a distribution and perform a maximum-likelihood estimation of the distribution parameters from the data through an iterative procedure.\\
EMJoint algorithm is the combination of both EMRead and EMSequence at the same time.}
\label{atac_seq_em}
\end{center}
\end{figure}
Most of the above mentioned algorithms and softwares deal with some of these issues. However, the probabilistic algorithm developed by \citep{nair_probabilistic_2014} is really interesting. For the sake of simplicity, let us call it EMRead. EMRead is a probabilistic partitioning method that softly clusters a sets of genomic regions represented as a vector of counts corresponding to the number of reads (ChIP-seq, DNase-seq) along them. The regions clustered based on their signal shape resemblance. To ensure proper comparisons between the regions, the algorithm allows to offset one region compare to the other to retrieve a similar signal at different offsets and to flip the signal orientation. Finally, it has been demonstrated to be really robust to sparse data.
This algorithm explicitly models the data has having being generated by sampling from $K$ distinct signal models that corresponds to $K$ different signal classes. Each class is represented by a vector of $L$ values that represented the aggregated signal for that class.
This algorithm models the data has having being generated, using a Poisson distribution, from $K$ signal models of lengths $L$. Each class is represented by a vector of $L$ values that represent the mean signal for that class. These values are then the mean parameters of the Poisson distributions
This algorithm models the signal over a region of length $L$ has having being sampled from a signal model, using $L$ independent Poisson distributions. The number of reads sequenced over this region is then the result of this sampling procedure. The entire set of regions is assumed to have been generated from $K$ different signal models (classes). Each class is represented by a vector of $L' \le L$ values that represent the expected number of reads at each position for that class. These values are thus the Poisson distribution parameters.
In order to discover the $K$ different chromatin signatures in the data, the algorithm proceed to a maximum likelihood estimation of the parameters using an expectation-maximization (EM) framework. Given a set of $K$ models, the likelihoods of each region given each class is computed. A posterior probability of each class given each region can, in turn, be computed. These probabilities can be understood as a soft clustering. The parameters of the classes are updated using a weighted aggregation of the signal. Since each region is computed a probability to belong to each class, it participates to the update of all the classes, with different weights.
If the length of the chromatin signature searched $L'<L$, then the algorithm slides a window along the regions and searched for this signature at each possible offset. This is how it deals with alignment issue. The signal orientation issue is tackled by also performing a searched with the flipped model. The procedure is depicted in Figure \ref{atac_seq_em}A.
\subsection{EMSequence : an algorithm to identify over-represented sequences}
This previous algorithm EMRead presented a really convenient point : it explicitly models the sequencing signal. Thus it can be adapted to search for different types of signals. EMRead is perfectly suited to discover different types of nucleosomes architectures or footprints. But because footprints reflects the binding of TFs, it is also critical to be able to identify the motifs within. To this end, I modified EMRead in order to discover over-represented sequence motifs. Let us called this new algorithm EMSequence. To this end, the following modifications have been applied to EMRead i) how the class signal is modeled and ii) the way data likelihood are computed and iii) the update of the class models. The process is illustrated in Figure \ref{atac_seq_em}B.
The input is composed of matrix $D$ of integer of dimensions $N \times L$ containing $N$ DNA sequences $d_{1}, d_{2}, ..., d_{N}$ of length $L$. Each sample $d_{i}=(d_{i1}, d_{i2}, ..., d_{il})$ is a vector of integers encoding the DNA sequence (A=1, C=2, G=3, T=4).
The $K$ classes profiles from which the data originate, instead of being modeled as signal profile, are modeled as sequence motifs $M_{1}, M_{2}, ..., M_{K}$ of expected base probabilities. A class motif $M_{j}$ is a matrix of dimensions $4 \times L'$ with the constrain $\sum_{i=1}^4 m_{i,j} = 1$.
\subsubsection{without shift and flip}
For the case where $L'=L$, the original equation (1) from \citep{nair_probabilistic_2014} to compute the probability of a sequence $d_{i}$ given a class $M_{j}$ is replaced by :
\begin{equation}
\begin{aligned}
P(d_{i}|m_{j}) & = \prod_{l=1}^{L} m_{j b,l} \\
\text{where } b & = d_{il}
\end{aligned}
\label{atac_seq_emseq_likelihood}
\end{equation}
Once the posterior probabilities $P(M_{j} | d_{j})$ have been computed, the original equation (3) in \citep{nair_probabilistic_2014}, to update the class models, is modified as follows :
\begin{equation}
\begin{aligned}
m^{*}_{j b,l} & = \frac{\sum_{i=1}^{N} (P(M_{j} | d_{il}) \times z}
{\sum_{k=1}^{N} (P(M_{j} | d_{il})} \\
\text{with } z & =
\begin{cases}
1, & \text{if $b = d_{il}$}.\\
0, & \text{otherwise}.
\end{cases}
\label{atac_seq_emseq_update_model}
\end{aligned}
\end{equation}
where $b$ takes the values $1,2,3,4$ for A, C, G and T respectively.
\subsubsection{with shift and flip}
For the case with shifting ($L'<L$) and flipping, the original equation (9) from \citep{nair_probabilistic_2014} to compute the probability of a sub-sequence of length $L'$ starting at offset $s$ in sequence $d_{i}$ given class $M_{j}$ is replaced by :
\begin{equation}
\begin{aligned}
P(d_{i} | M_{j}^{inv}, s) & = \prod_{l=1}^{L'} m_{j b,l}^{inv}\\
\text{with } b & =
\begin{cases}
d_{i,s+l-1} & \text{if $inv = 1$}.\\
4 - d_{i,s+l-1} + 1 & \text{if $inv = 2$}.
\end{cases} \\
\end{aligned}
\label{atac_seq_emseq_likelihood_shift_flip}
\end{equation}
where $inv$ is a notation indicating the orientation. If $inv=1$ we are searching in forward orientation (the forward strand) and $M_{j}^{1} = M_{j}$. If $inv=2$, we are searching in flipped orientation (reverse strand) and $M_{j}^{2}$ is the reverse complement motif of $M_{j}$.
Computing the reverse complement motif $M_{j}^{2}$ of a class motif $M_{j}$ is done using :
\begin{equation}
m_{j i,l}^{2} = m_{j 4-i+1, L'-l+1}
\label{atac_seq_emseq_reverse_motif}
\end{equation}
The computation of the posterior probabilities $P(M_{j}, s, inv | d_{j})$ remains the same as in \citep{nair_probabilistic_2014}. With the posterior probabilities, the model update can be undertaken. The original equation (12) in \citep{nair_probabilistic_2014} should be modified. The update of the model is made in a 2 step, by creating an intermediate motif for each strand separately and then by combining them, as follows :
\begin{equation}
\begin{aligned}
m^{*1}_{j b,l} = \sum_{s=1}^{S} \sum_{i=1}^{N} P(M_{j}, s, inv=1 | d_{il+s-1}) * z^{1} \\
\text{with } z^{1} & =
\begin{cases}
1, & \text{if $b = d_{il+s-1}$}.\\
0, & \text{otherwise}.
\end{cases} \\
m^{*2}_{j b,l} = \sum_{s=1}^{S} \sum_{i=1}^{N} P(M_{j}, s, inv=2 | d_{il+s-1}) * z^{2} \\
\text{with } z^{2} & =
\begin{cases}
1, & \text{if $4 - b + 1 = d_{il+s-1}$}.\\
0, & \text{otherwise}.
\end{cases} \\
m^{*}_{j b,l} = \frac {m^{*1}_{j b,l}} {\sum_{b'=1}^{4} m^{*1}_{j b',l}} +
\frac {m^{*2}_{j 4-b+1,L'-l+1}} {\sum_{b'=1}^{4} m^{*2}_{j 4-b'+1,L'-l+1}}
\label{atac_seq_emseq_update_model_shift_flip}
\end{aligned}
\end{equation}
where $m^{*1}$ is the partial motif for the forward strand, $m^{*2}$ is the partial motif for the reverse strand and $b$ takes the values $1,2,3,4$ for A, C, G and T respectively.
As in the original algorithm, the optimization process is then carried on for a given number of iterations.
\subsection{EMJoint : an algorithm to identify over-represented sequences and chromatin architectures}
Because EMRead and EMSequence algorithm computations are strictly identical with the exception of the likelihood computations and the model update, it is possible to design a third algorithm, let us call it EMJoint, that models at the same time the read coverage over a region and its sequence composition. To do so, I simply mixed both previous algorithms and applied to following modifications. For the sake of simplicity, I only expose the version with shift and flip as it is more general. This algorithm running is depicted in Figure \ref{atac_seq_em}.
The input is composed of two matrices of integers, $D$ and $R$, of dimensions $N \times L$. $N$ DNA sequences $d_{1}, d_{2}, ..., d_{N}$ of length $L$ and of $N$ vectors of read counts $r_{1}, r_{2}, ..., r_{N}$ of length $L$ are contained inside both matrices respectively. Each DNA sequence $d_{i}=(d_{i1}, d_{i2}, ..., d_{il})$ is a vector of integers encoding the DNA sequence (A=1, C=2, G=3, T=4) and each read count vector $r_{i}=(r_{i1}, r_{i2}, ..., r_{il})$ is a vector of integers containing the number of reads mapping over the sequences contained in $D$.
Each class is modeled by a vector of length $L'$ of expected number of reads $C_{j} = (c_{j1}, c_{j2}, ..., c_{jL'})$ and by a sequence motif $M_{j}$ of expected base probabilities $M_{j}$ of dimensions $4 \times L'$ with the constrain $\sum_{i=1}^4 m_{i,j} = 1$.
To compute the likelihood $P(r_{i}, d_{i}, s, inv |C^{inv}_{j}, M^{inv}_{j})$ of a region, equation \ref{atac_seq_emseq_likelihood} is modified as follows :
\begin{equation}
\begin{aligned}
P(r_{i}, d_{i}, s, inv |C^{inv}_{j}, M^{inv}_{j}) & =
\prod_{l=1}^L Poisson(r_{i,l}, \lambda=c^{inv}_{j,l})
\times
m^{inv}_{j b,l}\\
\text{with } b & =
\begin{cases}
d_{i,s+l-1} & \text{if $inv = 1$}.\\
4 - d_{i,s+l-1} + 1 & \text{if $inv = 2$}.
\end{cases} \\
\end{aligned}
\label{atac_seq_emjoint_likelihood}
\end{equation}
where $\lambda$ is the mean parameter of the $Poisson$ probability mass function.
The posterior probability $P(C_{j}, M_{j} | r_{j}, d_{j})$ computation remain unchanged. Once these values have been computed, it is possible to update both part $C_{j}$ and $M_{j}$ of a class using the original equation (11) from \citep{nair_probabilistic_2014} and equation \ref{atac_seq_emseq_update_model_shift_flip} respectively.
\subsection{Further generalizing EMJoint}
Based on the architecture of the EMJoint architecture, it is possible to design an more general algorithm that takes $Z$ different input matrices (layers) of dimensions $N \times L$ containing different types of signal (for instance DNA sequences, TF$_{1}$ ChIP-seq, TF$_{2}$ ChIP-seq, DNase-seq, ...) for a set of $N$ regions.
This only requires to adapt how the classes are modeled and equation \ref{atac_seq_emjoint_likelihood} to sum over the $Z$ different layers instead of only two. Additionally, care should be taken to use equation \ref{atac_seq_emseq_likelihood_shift_flip} for DNA sequence layer and equation (6) from \citep{nair_probabilistic_2014} for read count layers.
\section{Results}
Prior undertaking the chromatin organization study several pre-processing steps and checks have to been taken in order to ensure a proper treatment of the data.
\subsection{Open chromatin and nucleosome occupancy}
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{scATAC-seq/figures/fragment_lengths_classes.png}
\captionof{figure}{\textbf{Fragment size distribution : } the density of all the sequenced fragment was measured. The 3 peaks, from left to right, indicate i) the open chromatin fragments, ii) the mono-nucleosome fragments and iii) the di-nucleosome fragments. The 10bp oscillation reflect the DNA pitch.\\
A mixture model composed of 3 Gaussian distribution was fitted to the data in order to model the 3 classes. The classes fit are shown as dashed lines : open chromatin (red), mono-nucleosomes (blue) and di-nucleosomes (green). The violet dashed line show the sum of the three classes. The classes mean and standard deviation parameters were estimated using the Gauss-Newton algorithm.}
\label{atac_seq_fragment_size_dist}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{scATAC-seq/figures/fragment_lengths_class_prob.png}
\captionof{figure}{\textbf{Fragment size category probability mass functions : } probability that a fragment belongs to any of the three fragment classes, given its size : open chromatin (red), mono-nucleosomes (blue) and di-nucleosomes (green). The horizontal dashed lines indicates, for each class, the size limit at which the class probability drops below 0.9. The ranges spanned by each class are i) 30-84bp for open chromatin (red), ii) 133-266bp for mono-nucleosomes (blue) and iii) 341-500bp for di-nucleosomes (green). The upper limit of the di-nucleosome class was arbitrarily set to 500bp.}
\label{atac_seq_fragment_size_prob}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{scATAC-seq/figures/fragment_lengths_groups.png}
\captionof{figure}{\textbf{Final fragment size classes : } all fragments which sizes were overlapping the range spanned by a class, were assigned to that class. This ensured a high confidence assignment for more than 134 million fragments, leaving 46 millions of ambiguous and long fragments unassigned.}
\label{atac_seq_fragment_size_class}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{scATAC-seq/figures/ctcf_motifs_10e-6_aggregations.png}
\captionof{figure}{\textbf{Signal around CTCF motifs : } the human genome was scanned with a CTCF PWM and the densities of entire fragments (top row), entire reads (middle row), read edges or fragment centers (lower row) were computed around them. The densities were measured using bin sizes of 1 (left column), 2 (middle column) and 10bp (right column).}
\label{atac_seq_ctcf_all_data}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{scATAC-seq/figures/ctcf_sp1_myc_ebf1_footprint.png}
\captionof{figure}{\textbf{Signal around CTCF, SP1, myc and EBF1 motifs :} the human genome was scanned with using one PWM per TF. For each TF, the open chromatin architecture was measured by considering the corrected read edges (red) and the nucleosome occupancy by considering the center of the nucleosome fagments. The motif instances are located in the center.}
\label{atac_seq_ctcf_sp1_myc_ebf1_footprint}
\end{center}
\end{figure}
First of all, if a TF can protect a stretch of DNA against transposition and create a footprint, so can a nucleosome. As a matter of fact, both cases are biologically drastically different. Nucleosome compete with TFs to bind on DNA \citep{voss_dynamic_2014}. Thus nucleosome footprints represent regions of the genome that cannot be bound by TFs – if we except pioneering factors \citep{cirillo_opening_2002,zaret_pioneer_2011,soufi_pioneer_2015} – whereas TFs footprints actually represent TF binding events. Mixing both signals could drastically distort the analysis.
Nucleosomes fragments are expected to be large, as a nucleosome is usually wrapped around by 150bp of DNA whereas nucleosome free regions should be shorter. Long nucleosome free regions are unlikely because the longer an accessible fragment is, the most likely an insertion will happen. A fragment size analysis shows this clearly and allow to identify the different categories of fragments (Figure \ref{atac_seq_fragment_size_dist}). In this figure, open regions, mono-nucleosome anddi-nucleosome fragments are clearly visible. Morever, a 10bp periodicity oscillations reflecting the DNA pitch is also visible. This pattern is expected from [REFERENCE BUENROSTRO] and indicates a good data quality.
Rather than assigning arbitrary fragment size threshold to separate the categories, I preferred to use the data-driven approach developed by [REFERENCE Buenrostro]. The fragment sizes were fitted by a mixture of 3 Gaussian distributions (one per class, seen in dashed lines in Figure \ref{atac_seq_fragment_size_dist}) with their mean and standard deviation parameters learnt from the data. Then, the limits for each fragment class was defined as the size for which the probability of assignment to that fragment category was smaller 0.9 (Figure \ref{atac_seq_fragment_size_prob}).
Then each fragment was assigned to the class its size was overlapping. This method ensured the classification of 134 millions of fragments, leaving ~46 millions reads unassigned (Figure \ref{atac_seq_fragment_size_class}). However, this reduces drastically the risks of fragment mis-classification and protects the downstream analyses from a strong bias.
\section{Open chromatin and nucleosome occupancy signals}
We already discussed the importance of separating nucleosome fragments from nucleosome free fragments as they represent biologically different fragment categories.
The next question to ask was how should each category of fragment be represented along the genome?
First, for open chromatin fragment, it is clear that we want to know where the DNA is accessible. This information should be provided by the fragment edges – the transpositions sites. However, to account for the fact that the Tn5 transposase acts as a homo-dimer and inserts two barcodes side by side \citep{adey_rapid_2010}, the fragment edges positions were modified by +4bp for reads mapping the + strand and -5bp for reads mapping the - strand, as done in other studies \citep{buenrostro_transposition_2013,li_identification_2019}.
Second, for mono and di-nucleosome fragments, we are interested in knowing where the nucleosomes are sitting. For this, the fragment edges may not be the most informative. A better way to represent those fragments would be to use the center positions, which should correspond to the dyad for mono-nucleosomes [REFERENCE] or even to consider the entire reads or fragments.
To test these hypotheses I investigated the different signal aggregations predicted CTCF BS using a position weight matrix genome scan, using bin sizes of 1, 2 and 10bp. CTCF predicted binding sites were good candidates because CTCF is know to bind mostly through its motif [REFERENCE NEPH]. Additionally CTCF binding produces a really typical chromatin architecture with strongly positioned nucleosomes \citep{fu_insulator_2008} and leaves footprint \citep{neph_expansive_2012}.
As seen in Figure \ref{atac_seq_ctcf_all_data}, the open chromatin reads and fragments do not provide a signal fine enough to see a footprint (red lines). Both of them, nonetheless highlight open chromatin regions (red lines). The footprint becomes visible when considering the edges of the open chromatin fragments (bottom row). Increasing the bin size blurs it and eventually makes it disappear (10bp, lower right).
The nucleosomes can be revealed when using the centers of the mono-nucleosome fragments (bottom row, blue line).
Regarding the nucleosomes, considering the entire fragments blurs the signal (upper row, blue and green lines) and the entire reads reveal the region upstream and downstream of the nucleosomes (middle row, blue and green lines). The only way to obtain the precise location of the nucleosomes is to use the middle position of the mono-nucleosome fragments (bottom row, blue line). Interesting enough, the middle position of di-nucleosome fragments indicates the DNA linker between two adjacent nucleosomes but does not accumulate in open chromatin regions (bottom row, green line). This suggested that di-nucleosome fragments could be separated in two mono-nucleosome fragments. I tested this hypothesis by simply dividing a di-nucleosome fragment in two smaller ones, at its center position. I then pooled these new fragments with the mono-nucleosome fragments to create a nucleosome dataset. When looking at the middle of these fragments, they could perfectly reveal the nucleosomes directly adjacent to the CTCF motif. Additionally this nucleosome dataset was also able to reveal a second nucleosome in the arrays (bottom row, violet color). In the light of these results, I decided to use the edges of the open chromatin reads to investigate footprints and the fragment centers of the newly created nucleosome dataset to investigate nucleosome occupancy.
To further validate these choices, I also measured the chromatin organization around SP1, myc and EBF1 binding motifs as well. As shown in Figure \ref{atac_seq_ctcf_sp1_myc_ebf1_footprint}, the aggregation of the signal over the CTCF and SP1 motifs show an enhanced accessibility at the motif as well as a clear footprint. Moreover, the footprint is in a nucleosome free region. The situation is different for myc and EBF1. Both aggregations do not a clear nucleosome free region, neither a clear opening of the chromatin at the level of the motif. However, myc present an insert pattern that is compatible with a local protection of the DNA. On the contrary, EBF1 motifs shows an increased accessibility and a sligthly decreased accessibility on each sides.
CTCF and SP1 motifs are supporting the fact that footprints and nucleosome occupancy can be revealed using this method. Together with myc and EBF1, they clearly show an heterogeneity of chromatin organizations, at least at the aggregation level.
One possible cause for these differences could be from a differential TF activity within this cells. The most visible and obvious footprint could be the consequence of an increased TF activity, with more motif instances bound and leading to a clearer aggregation pattern. However, one should remain cautious on the interpretation of aggregation patterns as the individual sites may show several types of different chromatin organization that interfere with each other, creating an artificial aggregation that exist at none of the individual sites \citep{kundaje_ubiquitous_2012}.
\section{Setting up a method}
In order to
\chapter{Acknowledgements}
TODO
% bibliography
\bibliography{phd_thesis}
\bibliographystyle{abbrvnat}
\pagebreak
\begin{appendices}
TODO
\end{appendices}
\end{document}
c4science · Help