Etienne Routhier

Edgard Pierre

Ghazaleh Khodabandelou

Julien Mozziconacci

Genome-wide prediction of DNA mutation effect on nucleosome positions for yeast synthetic genomics

Keywords: deep learning, genomics, S. cerevisiae, mutation, synthetic biology, nucleosome, DNA motifs

Introduction

The first genetically modified organisms were created in the seventies, shortly after Cohen et al developed the DNA recombination technology [START_REF] Cohen | Construction of biologically functional bacterial plasmids in vitro[END_REF]. This has been the foundation of biotechnology which is now a flourishing domain both in fundamental research and industrial applications [START_REF] Russo | Learning how to manipulate DNA's double helix has fuelled job growth in biotechnology during the past 50 years[END_REF]. The recent development of game changing technologies such as CRISPR and DNA oligonucleotide de novo synthesis now open the way to # Contact: julien.mozziconacci@mnhn.fr 1 major genome re-writing projects [START_REF] Ostrov | Technological challenges and milestones for writing genomes[END_REF]). A first paradigmatic example of this effort, the S. cerevisiae 2.0 (Sc 2.0) project [START_REF] Richardson | Design of a synthetic yeast genome[END_REF] will soon deliver the first example of a complete synthetic eukaryotic genome. Several projects are now starting with the aim to design more synthetic genomes [START_REF] Ostrov | Technological challenges and milestones for writing genomes[END_REF]) that could reach even the scale of the human genome. In the field of genomic engineering, the first step is to design the DNA sequence of interest, either resulting from very few edits of the wild type sequence, or from a more extensive genome re-writing, or even from the introduction of DNA sequences coming from a different organism. When introduced in the cell, this sequence will be interpreted by the cellular machinery and the resulting activity can be unpredictable. To date, there is no way to know whether the nucleosomes will assemble and position themselves on the DNA as expected, whether they will be modified or not by enzymes, or whether the chromatin will fold in space in an appropriate way. Since experimentally testing a huge quantity of trial sequences is cumbersome, if not unfeasible, computational tools are a good alternative to optimize the design of synthetic sequences so that they can fold into a functional chromatin in vivo. While this is a complex problem to deal with, the solution could come from the recent uptake of deep neural networks.

In parallel to the evolution of experimental genome editing techniques, the explosion of the amount of data available together with algorithmic advances and the use of graphical processing units (GPUs) [START_REF] Shi | Benchmarking state-of-the-art deep learning software tools[END_REF]) enabled the development of deep neural networks in many different contexts. This led to several breakthroughs in domains such as computer vision [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF], speech recognition [START_REF] Hannun | Deep speech: Scaling up end-to-end speech recognition[END_REF]) and machine translation [START_REF] Wu | Google's neural machine translation system: Bridging the gap between human and machine translation[END_REF]. As a data driven domain, genomics followed the trend and pioneering studies [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF][START_REF] Zhou | Predicting effects of noncoding variants with deep learning-based sequence model[END_REF] have demonstrated the efficiency of deep neural networks to annotate the genome with functional marks directly by interpreting the DNA sequence. The application of deep neural networks to genomics is growing at a high pace and it can now be considered as a state of the art computational approach to predict genomics annotations [START_REF] Kelley | Basset: learning the regulatory code of the accessible genome with deep convolutional neural networks[END_REF][START_REF] Umarov | Recognition of prokaryotic and eukaryotic promoters using convolutional deep learning neural networks[END_REF][START_REF] Min | DeepEnhancer: Predicting enhancers by convolutional neural networks[END_REF][START_REF] Eraslan | Deep learning: new computational modelling techniques for genomics[END_REF][START_REF] Quang | DANN: a deep learning approach for annotating the pathogenicity of genetic variants[END_REF][START_REF] Kim | Opening up the blackbox: an interpretable deep neural network-based classifier for cell-type specific enhancer predictions[END_REF][START_REF] Zou | A primer on deep learning in genomics[END_REF][START_REF] Jones | Computational biology: deep learning[END_REF][START_REF] Kelley | Sequential regulatory activity prediction across chromosomes with convolutional neural networks[END_REF]. One of the advantages of deep neural networks is their ability to predict a learned annotation on a variation of the genome, i.e. to predict the effect of mutations.

In this work, we report the use of deep learning to estimate the effect on nucleosome positions of changing each single nucleotide of the S. cerevisiae genome into another nucleotide.

The nucleosome positioning in S. cerevisiae has been extensively studied in the past by MNaseseq. This protocol relies on the enzymatic digestion of the linker DNA between nucleosomes and the sequencing of the protected DNA [START_REF] Krietenstein | Genomic nucleosome organization reconstituted with pure proteins[END_REF][START_REF] Hughes | Comparative genomics reveals Chd1 as a determinant of nucleosome spacing in vivo[END_REF][START_REF] Zhang | A packing mechanism for nucleosome organization reconstituted across a eukaryotic genome[END_REF]. Several studies pointed towards a close link between gene regulation and nucleosome positions [START_REF] Tsankov | The role of nucleosome positioning in the evolution of gene regulation[END_REF][START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF]. The role of the DNA sequence in the nucleosome positioning process has been a long standing debate. To assess this question, a pioneering study [START_REF] Segal | A genomic code for nucleosome positioning[END_REF] showed evidence for the existence of motifs negatively correlated with nucleosome positions [START_REF] Widom | Poly (dA-dT) promoter elements increase the equilibrium accessibility of nucleosomal DNA target sites[END_REF][START_REF] Iyer | Poly (dA: dT), a ubiquitous promoter element that stimulates transcription via its intrinsic DNA structure[END_REF][START_REF] Raisner | Histone variant H2A. Z marks the 5 ends of both active and inactive genes in euchromatin[END_REF]. Kaplan et al. [START_REF] Kaplan | The DNA-encoded nucleosome organization of a eukaryotic genome[END_REF]) developed a statistical method to predict the nucleosome density from the DNA sequence, emphasizing the preferential positioning of nucleosomes on specific DNA statistical motifs. Numerous computational methods -reviewed by Teif [START_REF] Teif | Nucleosome positioning: resources and tools online[END_REF] -were developed afterwards to predict the positions of nucleosomes from the DNA sequence. Recently, deep neural networks were also applied to discriminate between 147 bp long sequences bound by a nucleosome and 147 bp long sequences without any nucleosomes [START_REF] Gangi | Deep learning architectures for prediction of nucleosome positioning from sequences data[END_REF][START_REF] Zhang | LeNup: learning nucleosome positioning from DNA sequences with improved convolutional neural networks[END_REF].

Building on these previous works, we use here convolutional neural networks (CNNs) to predict the experimental nucleosome density (i.e. the results of the MNase protocol) for every position on the S. cerevisiae genome from the raw DNA sequence. The model reproduces well the characteristic nucleosome depletion around transcription start sites (TSS) as well as the typical periodic nucleosomal pattern on gene bodies. We then use the model as an in silico model of the yeast machinery to predict the effect of every single mutation along the genome. In doing so we assign to every nucleotide a score representing its importance regarding the nucleosome positioning process. This genomic track is accessible at https://github.com/etirouthier/NucleosomeDensity/Results nucleosome/mutasome scerevisiae.bw and can be used freely by others, when designing genetically modified yeast, to anticipate the effect of induced mutations on nucleosome positioning. We also use this track to analyze the DNA motifs that present a high mutation score, corresponding to motifs that are important for nucleosome positioning.

Results

Quality of the prediction

The first goal of this study is to accurately predict the nucleosome density directly from the DNA sequence. We use a CNN model whose input is defined by a one-hot-encoded DNA sequence of a given length L and whose output is the nucleosome local density associated with the nucleotide found in center of the input sequence. Several approaches aiming at extracting nucleosome positions from the nucleosome density have been proposed (e.g. [START_REF] Chen | DANPOS: dynamic analysis of nucleosome position and occupancy by sequencing[END_REF][START_REF] Chen | Improved nucleosomepositioning algorithm iNPS for accurate nucleosome positioning from sequencing data[END_REF] but our goal here is to predict the experimental output, that is the continuous nucleosome density, rather than the nucleosome positions that can be inferred from this experimental density. We present in this section a quantitative evaluation of the prediction quality.

A typical experimental result such as the one of [START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF], exhibits a locally periodic signal, with depleted regions preferentially found in intergenic regions (Fig. 1A, red signal). We train a CNN (refer to Methods for details) using the experimental nucleosome density of all chromosomes but the Chr16, which is kept aside as a test set. A length L = 2000bp of the input sequence was chosen and all sequences of length L obtained with a 1 bp sliding window on each chromosomes are used for training (Chr1 to Chr13) , validation (Chr14 and Chr15) and test (Chr16). The prediction on Chromosome 16 matches the experimental density, both in genic and inter-genic regions (Fig. 1A, blue signal).

A quantitative comparison between the two signals on the whole chromosome is displayed on Fig. 1B. Our method reaches a Pearson's correlation of 0.68 between prediction and experiment, a value comparable with the results obtained by state of the art CNN-based-methods on other tracks annotating the human genome, such as DNase sensitivity or histone modifications [START_REF] Kelley | Sequential regulatory activity prediction across chromosomes with convolutional neural networks[END_REF].

To investigate the generalisability of our results, we trained four CNNs models independently on four experimental replicates from a different data-set [START_REF] Kaplan | The DNA-encoded nucleosome organization of a eukaryotic genome[END_REF]). The Pearson's correlation between predictions (0.85 ± 0.05) is in the same range as the correlation between experimental replicates (0.87 ± 0.07). These values are constantly higher than the correlation between predictions and experiments (0.58 ± 0.05) which is itself lower that the correlation of 0.68 we obtained with the Hughes et al. data-set above (AL [START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF]) (Fig. 1C). An important control is that the performance obtained by comparing the predicted density with the experimental density coming from the dataset used for training the model is not significantly higher than the performance obtained by comparing the predicted density with the experimental densities from datasets which were not used for training (i.e. correlation between rep1 pred, as defined on Fig. 1C, and rep1 is comparable to the correlation of rep1 pred with the three other experimental densities). The model thus filters a large part of the replicates variability, which indicates a good generalisability of the results.

Training a network is a non deterministic process. We trained four networks on the same dataset and looked at the variability in predictions. We found that this variability can vary locally and is correlated with the experimental variability observed between replicates (Supplemental Fig S1). This indicates that using CNN could also be valuable to find regions of higher or lower confidence in experimental data.

To further investigate the generalisability of our models, we trained independently CNNs on experimental densities obtained under different experimental conditions [START_REF] Kaplan | The DNA-encoded nucleosome organization of a eukaryotic genome[END_REF].

Those conditions are the growth medium (YPD, YPEtOH, YPGal) and the presence or absence of a formaldehyde cross-linking step in the experimental protocol. For each condition several technical replicates are available. All the model-predicted densities do not correlate significantly better with the experimental densities that were obtained with the same growth medium as compared to the density that used to train the model. This points towards an overall similar nucleosome positioning in these different growth conditions (Fig. 1D). We next focus on regions surrounding the GAL1-10 promoters which are known to exhibit a different The models trained on experimental results including or lacking a cross-linking step produce quite different predictions when applied on the Chromosome 16 sequence. The models that were trained with an experimental density lacking a cross-linking step are predicting densities that correlate on average better (0.59) with experimental densities obtained with no cross-linking step than with experimental densities obtained using a cross-linking step (0.49).

Similarly, predicted densities obtained with a model trained on cross-linked data correlate on average better with experimental densities obtained with with cross-link (0.45 vs 0.39). The globally lower correlation values obtained using experiments that include a cross-linking step show that this step generates modifications in the experimental nucleosome density that cannot be predicted from the sequence alone. This suggests that this step can alter the nucleosome profile in a non-reproducible way.

For the sake of comparison of our CNN based method with previously proposed CNN based method for predicting nucleosome positions from DNA sequences, we used two previously proposed networks (J [START_REF] Zhang | LeNup: learning nucleosome positioning from DNA sequences with improved convolutional neural networks[END_REF][START_REF] Gangi | Deep learning architectures for prediction of nucleosome positioning from sequences data[END_REF] to predict the nucleosome density over the entire Chromosome 16. We found a correlation between prediction and experiment of 0.43 and O.40 to be compared with 0.68 obtained with our model trained and tested on Hughes [START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF] dataset. The lower performance of previously published methods for this task is expected since both methodologies were designed as classifiers which discriminates between fragments of DNA containing a nucleosome and fragments of DNA devoid of nucleosomes, whereas our method is designed to predict directly the nucleosome density over a whole chromosome. All methods nevertheless reproduce accurately inter-genic nucleosome depleted regions but our method reproduce with more accuracy the locally periodic density observed between depleted regions (see Supplemental Fig S3).

et al. (AL
Studying the effect of input length L on the predicted nucleosome phasing at TSS While specific DNA binding proteins usually recognize short DNA motifs, nucleosomes are positioned by a combination of several others mechanisms, including DNA local flexibility and shape, as well as the presence of neighboring nucleosomes [START_REF] Tsankov | The role of nucleosome positioning in the evolution of gene regulation[END_REF][START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF][START_REF] Zhang | A packing mechanism for nucleosome organization reconstituted across a eukaryotic genome[END_REF][START_REF] Riposo | Nucleosome positioning and nucleosome stacking: two faces of the same coin[END_REF][START_REF] Mavrich | A barrier nucleosome model for statistical positioning of nucleosomes throughout the yeast genome[END_REF]. The length L of the input sequence in our model is thus an important parameter that can change the performance of the CNN. We display on figure 2 a comparison between the predicted and experimental metagene profiles (i.e nucleosome density averaged over TSS regions) for different input lengths (L = 151, 501, 1001 and 2001bp). The experimental nucleosome density exhibits a characteristic pattern when averaged over TSS regions: a region with a low density, known as the nucleosome depleted region (NDR), precedes the TSS position. Further nucleosomes are regularly spaced with a periodicity of 167 bp on the gene body [START_REF] Riposo | Nucleosome positioning and nucleosome stacking: two faces of the same coin[END_REF][START_REF] Mavrich | A barrier nucleosome model for statistical positioning of nucleosomes throughout the yeast genome[END_REF]. This pattern reflects the molecular mechanisms at work in nucleosome positioning. Nucleosomes are excluded from regions preceding TSS, which are enriched with RNA polymerase. These regions act as barriers around which nucleosomes tend to stack upon each other, either due to thermal noise [START_REF] Mavrich | A barrier nucleosome model for statistical positioning of nucleosomes throughout the yeast genome[END_REF], or to chromatin remodelers The CNN model is able to identify and predict the NDR for all values of the input length L. This result is expected since short DNA motifs such as poly(A) motifs are known to exclude nucleosomes from these regions. However the prediction of the periodical pattern has a strong dependence on the input length L. For L = 151bp (which corresponds to a single nucleosome) the model is not able to recover the periodical pattern (Fig. 2A), the correlation between the predicted and the experimental metagene profiles is nevertheless high, it reaches 0.9 (Supplemental Fig S4). For L = 501bp (approximately 3 nucleosomes) a periodical pattern starts to appear and the first nucleosome after the TSS is well positioned (Fig. 2B), this improvement can be quantified by the correlation between metagene profiles which increases to 0.93. For L = 1001bp (6 nucleosomes) the periodical pattern improves and the first three nucleosomes after the TSS are well positioned (Fig. 2C), the correlation between the metagene profiles reaches 0.95. The best prediction quality is obtained for L = 2001bp, which corresponds to 12 nucleosomes, the typical size of longer genes. For this particular length of the input, the characteristic nucleosome pattern in TSS regions is well reproduced by our network (Fig. 2D), the correlation between the metagene profiles is accordingly increased to 0.97. 3 where the predicted density averaged over TSS is compared with two different experimental densities: (Fig. 3A)

the nucleosome density on the Chromosome F of k.lactis, (Fig. 3B) the nucleosome density on the Chromosome F of K.lactis transfected in S.cerevisiae. The CNN is able to capture accurately the NDR in K.lactis while being trained on S.cerevisiae. We conclude that DNA sequence motifs that determine NDR are similar between those two species. If we consider now the periodical pattern on the gene body, we can see that the prediction on K.lactis displays a periodical pattern but the value of the period, called nucleosomal repeat length (NRL), is the same as in S.cerevisiae (167 bp) whereas it should be 176 bp, as observed in K.lactis.

The predicted density is indeed similar to the experimental nucleosome density obtained on the transfected K.lactis sequences in S.cerevisiae. Our model is able to predict the behavior of the cell machinery of S.cerevisiae for the task of positioning nucleosomes on an exogenous genome.

Having carefully characterized the behavior of our model across replicates, experimental conditions and DNA sequences from different species, we now wish to use it to predict the effect of single mutations on the nucleosome positions.

Predicting the effect of single mutations

With our model in hand it is now possible to predict the nucleosome positions resulting from a mutation in the genome. The rationale behind this is that the more important a nucleotide is regarding nucleosome positioning, the more the effect of a mutation of this nucleotide will modify the predicted nucleosome density. In order to find positions on the genome associated with such modifications we generate all the possible single mutations along the genome and assign to every position a mutation score. This mutation score represents the Z-normalized distance between the nucleosome density predicted with and without mutation. Training several CNNs by letting aside each time different chromosomes, we computed the mutation score across the whole genome (see Methods for details).

A typical example of the mutation score along a region of Chromosome 16 with representative peaks at specific positions is outlined in Fig. 4A. Those peaks often coincide with NDR, represented as red dotted lines in Fig. 4A. Aligned and averaged around every NDR start, the mutation score displays a peak centered on the NDR (Fig. 4B). This result highlights the fundamental role of the NDR in nucleosome positioning.

The distribution of the mutation score (Fig. 4C) exhibit a narrow peak with over 90% of the values falling between -1 and 1 standard deviations, to be compared with 68% expected for a normal distribution. The distribution also features a long tail towards positive values corresponding to mutations having a strong impact on nucleosome positioning. We focus in the following on mutations with a score above 5, representing 0.6% of the genome.

To investigate the existence of DNA motifs for nucleosome positioning [START_REF] Segal | A genomic code for nucleosome positioning[END_REF], we analyse the motifs found at these high mutation score positions. We collect all 15-bp sequences surrounding a nucleotide with a high mutation score and extract from them overrepresented motifs (see Methods). Those motifs can be separated into two groups. The first group corresponds to poly(A) and the second group corresponds to poly(CG) (Fig. 4C).

The first group, poly(A), has previously been shown to be over-represented in NDR and to have a role in nucleosome exclusion from these regions [START_REF] Raisner | Histone variant H2A. Z marks the 5 ends of both active and inactive genes in euchromatin[END_REF][START_REF] Widom | Poly (dA-dT) promoter elements increase the equilibrium accessibility of nucleosomal DNA target sites[END_REF][START_REF] Segal | Poly (dA: dT) tracts: major determinants of nucleosome organization[END_REF][START_REF] Suter | Poly(dA•dT) sequences exist as rigid DNA structures in nucleosome-free yeast promoters in vivo[END_REF]. This effect has been proposed to be in part due to the natural stiffness of the poly(A) stretches [START_REF] Iyer | Poly (dA: dT), a ubiquitous promoter element that stimulates transcription via its intrinsic DNA structure[END_REF] The second group, poly(CG) also correspond to the binding sequence of a submit of the RSC (RSC3, [START_REF] Badis | A Library of Yeast Transcription Factor Motifs Reveals a Widespread Function for Rsc3 in Targeting Nucleosome Exclusion at Promoters[END_REF]). These sequences are found preferentially ∼ 100 bp upstream of TSS (see Supplemental Fig S5 ). Mutation of the RSC3 protein has been shown to result in an increase in nucleosome occupancy at NDR which contained a poly(CG) motif, suggesting that these sequences can as well exclude nucleosomes [START_REF] Badis | A Library of Yeast Transcription Factor Motifs Reveals a Widespread Function for Rsc3 in Targeting Nucleosome Exclusion at Promoters[END_REF].

The role of ploy(A) and poly(CG) have previously been described only in NDR upstream TSS and our findings are in line with these previous results. Indeed, only a fraction of those core motifs occurrences result in a significantly high mutation score (21% for poly(A), 13%

for poly(T) and 29% for poly(CG) Fig. 4D). These sites are enriched in NDR: when the motifs are present within gene bodies, their predicted impact on nucleosome positions is weaker.

We next ask whether all NDR have poly(A) or poly(CG) motifs and find that almost 60% of the NDR harbor at least one of these motifs (Fig. 4E). ∼ 40% of the NDR harbors motifs coming from one group only but the two groups are not mutually exclusive, since ∼ 20% of the NDR harbor motifs coming from both kind, i.e. one poly(A) and one poly(CG) (Fig. 4E ). We next investigate the relative position of these motifs relative to the TSS. Poly To investigate more quantitatively the effect of disrupting these motifs within the NDR we compute the averaged predicted nucleosome density in a 200 bp region centered on all motif instances within NDR with and without mutations in the motif. A mutation of a nucleotide in either a poly(A) or poly(CG) motif results in an increase of the nucleosome density in the vicinity of the mutation (Fig. 4F). A similar effect is seen for the complementary motifs poly(T) and poly(GC). This effect does not depend on the fact that one or two different motifs are found within a NDR. We conclude that, in agreement with previously reported experimental results, these two motifs are involved in the depletion of nucleosome. Using this methodology we do not find any motifs that would position nucleosomes by attracting them , i.e. motifs for which a mutation would locally reduce the nucleosome density.

Predicting the effect of multiple mutations

When designing synthetic genomes, one often needs to make several mutations in a given region. We therefore set out to investigate qualitatively the effect of having two or more mutations on the prediction of the nucleosome occupancy. We chose for illustration purposes a region of Chromosome 16 displaying two high mutation score positions in an inter-genic region (Fig. 5A, top). The two sites, numbered 1 and 2, fall into two NDR that flank a well positioned nucleosome. We computed the variation in predicted nucleosome occupancy that resulted in mutating each one of the sites, or mutating both sites. Mutation of site 1 results in a partial loss of the corresponding NDR. Nucleosome occupancy decreases at nucleosomal peaks in the vicinity of the mutation (indicated with dotted lines on Fig. 5A) and increases in linker regions. Mutation of site 2 induces an wider opening of the corresponding NDR as well as a sliding of the two neighboring nucleosomes away from the NDR. When mutating both sites, this results in a non-trivial combination of the two variations leading to an overall higher perturbation of the nucleosomal occupancy than for one mutation alone.

We then set out to quantify the average effect of having two mutations at high scoring mutation sites. We reasoned that the combination of two mutations may depend on the distance between these mutations. Based on the auto-correlation of the mutation score (Fig. 5B), we defined three different types of co-mutations: the first one for mutations that are closer than 5 bp corresponds to mutations in the same motif, the second one for mutations which are found between 5 and 90 bp away corresponds to mutations within a cluster of motifs, the third one corresponds to mutations which are found between 90 and 500 bp. Note that by construction of our network, mutations which are found at a distance greater than two input sequence length (here 4000 bp) will be independent so that the mutation score for the two mutations will be the sum of mutation score for each mutation taken independently. In order to evaluate the effect of two mutations and compare this effect with the effect of single mutations, we selected 1000 loci with high (> 5) mutation scores on each chromosome and computed their average mutation score (Fig. 5C, light grey). For all these mutations, we investigated the effect of high scoring mutations that were found within 5 bp (Short) by computing their average non-standardized mutation score (Fig. 5C, grey). We then compared these two values with the average non standardized score obtained by mutating both the primary and secondary mutations (Fig. 5C, red). We used here non-standardized scores which are always positive and additive by construction whereas standardized scores are not. This procedure was repeated for mutations found between 5 and 90 bp away (red Medium) as well as for mutations found between 90 and 500 bp away (red Long). We conclude that the effect of having two mutations of high mutation score within a region can be on average approximated by the sum of the effect of each mutation taken independently and that this effect does not depend on the distance between mutations. The sum of the mutation score of each mutation and the mutation score associated to their simultaneous mutations are strongly correlated (0.75 for short-, 0.70 for medium-and 0.80 for long-distance co-mutations, Fig. S6A,B,C). As a recommendation for genome design, we thus advice changing as few as possible nucleotides with high mutation scores.

When repeating a similar analysis for low mutation score nucleotides (score < 1) (Fig. 5D), we also found that mutating two nucleotides with low mutation score impacts on average the nucleosome occupancy in the proportions one would expect by adding the mutation scores for each mutation. The sum of the non-standardized mutation score of each mutations and the mutation score of the co-mutations are also strongly correlated (0.73 for short-, 0.82 for medium-and 0.97 for long-distance co-mutations, Fig. S6D,E,F). On average, mutating two such nucleotides leads to a standardized mutation score of 1, still much below the highest scores that can be obtained by changing some specific nucleotides with high mutation scores.

Nevertheless, when adding more and more mutations within the same region, the mutation score can reach values as high as 15, i.e. scores obtained when mutating two high score nucleotides (Fig. 5E). As a guideline for design, our analysis shows that on average 20 low score mutations have a similar effect as compared to one high score mutation.

These figures can serve as a baseline to evaluate the effect of a specific sequence design on 13 nucleosome positions, but in the case of a massive editing of many nucleotides, we recommend running a full prediction of nucleosome occupancy on the designed sequence in order to check for the potentially unwanted effects on nucleosome positioning. randomly sampled single nucleotides presenting a high mutation score (> 5, light grey), by mutating one by one all nucleotides that are found closer than 500 bp to these mutations (grey) and by mutating all pairs of nucleotides that are closer than 500 bp (red). Mutation scores are separated in 3 categories based on the distance between the two mutated nucleotides: less than 5 bp (Short), between 5 bp and 90 bp (Medium) and between 90 bp and 500bp (Long). (D) Same as (C) but for nucleotides with a low mutation score (< 1). (E) Evolution of the mutation score with the number of mutations with a low mutation score. All mutations considered here -taken individually -have a score smaller than 1. The solid line represents the average mutation score, the width of the line represents the standard deviation of the distribution of mutation scores.

Non

Discussion

In this study we used deep learning to generate a genomic track that leverages MNase-seq experimental results in order to predict the potential changes of nucleosome positioning resulting from mutating any bp in the S.cerevisiae genome. A similar procedure can in principle be done for any genomic track. The benefits of this track are two fold. Firstly, it can gives some guidelines to create a synthetic genome without modifying nucleosome positions in an unwanted manner. Secondly the results can be used to better understand how nucleosomes are positioned by the underlying DNA sequence.

The study of the predicted nucleosome phasing at TSS as a function of the input length of the network gives us more precise information. A network whose inputs are too short is able to capture the position of the NDR whereas is it not able to capture the periodicity in the pattern of nucleosome positions away from the NDR (Fig. 3). This observation leads to two conclusions : NDR are hard coded by DNA motifs whereas the DNA sequence wrapped around nucleosomes is not sufficient to precisely set their positions. In a cross species context, in which the network is trained in S. cerevisiae and the predictions are made on K.lactis genome, the periodicity is wrongly predicted to be the one of S. cerevisiae whereas the NDR are well predicted. These observations reveal the importance of the conserved DNA motifs residing in the NDR in the process of nucleosome exclusion upstream of the TSS. Studying the influence of single mutations all along the genome allows us to confirm this mechanism and to point out these specific motifs. Poly(A) and poly(CG) are, in agreement with earlier experimental studies [START_REF] Raisner | Histone variant H2A. Z marks the 5 ends of both active and inactive genes in euchromatin[END_REF][START_REF] Widom | Poly (dA-dT) promoter elements increase the equilibrium accessibility of nucleosomal DNA target sites[END_REF][START_REF] Suter | Poly(dA•dT) sequences exist as rigid DNA structures in nucleosome-free yeast promoters in vivo[END_REF][START_REF] Badis | A Library of Yeast Transcription Factor Motifs Reveals a Widespread Function for Rsc3 in Targeting Nucleosome Exclusion at Promoters[END_REF][START_REF] Krietenstein | Genomic nucleosome organization reconstituted with pure proteins[END_REF]), shown to be the core motifs preventing nucleosomes from binding in the NDR. While our study confirms the role of these core motifs, it also outlines that not all of these motifs in the genome are important for nucleosome positioning. We also show that other positions along the genome can play a role in this process. A general guideline for designing a synthetic yeast genome would be to preserve nucleotides that present a high mutation score.

For a quantitative anticipation of effects of mutation, the mutation score track is available at https://github.com/etirouthier/NucleosomeDensity/Results nucleosome/mutasome scerevisiae.bw.

Of course, the procedure used here in yeast for nucleosome positioning can be extended to other genomic tracks and other species. While we have validated here the high potential of this approach by studying how the DNA sequence drives nucleosome positioning in yeast, we anticipate that it will be very valuable tool to study nucleosome positioning rules in more complex organisms and ultimately in human. This would for instance empowers the community to ask whether some mutations frequently associated with diseases have a predominant role in positioning nucleosomes. Several issues will need to be solved to achieve this aim. The first is mappability, since many genomic regions are repetitive, the nucleosome density can not be measured on these sequences and this needs to be explicitly taken into account during training. The second is the size of the genome. The human genome is more than 200 times longer than the yeast genome. This has two impacts. The first is the coverage of the MNase experiment. The best coverage achieved for human cells is about 30 reads per nucleosome whereas it is usually 10 times higher in a standard yeast experiment [START_REF] Valouev | Determinants of nucleosome organization in primary human cells[END_REF]vs Kaplan et al. 2009). Last but not least, nucleosome spacing in human depends both on the cell type considered as well as on the location on the genome. The nucleosomal spacing is 167 bp in yeast, and the nucleosomal density shows only minor changes in different growth conditions (Fig S2). In human the spacing can change along the genome as well as in different cell types, taking values ranging from 178 to 205 bp [START_REF] Valouev | Determinants of nucleosome organization in primary human cells[END_REF]. We expect that these issues can be solved by using more sophisticated network architectures as well as increasing the computing power and that future developments in deep learning algorithms will become a game changing technology for genome writing.

Methods

Data accession and preprocessing

We use the reference genome sacCer3 of Saccharomyces cerevisiae, available at: http://hgdownload.soe.ucsc.edu/goldenPath/sacCer3/bigZips/.

MNase-Seq experimental results are available with accession numbers GSM552910 and GSE13622 (kaplan et al. [START_REF] Kaplan | The DNA-encoded nucleosome organization of a eukaryotic genome[END_REF]). Data for klyuverii Lactis used in the present study has been obtained from Tsankov et al. [START_REF] Tsankov | The role of nucleosome positioning in the evolution of gene regulation[END_REF]) and is available through accession number GSE21960. The experimental averaged nucleosome density in the TSS regions was obtained from the study of [START_REF] Hughes | A functional evolutionary approach to identify determinants of nucleosome positioning: a unifying model for establishing the genome-wide pattern[END_REF]) and is available using accession number GSM953761.

To obtain the nucleosome density from single end reads, we take the beginning of each read and add one count for each bp in a region of 100 bp in the direction of the read. We then truncate the obtained nucleosome density to a threshold -the 99 th percentile of the distribution of density scores -and divide the density by the threshold. This finally yields a density signal comprised between zero and one. We prepare input sequences of 151, 301, 501, 1001, 1501, 2001 bp for every positions in the genome except for Chromosome 16 to define the training (Chromosome 1 to 13) and validation (Chromosome 14 and 15) sets. We thus generate 10 613 042 input sequences among which those corresponding a nucleosomal density equal to zeros are excluded (as they correspond to non-uniquely mappable sequences). Each input sequence is then labeled with the nucleosome density value found at its central position.

Model Architecture and training

We implement the CNN using the Keras [START_REF] Chollet | Keras[END_REF] library and Tensorflow (Martin [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF] as back-end. A RTX 2080 Ti GPU is used to improve training speed. We use the adaptive moment estimation algorithm (Adam) to compute adaptive learning rates for each parameter and a batch size of 512.

Our CNN architecture (see Fig. 6A) consists of three convolutional layers with respectively 64, 16 and 8 kernels of shape (3 × 1 × 4), (8 × 1 × 64) and (80 × 1 × 16). The stride is equal to 1.

Our model takes inputs of shape (L, 1, 4), the last dimension representing the four nucleotides.

The first and second layer kernels identify 20 bp long motives which will play a role in the local affinity of the DNA sequence for nucleosomes. It is known for instance that poly(A)

will disfavor nucleosome formation while a ∼ 10bp periodic enrichment of AA/TT/TA dinucleotides that oscillate in phase with each other and out of phase with GC di-nucleotides will increase the affinity of the sequence for nucleosomes [START_REF] Segal | A genomic code for nucleosome positioning[END_REF]. The third layer is designed to capture long range information (coming from several nucleosomes) in order to grasp the nucleosome stacking one against another [START_REF] Riposo | Nucleosome positioning and nucleosome stacking: two faces of the same coin[END_REF].

The ReLU function is applied to the outputs of convolutional layers which are then entered into a max-pooling layer with pooling size (2 × 1). After each max-pooling layer batchnormalization is applied as well as a dropout with a ratio of 0.2. Finally, the output of the last layer is flattened and connected to the output layer containing one neuron through a single perceptron and using a linear activation function to make predictions. We tested several other architectures before choosing the one described above and a full recapitulation of hyper-parameters values which were tested is presented in Supplemental TableS1.

The loss function combines the Pearson's correlation (corr) between the prediction and the target and the mean absolute error (MAE) between them (loss = M AE(ŷ, y) + 1 -corr(ŷ, y),

with ŷ being the model prediction and y the target). The rationale for using this combined loss function is that we get a faster convergence and better final values both for the MAE and the correlation than using only one of them as a loss function (see Supplemental Fig S7).

An early stopping procedure is applied during training to prevent models from overfitting.

TSS alignment

For figure 2, genes positions of the studied species are download from the Ensembl fungi browser (ftp://ftp.ensemblgenomes.org/pub/fungi/release-46/gff3/saccharomyces_cerevisiae/ Saccharomyces_cerevisiae.R64-1-1.46.gff3.gz). The alignment on the TSS is simply made by taking a window of [-500, 1000] bp around every TSS and by averaging the signal. mutation score greater than 5 is considered to be the center of a 15-bp important motif. All those loci are collected (64610 loci) and aggregated when they intersect (23585 loci). We then use MEME [START_REF] Bailey | The MEME suite[END_REF] to extract significant motif logos from those loci. We used the following options: meme -oc outdir -nmotifs 10 -dna sacCer3peakseq.fa MEME is commonly used to extract binding sites logos from the DNA windows underlying peaks of ChIP-seq data, we use it to extract meaningful logos from the DNA windows containing nucleotides with high mutation scores.

  nucleosome occupancy profile in YPD vs YPGal (Supplemental Fig S2). Our prediction captures the nucleosome depletion at these gene promoters in YPGal (highlighted in light blue on Fig S2) but fails to reproduce the strong positioning of specific nucleosomes neighboring these promoters in YPD (highlighted in light grey on Fig S2). The model thus learns some of the the growth medium specific patterns at locations where the nucleosome density changes significantly between experimental protocols.

Figure 1 .

 1 Figure 1. Evaluation of the predicted nucleosome density. (A) Comparison between experimental (red) and predicted (blue) nucleosome densities in a region of the Chromosome 16. Genes are shown in blue on top of the two tracks. Data from AL Hughes, Jin, et al. 2012. (B) Density plot of the predicted nucleosome density in function of the experimental nucleosome density. The correlation between the two signals is 0.68. The distributions of the values of these two tracks on Chromosome 16 are also shown by side. (C) Cross-correlation between nucleosome densities on Chromosome 16 for 4 technical replicates (data from[START_REF] Kaplan | The DNA-encoded nucleosome organization of a eukaryotic genome[END_REF]) and 4 predictions obtained with models trained on each of the 4 replicates (e.g. rep1 pred is obtained with a model trained the rep1 nucleosome density). (D) Cross-correlation between nucleosome densities for 13 experiments and 13 predictions with models trained on each of the 13 experimental densities (experimental densities are on the horizontal-axis, predicted densities on the verticalaxis). The 13 experiments were carried out using different growth medium (namely YPD, YPEtOH and YPGal). Two different cross-linking conditions were used (cross-linking of nucleosomes on DNA prior to MNase digestion: CL or no cross-link: NOCL).

  Further increasing the input length does not change the performances significantly, while it penalizes the training process due an increasing need for computational memory. The global correlation between experimental and predicted densities over the whole chromosome increases from 0.63 to 0.68 (Supplemental Fig S7)when the input length L increases.

Figure 2 .

 2 Figure 2. Influence of the input length L on the pattern of the predicted nucleosome density in TSS region. Average predicted (blue) and experimental (red) nucleosome density in TSS regions. The predicted nucleosome density is here obtained with CNN models trained with different values of L ((A) : 151 bp, (B) : 501 bp, (C) : 1001 bp, (D) : 2001 bp). The other hyperparameters of the network are the same. Nucleosomes positions (in red for the experimental and blue for the predicted densities) are sketched below the curves. +1 and -1 refer to the first nucleosomes before and after the NDR.

Figure 3 .

 3 Figure 3. Predictions on the K.lactis genome compared with experiments. The prediction of the model, trained on S.cerevisiae and applied on the Chromosome F of K.lactis , is compared with the experimental data obtained by Hughes et al. (AL Hughes, Jin, et al. 2012). The signal around every TSS is aligned with respect to the first nucleosome downstream the TSS and averaged. (A) Endogenous context (Chr F of K.lactis in K.lactis) (B) Transfected context (Chr F of K.lactis in S.cerevisiae).

  and is enhanced and modulated by active nucleosome remodeling (Z[START_REF] Zhang | A packing mechanism for nucleosome organization reconstituted across a eukaryotic genome[END_REF][START_REF] Boer | Poly-dA: dT tracts form an in vivo nucleosomal turnstile[END_REF]). An important player in this process is the Remodeling the Structure of Chromatin Complex (RSC) which has been shown in vitro to clear promoters by removing nucleosomes from poly(A) sequences[START_REF] Krietenstein | Genomic nucleosome organization reconstituted with pure proteins[END_REF]).

  (A) and poly (CG) motifs are typically located 120 and 140 bp upstream of the TSS (Supplemental Fig S5 ). When both sites are present within a NDR, the poly(A) is on averaged moved further away from the TSS (160 bp). The position of the start site of the NDR does not depend on the group of motifs present: the NDR always starts on average 75 bp upstream of the TSS (Supplemental Fig S5 ).

Figure 4 .

 4 Figure 4. Effect of single mutations on the nucleosome density. (A) The mutation score on a region of Chromosome 16. NDR are shown with red dotted lines. (B) The average of the mutation score aligned on all NDR starts. On average the mutation score is peaking in the NDR, showing the major role of those regions in the nucleosome positioning process. (C) The distribution of the mutation score, as well as the 3 motifs enriched in the DNA sequences found in peaks with high mutation scores. (D) Proportion of poly(A),poly(T), and poly(CG) motifs in the genome that correspond with a high mutation score. (E) Proportion of four groups of NDR: NDR containing only poly(A/T) motifs (referred to as A-T), containing only poly(CG) like motifs (CG), containing both poly(A) and poly(CG) like motifs (A-T-CG) and NDR harboring none of these motifs. (F) Effect of a mutation in the poly(A/T) and poly(CG) motifs found in NDR on the nucleosome density in A-T and CG NDR respectively (top), and in A-T-CG NDR (bottom).

Figure 5 .

 5 Figure 5. Effect of mutating multiple nucleotides. (A) Illustration of the effect of two mutations on nucleosome occupancy on a locus of Chromosome 16. (top) mutation score over the region. (blue) predicted nucleosome occupancy for the wild type sequence. (light grey and grey) Local variation of the prediction obtained when mutating position 1 or 2 respectively. (red) Local variation of the prediction obtained when mutating both positions. Nucleosome occupancy peaks are highlighted with dashed lines.(B) Auto-correlation of the mutation score (semi-log-10 plot). Three different regimes can be identified, separated by dashed lines. (C) Average non-standardized mutation score obtained by mutating 16000 randomly sampled single nucleotides presenting a high mutation score (> 5, light grey), by mutating one by one all nucleotides that are found closer than 500 bp to these mutations (grey) and by mutating all pairs of nucleotides that are closer than 500 bp (red). Mutation scores are separated in 3 categories based on the distance between the two mutated nucleotides: less than 5 bp (Short), between 5 bp and 90 bp (Medium) and between 90 bp and 500bp(Long). (D) Same as (C) but for nucleotides with a low mutation score (< 1). (E) Evolution of the mutation score with the number of mutations with a low mutation score. All mutations considered here -taken individually -have a score smaller than 1. The solid line represents the average mutation score, the width of the line represents the standard deviation of the distribution of mutation scores.

Figure 6 .

 6 Figure 6. CNN architecture and mutation score computation. (A) The model is trained to predict the nucleosome density at the center position of a 2001 bp long DNA sequence. It contains three convolutional layers with maxpooling, batch-normalization and dropout. The final convolutional layer output is flattened and fed to a single output neuron. (B) We test all the possible mutations at the position indicated in blue, here a T, and predict the nucleosome density around this position with and without mutations. The mutation score is the Z-normalized sum of the distance between the wild type density and all the mutated type. The loss function used to train the network is used to compute the distance.

Figure 3

 3 Figure 3 displays the average nucleosome density in TSS region realigned on the first nucleosome downstream the TSS as previously performed by Hughes et al. (AL Hughes, Jin, et al. 2012).

Acknowledgements

We wish to thank Ayman Bin Kamruddin and Thomas Haschka for comments on the manuscript and Jean Baptiste Boule and Jean Baptiste Morlot for discussions. We also thank Michel Quaggetto for technical support. This work was supported by the IUF and the ANR-15-CE11-0023

Software availability

The mutation score track is available at GitHub (https://github.com/etirouthier/NucleosomeDensity/Results nucleosome/mutasome scerevisiae.bw) and as Supplemental Code. All of the code necessary to reproduce the results is accessible at GitHub (https://github.com/etirouthier/NucleosomeDensity.git) and as Supplemental Code

NDR determination

The NDR positions are defined as in [START_REF] Tsankov | The role of nucleosome positioning in the evolution of gene regulation[END_REF]. Briefly, the signal is firstly modified by setting to zero all the position with a value lower than 40% of the mean value, so that DNA linkers appear as a series of zeros. NDRs are then defined as the first linkers longer than 60 bp upstream of each TSS. If no sufficiently long linker is found closer than 1000 bp away from the TSS, the first zero is set as the beginning of the NDR.

Mutation score

To assign a mutation score to every position on the genome we use the methodology displayed on Fig. 6B. All the three possible mutations at a given position (highlighted in blue) are performed. The wild type and the three mutated genomes are used to predict the the nucleosome density. Predictions are made on the complete range in which the mutation can affect the model i.e. ±1000bp around the mutated position. Then -using the training loss functionwe compute the distance between the wild type local density and all the three mutated type local densities. By summing these three distances we assign a score to the mutated position.

This score is then Z-normalized within each chromosome to give the mutation score. This score reflects how the nucleosome density was perturbed by the mutation at the chosen position. Knowing that the nucleosome positions are not directly encoded in the underlying DNA sequence, we take long range perturbations into account using this methodology.

For the track presented along with the manuscript, the mutation score on Chromosome N is the average of three mutation scores obtained with three models independently trained on all the chromosomes with the exclusion of Chromosome N . While only a small fraction of the training set is sufficient to reach the maximum of accuracy (Supplemental Fig S8), we chose to use all sequences for training to improve the reproducibility of the mutation scores.

We finally obtain a robust mutation score, as two independently computed scores reach a correlation of 0.88 (Supplemental Fig S9). In this regard, it is good practise to train several models independently and to use the average of the prediction scores to assess the effect of a mutation.

Motif analysis

We make the assumption that a nucleotide assigned a high mutation score belongs to a motif that plays a role in the nucleosome positioning process. Every nucleotide assigned with a

Competing interest statement

The authors declare no competing interests.